{'fraction_non_alphanumeric': 0.04079206979342478, 'fraction_numerical': 0.01366867140513942, 'mean_word_length': 4.3243893178078725, 'pattern_counts': {'":': 1, '<': 0, '<?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': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The mechanics of turbulence formation via an interaction potential: A simple two-particle model

In a recent work, we proposed the hypothesis that the turbulence in fluids could be produced by fluid particles interacting via a potential (for example, an interatomic potential at short ranges, and the electrostatic potential at long ranges). Here, we examine the basic mechanics of turbulence formation by studying a simple model which consists of only two particles, which interact via a potential. We start with the equations of motion for this pair of particles, then consider the corresponding Liouville equation for these dynamics, and subsequently derive the velocity moment transport equations. We close these equations using the hydrostatic and thermostatic approximations, and arrive at the closed equation for the velocity variable alone. We show that the velocity equation in the Fourier space can form time-periodic patterns under the potential forcing, coupled to a strong large scale flow. We find that the scaling of the equilibrium states of these patterns is consistent with the Kolmogorov energy scaling.


Introduction
In our recent work [1], we proposed the hypothesis that the turbulence in fluids could be produced by their particles interacting via a potential, which decays inversely proportionally to the distance between the particles. At short ranges, it could be the interatomic potential (such as the Thomas-Fermi [8,18], or the Ziegler-Biersack-Littmark [19] potentials), while at long ranges it could be the electrostatic potential.
The reason why we suggested such a radical hypothesis is the following. In [1], we investigated the full Liouville equation for N particles interacting via a generic potential, without reducing it to a single-particle equation as typically done in the conventional kinetic theory. We found that, due to the presence of the potential, the strong large scale flow creates the forcing in the 3-dimensional bundles of the full 3N-dimensional coordinate space, with each bundle belonging to an unordered pair of particles. We also found that these bundles are destroyed in the course of the standard Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) formalism [3,5,11], so that no such forcing manifests in the Boltzmann equation [4,7], and, subsequently, in the Euler [2] and Navier-Stokes [9] equations.
However, we noted that the direct observations and measurements of a turbulent fluid can register some bulk properties of the induced flow in these particle-pair bundles, such as the power scaling of the Fourier coefficients of the kinetic energy of the flow. In [1], we made crude estimates of the power scaling of the kinetic energy in the inertial range, Department  induced by a strong large scale flow via the interaction potential. These estimates suggested that the kinetic energy of such flow should scale in the Fourier space as k −3 at large scales (if the forcing was driven by the long-range electrostatic potential), and as k −5/3 at small scales, where k is the three-dimensional Fourier wavevector. Remarkably, we found our crude estimates to be consistent with direct measurements and observations [6,15]. The results in [1] further motivated us to look for a more detailed explanation of how turbulence could possibly be produced by the interaction potential. In the current work, we study the proposed mechanism of the turbulence formation in more detail; to simplify the setting, we consider only a single pair of particles which interact via the potential. The paper is organized as follows. In Section 2, we start with the equations of motion of the two particles, and change the variables to those which correspond to the motion of the center of mass of the pair of particles, and the differences between the two particles. The dynamics for these two variables are completely decoupled, which allows us to write the Liouville equation for the difference variables alone, and subsequently derive the corresponding velocity moment equations (just as done in the conventional kinetic theory for a single particle). We then combine the equations for the mass, momentum and energy transport to obtain a single equation for the velocity variable.
However, the equation for the velocity involves other variables, such as the density, temperature, shear stress and heat flux. In order to exclude these variables and obtain the closed equation for the velocity alone, in Section 3 we use the hydrostatic balance approximation for the transport of the momentum, and the thermostatic balance approximation for the transport of the energy. The hydrostatic balance approximation assumes that, for a nonstationary solution, the potential forcing term for the momentum transport remains balanced by the pressure gradient (the so-called hydrostatic equation). The thermostatic balance approximation assumes that, for a nonstationary solution, the emerging advective heat flux is automatically balanced by the conductive heat flux, so that the net heat flux of the system remains preserved. These two closures, along with some minor simplifications, result in a fully closed turbulent velocity equation.
To examine the scaling of the velocity in the Fourier space, in Section 4 we express the turbulent velocity equation in terms of the Fourier coefficients. This results in the infinite-dimensional system of nonlinear ordinary differential equations of the second order in time. We simplify this system by removing the majority of the coupling terms (assuming that their overall effect is akin to a random noise), and only leaving the coupling to the strong large scale flow, which is presumed to be the triggering mechanism for the turbulence formation. The latter simplification exposes what we presume to be the fundamental mechanism of turbulence formation -the cubic oscillator with a linear term, coupled to the large scale velocity forcing via the interaction potential.
In Section 5, we proceed to study the behavior of the discovered oscillator. Assuming that the strong large scale flow is a jet, we identify the following regimes of its behavior: a) The constant acceleration regime, where the wavevector of the chosen Fourier coefficient is precisely orthogonal to the direction of the strong large scale flow; b) The linear oscillator regime, where the Fourier wavevector is parallel to the strong large scale flow; c) The cubic oscillator regime, where the Fourier wavevector is nearly orthogonal to the strong large scale flow. In the latter case, we find that the Fourier coefficients of the velocity, whose direction is parallel to the large scale flow, develop periodic motions with large amplitude. This amplitude scales as the inverse two-thirds of the length of the Fourier wavevector.
In Section 6, we show that the corresponding time-windowed averages of the kinetic energy scale as the inverse five-thirds of the length of the Fourier wavevector. This coincides with the famous Kolmogorov scaling [12-14, 16, 17]. The results of this work are summarized in Section 7.
2. The transport equation for the turbulent velocity 2.1. The equations of motion. We start with the dynamical system consisting of two identical particles, with coordinates x 1 and x 2 , and velocities v 1 and v 2 , respectively. These particles interact with each other via the potential φ(x). Additionally, we assume that the particles may optionally be forced by an external acceleration g (such as the Earth gravity, for example). The corresponding equations of motion for the coordinates and velocities of these particles are given via Above, observe that the forcing terms in both equations depend only on the differences in coordinates of the two particles. Thus, it is convenient to change the variables of the dynamics so that the center of mass of the system is one of the variables, and the difference in coordinates is the other one. First, we observe that which allows us to write (2.3a) d dt Next, we note that the motion of the center of mass of the system is affected solely by the external acceleration. At the same time, the dynamics for the differences of coordinates and velocities are fully decoupled from the dynamics of the center of mass of the system, and are affected solely by the interaction potential. Thus, denoting x = x 2 − x 1 and v = v 2 − v 1 , we can consider the dynamics for x and v separately from the dynamics of the center of mass: Throughout the work, we will assume that φ(x) ∼ 1/x, and, on the scales of interest, its Fourier coefficient has the form where φ 0 is a real constant. These are the same assumptions as in [1]; in particular, the Thomas-Fermi [8,18] and Ziegler-Biersack-Littmark [19] interatomic potentials, as well as the electrostatic potential, have the Fourier coefficients of the form (2.5) in the range of interest (for details, see [1]).

The Liouville equation and the moment equations.
Let f (t, x, v) describe the distribution density of the differences of the coordinates and velocities, irrespectively of the distribution of the coordinate and velocity of the center of mass. The corresponding Liouville equation for this probability density f is given via Observe that (2.6) is, effectively, the single-partcie transport equation with a forcing term, even though x and v are not the actual coordinate and velocity, but rather the difference between the coordinates and velocities of the two particles. In what follows, we use the old-fashioned divergence notation for the trace of a vector field derivative whenever stylistically appropriate: We will also use the following notation for the Jacobian of a vector field: that is, in contrast to the divergence, the Jacobian is the "outer product", with a being the column vector, and the differentiation operator being the "row vector". These operations are easily extended onto matrices -for example, the divergence of a matrix is the contraction of the differentiation operator along each column of the matrix. As usual, we denote the velocity moments via the density ρ, average velocity u, and energy E, where v 2 = vv T is the outer product of v with itself. Then, for the moments of (2.6), we integrate by parts the forcing terms with v-derivatives, and obtain where the symmetric 3-rank tensor v 3 is the outer product of v with itself, computed twice.

The equation for the turbulent velocity.
According to our hypothesis in [1], the turbulent velocity motions arise when a strong large scale flow communicates its presence onto the inertial scales via the interaction potential φ. Thus, in order to study this process in the form of a dynamical system, it is best to have the equation which governs the evolution of the velocity, and, at the same time, includes the terms which couple the velocity itself and the interaction potential. The moment equations in (2.10) do not, however, satisfy this requirement -indeed, observe that the potential term in the momentum equation (2.10a) is not coupled to the velocity u; instead, such a coupling is present in the energy transport equation (2.10b). In order to combine the time-derivative of the velocity and its coupling to the potential into a single equation, we apply the following series of manipulations: 1. Differentiate the momentum equation in (2.10a) in t, where in the right-hand side we used the mass equation in (2.10a) to replace the time-derivative of ρ with the divergence of the momentum. 2. Differentiate the energy equation (2.10b) in x: where ":" denotes the "Frobenius product" of the double-differentiation operator with the two indices of a symmetric 3-rank tensor. 3. Combine the resulting equations using their common mixed derivative of the energy: 4. Expand the second time-derivative of the momentum as (2.14) and express the time derivatives of the density and velocity via the transport equations for the mass and momentum in (2.10a): The substitution of the above expressions into (2.13) leads to The transport equation above has the desired form -the only time derivative it contains is the one for the velocity (which happens to be of the second order), and it also has the forcing terms which couple the velocity to the potential. Hence, we refer to the equation in (2.16) as the turbulent velocity equation.

The closure for the turbulent velocity equation
Observe that the turbulent velocity equation in (2.16) is not closed with respect to its prognostic variable u, since it is also coupled to the density ρ and other velocity moments. In what follows, we introduce suitable approximations which allow us to discard, for the most part, the unwanted couplings.
Recall that the energy and the cubic moment can be expressed, respectively, via where T and Q are the kinetic temperature matrix and the skewness tensor, respectively: Above, the subscripts "T" and "TT" denote the two possible permutations of indices of a 3-rank tensor (just as for a matrix, where the transposition is a single permutation of its two indices). The quantity ρT is known as the pressure tensor. Additionally, we separate T into the sum of its own trace, and the remainder of the matrix S: Above, θ is the kinetic temperature, and S is the shear stress, whose trace is zero, by construction. The product ρθ is known as the pressure. In a typical fluid dynamics scenario, the standard closure of the moment equations is either the Euler closure (S and Q set to zero), or the Navier-Stokes closure (S and Q are expressed via the Newton law of viscosity, and the Fourier law of heat conduction, respectively). However, these laws are the consequences of the presence of the collision integral in the Boltzmann equation (for details, see, for example, [10]). This collision integral is time-irreversible, and, in effect, describes the irreversible dissipation of the kinetic properties of the flow into its heat properties (particularly, the Newton law of viscosity describes the momentum dissipation, while the Fourier law of heat conduction describes the energy dissipation). In the current setting, the irreversible collision integral is not present, and, thus, the standard closures of molecular kinetics cannot be directly applied. Therefore, we will need to close the turbulent velocity equation in (2.16) based on different practical considerations.
3.1. The closure for the shear stress. As the first step towards the closure, we are going to assume that the shear stress is, in practice, much smaller than the square of the velocity, that is, S ≪ u 2 . While it is true that our dynamical system lacks the Newton's law of viscosity to dissipate the shear stress, we will simply assume that, in a turbulent scenario, S does not become large enough to match the magnitude of u 2 (at least not on the time scale of interest). With this, we discard S from the expression for the energy and the cubic moment, which subsequently become 3.2. The hydrostatic balance approximation. The second approximation we introduce is known as the "hydrostatic balance" -we assume that, in the momentum equation in (2.10a), the gradient of the pressure ρθ is balanced by the potential forcing in the momentum equation: Observe that this relation is exact in the momentum equation in (2.10a) in the case of the steady state with u = 0, and all we do is extend it onto an unsteady flow near this steady state. This approximation is known to hold very well in practice for typical subsonic flows (many operational atmospheric models are, in fact, hydrostatic). Together with the approximation of the negligible shear stress, the hydrostatic balance approximation leads to and, therefore, the turbulent velocity equation in (2.16) becomes 3.3. The thermostatic balance approximation. In order to close the system with respect to the skewness tensor, we will use a similar principle as the hydrostatic balance above, except that it will be applied to the net heat flux of the system. The conductive heat flux q is defined via the following velocity-centered cubic moment: Observe that 2q is the contracted, along any pair of indices, skewness tensor Q: It is easy to see that the full skewness tensor Q can be written in the form where R is a fully symmetric 3-rank tensor whose contractions along any pair of indices are zero. The skewness deviator R is related to Q in the same way as the shear stress S is related to the kinetic temperature matrix T -namely, R quantifies the deviation of Q from its own trace along any pair of its indices. In the current work, we are going to make the same assumption for R as we made above for the shear stress S, that is, R ≪ u 3 , and thus it can be discarded from the cubic moment in (3.4): Next, we are going to introduce the appropriate closure with respect to the heat flux. Recall that, in practical real-world settings, the moment transport equations for the probability density of a single particle can be endowed with external heat fluxes. In particular, this can be seen in the atmospheric models, which exchange the heat energy with the surface of the Earth, absorb the energy arriving from the Sun, and also drain the heat energy into the outer space in the form of the electromagnetic emission. Thus, the dynamics of the atmosphere are closed with respect to the mass transport (as the mass is not exchanged with the outside medium), largely closed with respect to the momentum transport (as the momentum is not exchanged with the outside medium, except for the fictitious Coriolis force due to Earth rotation), but completely open with respect to the heat exchange. Thus, the balance of the heat fluxes is the important benchmark of any realistic model of the atmosphere.
Here, we are also going to assume that our system is thermally open, and thus must balance the heat fluxes. For simplicity, we will assume that the net external heat flux in the system is zero. Now, assume that, under the above condition, together with the hydrostatic balance approximation in (3.5), the system of moment equations has a steady solution with u = 0 (that is, the gas is at rest). If this is possible, it implies that q = 0, otherwise there would be a nonzero forcing term in the energy transport equation. Now, if a nonzero velocity u is introduced into the system, it creates the advective heat flux q a , given via where c p = 5/2 is the heat capacity at constant pressure. Observe that c p depends only on the number of degrees of freedom which can store the kinetic energy, and since, in our variables, the pair of particles has three degrees of freedom, its c p is the same as that of the usual ideal gas. The advective heat flux q a is the amount of the heat energy mechanically transported by the flow with the velocity u and temperature θ. In the current work, our assumption is that the heat flux remains balanced, which means that the system responds to the change in velocity via the "automatic" creation of the conductive heat flux q, which counters q a so that the net heat flux in the system remains preserved: Thermodynamically, it means that, even though the advective heat flux is expected to carry the heat energy by transporting the mean temperature along the mean flow, the distribution f "automatically" adjusts its skewness in the opposite direction, so that the net heat flux of the system remains preserved. We call the relation in (3.13) the thermostatic balance approximation. Upon the substitution of (3.12) and (3.13) into (3.11), the expression for the cubic moment simplifies to and the turbulent velocity equation in (2.16) subsequently becomes 3.4. The closed form of the turbulent velocity equation. At this point, what remains to be done is to exclude the density ρ from the turbulent velocity equation. For that, first observe that which leads to The left-hand side is now a multiple of ρ. Next, in the right-hand side, we express where ∆ is the Laplace operator (that is, the trace of the second derivative). We approximate the derivatives of ρ in the potential forcing term via the hydrostatic balance relation, assuming that the temperature θ is a constant (i.e. the system is near the thermodynamic equilibrium): This leads to Observe that now all terms of the turbulent velocity equation are multiples of ρ, which thus can be canceled out: For the final simplification, we will assume that, on the scales of interest, the potential is sufficiently weak compared to the kinetic temperature θ, and thus the last term in the right-hand side is small in comparison to the rest and can be discarded. This leads to the turbulent velocity equation in the following closed form: Observe that, in the present form, the closed turbulent velocity equation above in (3.22) is a nonlinear PDE, which was derived purely from the "physical principles", without taking into account its mathematical properties. As a result, it may have solutions which blow up in finite time, or otherwise behave irregularly. However, for the purpose of this work, we will presume that finite, physically meaningful solutions of (3.22) also exist, at least on the time scale of interest.

Simplified dynamics for the Fourier coefficients of turbulent velocity
In terms of the Fourier coefficientsû k , the closed turbulent velocity equation in (3.22) is given via where " * "denotes a convolution, and "⊛" is the combination of the convolution with an outer product. In what follows, we make the same assumptions as in [1], namely that, initially, the solution of (3.22) is placed near its "unexcited" state, with u = 0, perhaps with small random fluctuations (which are always present in practice). At the initial time, a strong jet u 0 is instantaneously created at large scales -that is, it is confined to a small subset of wavenumbers k ∼ 1. By "strength", here we refer to the average energy of u 0 , which we assume to be much greater than the energy φ 0 of the interaction potential in (2.5): where V is the physical volume of the domain, and the equality is due to Plancherel's theorem. This condition ensures that u 0 effectively obeys the moment equations in (2.10) without the potential forcing, and thus behaves in a "normal way" -that is, as a typical large scale flow in a fluid dynamics setting. Our wavenumber of interest k is located at the upper range, and corresponds to small scales, whose motions change rapidly in time. We assume that, on the characteristic time scale of k, the Fourier coefficientsû 0 k of the strong large scale flow can be regarded as constants. With this setup, we further make the following assumptions: a) Only the Fourier coefficientsû 0 k of the large scale flow can create substantial forcing at the inertial range via the potential in the right-hand side of (4.1), and b) Due to the fact that the initial conditions for the inertial range wavenumbers consist essentially of small random decorrelated fluctuations, the coupling ofû k to the surrounding coefficients via the convolutions in (4.1) can be regarded as chaotic and random forcing. Thus, in the terms with convolution, we assume that the only coupling that should be taken into account is • that which self-couplesû k to itself and/or its complex-conjugateû * k ; • that which couplesû k to the Fourier coefficientsû 0 k of the strong large scale flow. While the chaotic/random forcing from other inertial range velocity coefficients should generally have an effect on the dynamics ofû k , here our goal is to isolate the fundamental mechanics of turbulence formation from other effects, and discarding this random "noise" serves the stated purpose.
Under these assumptions, the two terms with the potential forcing in the right-hand side of (4.1) are given via where we observe that m ∼ 1, and take into account the form ofφ k in (2.5). Thus, leaving only the dominant term in the right-hand side of (4.1), we arrive at Next, it is straightforward to show that the terms with convolutions in (4.1) become (4.5a) 2(û k k T ) * (û k k T ) * û k ≈ −2(k ·û k ) 2û * k + For the sum of the two terms, we subsequently have Clearly, the second sum in the right-hand side above is much smaller than the first sum, because the latter is of order k 2 , while the former is of order k . Thus, we further drop the second sum from the approximation. Finally, we note that where V is the volume of the physical domain, and the last identity is a consequence of Plancherel's theorem. Denoting the bulk energy matrix of the large scale flow viaĒ 0 , we have (4.9) 2(û k k T ) * (û k k T ) * û k + (û k ⊗ k 2 ) * (û k ⊛û k ) ≈ k 2 :Ē 0 + 2|k ·û k | 2 û k − (k ·û k ) 2û * k . The simplified form of the turbulent velocity equation in (4.1) is thus given via (4.10) d 2û As the equation above is fully closed with respect to u k and does not involve coupling with the Fourier coefficients of other wavevectors, we henceforth drop the subscript "k" fromû k , and denote it asû.

Behavior of turbulent velocity in the simplified model
Recall that, under our prior assumptions, the initial stateû and its time derivative are zero (at the time when the strong large scale flow u 0 is introduced into the system). If so, initially the second time-derivative ofû will be driven exclusively via the forcing in the right-hand side of (4.10). From the structure of this forcing, we can see that it generally separates into two components: • One component is in the direction of k (that is, aligned with the wavevector ofû); • Another component is in the direction ofφ * û 0 (that is, aligned with u 0 , since the latter is a large scale jet). If these two directions are arbitrary and unrelated to each other, then the dynamics of u appear be rather complicated for an easy examination. Thus, here, for simplicity, we will restrain ourselves to few special cases, where the dynamics of (4.10) are relatively easy to examine.
Observe that the left-hand side of (4.10) is a combination of the linear and nonlinear (that is, cubic) oscillators. Thus, we will examine three standard patterns of behavior: a) Constant acceleration, where only the time-derivative ofû is driven by the potential forcing; b) Linear oscillator, where the cubic term is negligible, and the dynamics are dominated by the linear term; c) Cubic oscillator, where the cubic term is non-negligible.

Constant acceleration.
Here we consider the case where the wavevector k in the solutionû of (4.10) is orthogonal toφ * û 0 , and belongs to the null space ofĒ 0 . If, initially,û, as well as its time derivative, are zero, observe that the forcing in the righthand side of (4.10) is also strictly orthogonal to k. Thus, the solution that emerges under this forcing is also orthogonal to k, which transforms (4.10) into the constant acceleration scenario: However, these solutions are likely difficult to observe, since, first, the initial state ofû and its time derivative must be perfect zeros (which is likely only achievable in laboratory conditions, if at all), second,Ē 0 must have an actual null space, and not just a small eigenvalue (which means that the large scale flow must be a perfect one-dimensional jet), and, third, the direction of measurement sampling must lie strictly in the null space ofĒ 0 . As we will see below, even if k only slightly non-orthogonal to the direction of the strong large scale flow, solutions of different kind emerge.

Linear behavior.
If (4.10) can at all behave in the regime of a linear oscillator (that is, the cubic term can be ignored), it is surely achievable when the term k 2 :Ē 0 is maximized, for a given magnitude k . Now, if u 0 is a one-dimensional jet, thenĒ 0 has a large null space, which is comprised of the plane orthogonal to the direction of the jet. Thus, if we intend to maximize k 2 :Ē 0 for a fixed k , the direction of k should be orthogonal to the null space ofĒ 0 , and thus collinear to the direction of the jet. As a result, the two components of the right-hand side forcing in (4.10) in such a regime are collinear. Now, ifû and its time derivative are initially zero, it means that the whole solution of (4.10) is one-dimensional and collinear to k. Thus, to fully describe such a solution, it suffices to examine its projection onto k.
For convenience, let n denote the unit vector in the direction of k, and letû denote the projection ofû onto n: (5.2) n = k/ k ,û = n ·û,û =ûn.
Next, scalar-multiplying (4.10) by n/ k 2 on both sides, we obtain a closed equation for u: where we note that −û 2û * = −|û| 2û . Now, recalling thatφ is real, it is easy to see that, for a purely real, or a purely imaginary strong large scale flowû 0 , the solutionû is also purely real or purely imaginary, and the equation in (5.3) in either of these two cases is a fully real equation for a nonlinear forced oscillator. While it is possible to consider the solutions of (5.3) for general complex large scale flowû 0 , it leads to unnecessary technical complications. Thus, as a proof-of-concept, below we will treat the solutions of (5.3) as purely real. So far, we have chosen the direction of k to maximize the likelihood of a linear behavior, however, it does not automatically mean that (5.3) behaves like a linear oscillatorfor example, even though the linear term in the left-hand side is maximized, the forcing in the right-hand side can still be large enough to cause the cubic term to dominate the linear term. Thus, we need to examine the relationship between the terms by looking at the equilibrium stateû eq of (5.3), by setting the second t-derivative ofû to zero: This is a depressed cubic equation, which has a unique real root via Cardano's formula; namely, denoting (5.5) p = 1 3 n 2 :Ē 0 , q = 2n · (φ * û 0 ), we have, via Cardano's formula, . From Cardano's formula, it is clear that the linear behavior is achieved when p 3 ≫ q 2 ; indeed, in such a case we can expand the formula in Taylor series as which is precisely the equilibrium state for (5.4) without the cubic term. Now observe that, sinceû 0 occupies only a few low-order wavenumbers, for a fixed direction n, we have p as a constant, while q is roughly proportional toφ: 8b) |n · (φ * û 0 )| ∼ |φ| n 2 :Ē 0 , or q 2 ∼ 3pφ 2 .
Now, since n is aligned with the direction of u 0 , we have n 2 :Ē 0 ≈ Ē 0 , and, due to (4.2), the relation above holds automatically. Thus, we conclude that, if the wavevector k above in (4.10) is chosen in the direction of u 0 , then its solution is approximately that of a linear oscillator with the equilibrium state given via In particular, observe that the stronger the large scale flow, the weaker the equilibrium stateû eq . If the initial conditionû(0) and its time derivative are zero, it is easy to see that the solution of (4.10) in the linear regime is given via The period of oscillations is given via Observe that the period decreases (and, therefore, the oscillations become faster) both with increasing k , and with increasing u 0 .

Cubic regime.
Above, we found that choosing k in the direction of the strong large scale flow u 0 resulted in the linear regime. Here, let us choose k which is transversal to u 0 . In this case, the forcing in the right-hand side of (4.10) separates into two components -one collinear to u 0 , and another one collinear to k. In this situation, it is natural to also decompose the solutionû into two (nonlinearly coupled) components -one,û c , collinear to u 0 , and another,û t , collinear to k (and transversal to u 0 ): Somewhat fortunately, the equation forû t alone is fully closed; indeed, observe that the equation in (5.3) forû works irrespectively of the direction of k, and thus we can also use (5.4) to examine its equilibrium state. Now, since k is transversal to u 0 , let us make it almost orthogonal to it; particularly, let us choose k so that (5.14) This leads, according to (5.5) and (5.8), to p 3 ≪ q 2 , which signifies that the cubic term in the oscillator dominates the linear term; subsequently, Cardano's formula in (5.6) becomes (5.15)û t eq = −q 1/3 1 + 1 + and thus (5.16)û t eq =û t eq n ≈ − 4n · (φ * û 0 ) 1/3 n, û t eq ∼ |4φ 0 | 1/3 (n 2 :Ē 0 ) 1/6 k 2/3 .
Here, observe that, for a fixed direction n, the magnitude of the strong large scale flow u 0 amplifies the magnitude of the equilibrium state. Let us also compute an estimate for the period of the oscillations ofû t (t) in the regime of the cubic oscillator. The direct integration of such a nonlinear oscillator, of course, leads to elliptic integrals, and thus we will resort to a crude approximation of the period. First, we center the cubic oscillator at zero, (5.17) d 2ût dt 2 + k 2 (û t ) 3 = 0, and instead assume that the initial condition is offset from zero byû t eq , while its starting time derivative remains zero. This, of course, is not an invariant transformation (unlike the one for the linear oscillator), because, even though the oscillations around zero are symmetric, the oscillations aroundû t eq , which occur between zero and some maximal state, are asymmetric due to nonlinearity. Nonetheless, we disregard this observation to obtain a crude estimate of the period.
The equation above can be integrated at least once, and the period of oscillations can be crudely estimated from there. Multiplying both sides by the time derivative ofû t , we have Integrating, we obtain where the constant of integration is chosen so as to ensure that the time derivative ofû t at the initial time is zero. Now let us make the substitution (5.20a)û t =û t eq cos w, dû t dt = −û t eq sin w ·ẇ, Upon substitution, the sines cancel out on both sides, and we obtain (5.21) (ẇ) 2 = ( k û t eq ) 2 2 (1 + cos 2 w), orẇ = k |û t eq | Here, the variables can be separated, however, that would lead to an elliptic integral. Instead, we obtain the crude approximation to w by observing that Assuming that the cosine part approximately cancels out over time, we approximate The coefficient in front of t is the frequency, thus, the estimate for the period is given via Like in the linear regime, for a fixed n, here the period of oscillations decreases with increasing k , and increasing magnitude of u 0 . However, the rate of increase is slowed down by the relatively low powers of k and u 0 . Thus, while the oscillations keep becoming faster, they do so at a much decreased rate. Now, let us examine the behavior ofû c , which is collinear toφ * û 0 . For that, let us write (4.10) in the form (5.25) d 2ûc dt 2 + d 2ût dt 2 n + k 2 :Ē 0 + (k · (û c +û t n)) 2 (û c +û t n) = −2( k 2 I + k 2 ) φ * û 0 , where we, as before, restrict ourselves to either purely real, or purely imaginary solutions. Subtracting the closed dynamics forû t n (given via (5.3)), we arrive at or, after collecting the terms with same powers ofû c , Here, observe that the component of forcing transversal to the direction ofφ * û 0 is negligible in magnitude (in comparison with the forcing component collinear to the same direction). Thus, we can expectû c to be primarily aligned with the direction of φ * û 0 . If so, let us denote (5.28) n 0 =φ * û 0 φ * û 0 û c =û c · n 0 , and multiply the equation forû c by n 0 / k 2 : (5.29) 1 k 2 d 2ûc dt 2 + n 2 :Ē 0 + 1 + (n · n 0 ) 2 (û t ) 2 û c + 2(n · n 0 ) + (n · n 0 ) 3 û t (û c ) 2 + Noting that, since n and n 0 are almost orthogonal, we have |n · n 0 | ≪ 1, and thus the above equation can be simplified to The equation above is not autonomous with respect to t, sinceû t is, essentially, an external time-dependent periodic forcing, given by the solution of (5.3). Because of this, looking for the explicit solutions of (5.30) is a complicated problem. At the same time, observe that the solutions of (5.30) are "guarded" from above by the cubic oscillatorindeed, since (n · n 0 ) 2 > 0, wheneverû c becomes large enough, the negative acceleration from the cubic term dominates the rest of the terms. Hence, here we will merely estimate the equilibrium stateû c eq of (5.30) under the assumption thatû t =û t eq . Since the latter is given via (5.16), we can estimate (5.31) n 2 :Ē 0 + (û t eq ) 2 ∼ (n 2 :Ē 0 ) 1/3 (n 2 :Ē 0 ) 2/3 + (φ) 2/3 ∼ (û t eq ) 2 , due to (5.14). Thus, for the equilibrium stateû c eq of (5.30), we can estimate, without computing the cubic root exactly, (5.32) (û t eq ) 2ûc eq + 2(n · n 0 )û t eq (û c eq ) 2 + (n · n 0 ) 2 (û c eq ) 3 ∼ 2 |φ 0 | Ē 0 1/2 k 2 .
Observe that the magnitude ofû c eq can be extremely large -for a given potential φ, the right-hand side is multiplied by Ē 0 1/2 (which is large), and divided by (n 2 :Ē 0 ) 1/3 (which is small). In particular, due to (5.14), we can estimate the magnitude ofû c eq below as (5.36) û c eq ≫ Ē 0 1/2 k 2/3 , that is, we can expect the Fourier coefficientsû of (4.10), whose wavevector k is transversal to the large scale flow u 0 , butû itself is aligned with u 0 , to scale, at the very least, as the magnitude of the large scale flow u 0 divided by k 2/3 . Since all linear and nonlinear terms in (5.30) are comparable, computing the period of oscillation is difficult; here we can crudely estimate its scaling by looking at the square root of the equilibrium state of the coefficient at the linear term, |φ 0 | 1/3 k 1/3 (n 2 :Ē 0 ) 1/6 , and observe that the scaling of this crude estimate actually coincides with that of (5.24). a) The scaling estimate in (5.10) corresponds to both k andû k aligned with u 0 ; b) The scaling estimate in (5.16) corresponds to both k andû k nearly orthogonal to u 0 ; c) The scaling estimate in (5.35) corresponds to k nearly orthogonal to u 0 , andû k aligned with u 0 . This combination has the strongest magnitude ofû k , estimated below via (5.36). However, recall that, in the experiments, the measured quantity is the windowed time average of the energy instead. Here, we are going to relate the latter to the former.
To estimate the scaling of the windowed time average of the energy, we refer to the momentum equation in (2.10a), where we additionally use the hydrostatic balance approximation in (3.5). The momentum equation thus becomes Let the time averaging window be denoted via τ. Integrating the above equation between t and t + τ, and dividing by τ, we have In the right-hand side, we use the mean value theorem on ρ: In terms of the Fourier coefficients, we have To estimate the effect of the density above, let us assume that the dynamics occur near the thermodynamic equilibrium with constant θ, so that ρ can be estimated by integrating (3.5): The rate of decay ofρ k can now be estimated from the regularity of the above expression: observe that ρ, as well as all its derivatives, are continuous everywhere, even at zero. Indeed, recalling that φ(x) ∼ 1/ x , the n-th derivative of ρ is given via (6.6) D n ρ(x) = ρ 0 e −2φ/θ C 1 x n+1 + . . . + C n x 2n , for some constants C 1 , . . . , C n .
Clearly, as x → 0, the exponential function decays to zero faster than any negative power of x grows to infinity; thus, all derivatives of ρ(x) above are continuous, and ρ k decays faster than any given negative power of k as k → ∞. Observing thatû k decays as a relatively low negative power of k , we conclude that, near the thermodynamic equilibrium,ρ k decays significantly faster thanû k , and thus can be "factored out" from the convolutions. Further assuming that ρ does not change much over time, that is,ρ k (t) ∼ρ k (t 1 ) ∼ρ k (t + τ), we can cancel out the density and obtain the following estimate: If τ > T (where T is the period of the oscillations for the given scale), then we can expect that, on average, the registered magnitude of the right-hand side will be (6.8) u 2 k ∼ τ û eq k .

Summary
In the current work, we study the dynamics of a pair of particles which interact via a potential. Upon formulating the equations of motion for this pair of particles, we decouple the dynamics by considering separately the motion of the center of mass of the system, and the difference of the coordinates of the particles. We formulate the Liouville equation for the coordinate and velocity difference variables, and then derive the velocity moment transport equations in the same manner as in conventional fluid mechanics. We close these moment equations using the hydrostatic and thermostatic approximations, and arrive at the closed equation for the velocity variable. We show that the velocity equation in the Fourier space possesses a mechanism of formation of time-periodic patterns under the potential forcing, coupled to the strong large scale flow. We find that the scaling of the equilibrium states of these patterns is consistent with the Kolmogorov energy spectrum.
In the course of the derivation, we use the following assumptions and approximations: i) The shear stress and skewness deviator in the transport of the momentum and energy are presumed to be small, in comparison to the corresponding power of the velocity; ii) We use the hydrostatic balance approximation for the transport of the momentum; iii) We use the thermostatic balance approximation (the balance of advective and conductive heat fluxes) for the transport of the energy; iv) We assume that the kinetic temperature of the system is much greater than the energy of the interaction potential; v) For a given Fourier coefficient of the turbulent velocity, the coupling to other Fourier coefficients is discarded, with the exception of the self-coupling, the coupling to the complex-conjugate Fourier coefficient, and the coupling to the strong large scale flow. The strong large scale flow is treated as constant in time. What is left after the simplifications above are applied, reveals what we believe is the fundamental mechanism of the turbulence formation. This mechanism is a cubic oscillator with a linear term, forced by the strong large scale flow. If the strong large scale flow is a jet, we identify several regimes of behavior of this oscillator, depending on the orientation of the wavevector of the given Fourier coefficient relative to the direction of the jet. We find that, in the regime with the largest amplitude of oscillations, the scaling of the time-windowed averages of the kinetic energy corresponds to Kolmogorov's five-thirds law.