Homogeneous shear turbulence – bypass concept via interplay of linear transient growth and nonlinear transverse cascade

We performed direct numerical simulations of homogeneous shear turbulence to study the mechanism of the self-sustenance of subcritical turbulence in spectrally stable (constant) shear flows. For this purpose, we analyzed the turbulence dynamics in Fourier/wavenumber/spectral space based on the simulation data for the domain aspect ratio 1 : 1 : 1. Specifically, we examined the interplay of linear transient growth of Fourier harmonics and nonlinear processes. The transient growth of harmonics is strongly anisotropic in spectral space. This, in turn, leads to anisotropy of nonlinear processes in spectral space and, as a result, the main nonlinear process appears to be not a direct/inverse, but rather a transverse/angular redistribution of harmonics in Fourier space referred to as the nonlinear transverse cascade. It is demonstrated that the turbulence is sustained by the interplay of the linear transient, or nonmodal growth and the transverse cascade. This course of events reliably exemplifies the wellknown bypass scenario of subcritical turbulence in spectrally stable shear flows. These processes mainly operate at large length scales, comparable to the box size. Consequently, the central, small wavenumber area of Fourier space (the size of which is determined below) is crucial in the self-sustenance and is labeled the vital area. Outside the vital area, the transient growth and the transverse cascade are of secondary importance - Fourier harmonics are transferred to dissipative scales by the nonlinear direct cascade. The number of harmonics actively participating in the self-sustaining process (i.e., the harmonics whose energies grow more than 10% of the maximum spectral energy at least once during evolution) is quite large - it is equal to 36 for the considered box aspect ratio - and obviously cannot be described by low-order models.


Introduction
In the 1990s, the nonnormal nature of shear flows and its consequences became well understood by the hydrodynamic community (see e.g., Refs. [1,2,3,4,5,6,7]). As a result, a new, socalled bypass concept emerged for explaining the onset and sustenance of turbulence in spectrally stable shear flows [8,9,10,11,12,13,14,15,16]. This concept is based on the linear transient, or nonmodal growth of vortical perturbations due to the nonnormality/shear, which is the

Physical model and equations
The motion of an incompressible fluid with constant viscosity, ν, is governed by the Navier-Stokes equations ∂U ∂t where ρ is the density, U is the velocity and P is the pressure. Equations (1)-(2) have a stationary equilibrium solution -the flow in the x-direction with the linear shear of mean velocity in the vertical y-direction, U 0 = (Sy, 0, 0). Without loss of generality, the constant mean shear rate S is chosen to be positive. Equilibrium density ρ 0 and pressure P 0 are spatially constant. Such an idealized configuration of a flow with a linear shear, despite its simplicity, allows us to grasp key effects of shear on the perturbation dynamics and ultimately on a resulting turbulent state in shear flows.
Consider 3D perturbations of the velocity and pressure, u and p, about the equilibrium. Representing the total fields as the sum of the equilibrium and perturbed values, U = U 0 + u, P = P 0 + p, substituting these into Eqs. (1)- (2), we obtain the following system of equations governing the dynamics of perturbations with arbitrary amplitude From Eqs. (3)-(6) we derive a dynamical equation for the perturbation kinetic energy in order to investigate the energy balance in the self-sustained turbulent state. In the presence of the imposed shear-periodic boundary conditions (see Sec. III below), the evolution of the volume-averaged kinetic energy density, E = ρ 0 u 2 /2, is governed by the equation: where the angle brackets, ... denote a volume average over an entire domain at a given time.
The first term on the right hand side of Eq. (7) is the Reynolds stress and describes the exchange of energy between perturbations and the background flow. Note that this term originates from the linear term proportional to shear (−Su y ) on the right hand side of Eq. (3). The second term is negative by definition and describes energy dissipation due to viscosity. It is evident that there is no contribution from the nonlinear terms in the total energy evolution Eq. (7). Thus, only Reynolds stress, if negative, can supply perturbations with energy, extracting it from the mean flow due to shear. In turbulence studied below, the Reynolds stress ensures energy injection into turbulent fluctuations, balancing viscous dissipation. The nonlinear terms, not directly tapping into the shear flow energy and therefore not changing the total perturbation energy, only redistribute energy gained by means of the stress among Fourier harmonics of perturbations with different wavenumbers (see below). Obviously, in the absence of shear (S = 0), the contribution from the Reynolds stresses disappears in Eq. (7) and hence the total perturbation energy cannot grow, gradually decaying due to viscosity.

Spectral representation of the dynamical equations
We derive dynamical equations for the quadratic forms of velocity perturbations, (u 2 x , u 2 y , u 2 z ), and kinetic energy density, E, in Fourier space to gain a deeper insight into the interplay of the flow shear and nonlinearity in the self-sustenance of perturbations (turbulence). For a spectral representation of the main equations, we decompose the perturbations into spatial Fourier harmonics where f ≡ (u, p) denotes the perturbations andf ≡ (ū,p) is their corresponding Fourier transforms (d 3 k ≡ dk x dk y dk z ). Substituting decomposition (8) into Eqs. (3)-(6), we obtain the following equations governing the dynamics of perturbation harmonics in spectral space (the constant density is set to unity, ρ 0 = 1) where k 2 = k 2 x + k 2 y + k 2 z . These spectral equations contain the linear as well as the nonlinear, Q(k, t) = [Q x (k, t), Q y (k, t), Q z (k, t)], terms that are the Fourier transforms of corresponding linear and nonlinear ones in the original Eqs. (3)- (6). The latter are given by and describe nonlinear triad interactions among harmonics with different wavenumbers in Fourier k-space. From Eqs. (9)-(12) one can eliminate the pressurē and rewrite them in the following form: where H k = −(ū xū * y +ū * xūy )/2 is the spectrum of the Reynolds stress. Using Eqs. (16)-(18) together with the incompressibility condition (12), we obtain the evolution equation for the spectral energy density where D k = 2νk 2 E k describes energy dissipation and N x , N y , N z are the modified nonlinear transfer functions, Equations (16)- (19) fully determine the nonlinear dynamics of the considered system in Fourier space and are the basis for subsequent analysis. According to them, four basic processes underly the dynamics of the perturbations and spectral energy: (i) The second terms with Sk x ∂/∂k y in the brackets on the lhs of Eqs. (16)- (19), are the fluxes of the corresponding quantities parallel to the k y -axis. These terms are of linear origin, coming from the convective derivative on the lhs of the main Eqs. (3)-(5) and therefore correspond to the advection by the mean flow. In other words, due to the background shear flow, the spectral quantities (Fourier transforms) "drift" in k-space: harmonics with k x > 0 and k x < 0 drift in opposite directions along the k y -axis at a speed S|k x |, whereas the ones with k x = 0 are not advected by the flow. Since d 3 k∂(k x E k )/∂k y = 0, the drift only transports harmonics parallel to the k y -axis without changing the total kinetic energy. (ii) The first terms on the rhs of these equations are associated with shear, i.e., they originate from the linear term proportional to the shear rate on the rhs of Eq. (3) and describe influence of the mean flow on the dynamics of harmonics. Therefore, the first term on the rhs of spectral energy Eq. (19) describes exchange of energy between the mean flow and individual harmonics. This term in the spectral energy equation is related to the volumeaveraged Reynolds stress entering Eq. (7) through and hence serves as a main source of energy for harmonics at the expense of which they can undergo amplification. This shear, or nonnormality induced process is in fact linear by nature and has a transient character due to the drift in k-space [28,29,30,18]. In the case of turbulence studied in this paper, H k describes injection of kinetic energy into turbulent fluctuations as a function of wavenumbers. The first terms on the rhs of Eqs. (16)- (18) describe the same process for the quadratic forms of perturbation velocity components. The harmonics, drifting parallel to the k y -axis, go through dynamically important regions in spectral space at small and intermediate wavenumbers, where energy-supplying linear terms and redistributing nonlinear terms cooperate to sustain turbulence (see e.g., Fig. 8). (iii) The second terms on the rhs of these equations, N x , N y and N z , as mentioned above, describe, respectively, nonlinear transfers, or redistributions of the spectral amplitudes of the velocity components |ū x |, |ū y | and |ū z | in k-space. Their net effect in the spectral energy budget over all wavenumbers is zero at all times: which is, in fact, a consequence of vanishing of the nonlinear terms in total energy equation in real space (Eq. 7). Although nonlinearity does not directly change the total (i.e., summed over all wavenumbers) spectral energy, it plays a central role in shear flow turbulence. N k determines transfer, or cascade of the spectral energy in k-space. This energy is drawn from the mean flow by the linear stress term, leading to the development of its specific spectra. These transfer functions are one of the main focus of the present study. We aim to explore how they operate in the presence of shear by adopting the approach developed and used in [17,19,20], the main idea of which is performing a full 3D Fourier analysis of individual terms in the dynamical Eqs. (16)- (18), thus allowing for anisotropy of spectra and cascades. In particular, we will show that, like that in 2D shear flows considered in those studies, nonlinear transfers in the 3D case result in the redistribution of spectral energy among wavevector angles in k-space, which we refer to as NTC. We demonstrate that the interplay of NTC with linear phenomena (i.e., with the first rhs terms) becomes intricate: it realizes positive feedback, that is, repopulates those harmonics that have the potential of transient growth, thus contributing to the turbulence's self-sustenance. (iv) The third terms on the rhs describe viscous dissipation of velocity. Comparing them with the energy-supplying linear terms, we can see that viscous dissipation becomes important at large wavenumbers k ≥ k D = S/ν, where k D denotes the effective wavenumber for which dissipation effects start to play a role.

Results
Now we turn to an analysis of nonlinear evolution of perturbations employing numerical methods. Main emphasis will be on the spectral aspect of the dynamics using the mathematical formalism outlined in the previous section. We solve full Navier-Stokes equations (3)-(6) in a rectangular box with size 0 ≤ x ≤ L x , −L y /2 ≤ y ≤ L y /2, 0 ≤ z ≤ L z , in the streamwise (x), vertical/shearwise (y) and spanwise (z) directions, respectively, using N x × N y × N z number of grid points. Numerical simulations are performed with the code developed by Sekimoto et al. [27], which solves the Navier-Stokes equations in the velocity-vorticity representation [31]. The streamwise and spanwise directions are periodic so that the spectral method is applied. The vertical direction is shear-periodic, realized by adding "shift" factors in the compact-finite-difference matrices which guarantee spectral like resolution [32]. The method does not require the remeshing procedure leading to constant loss of enstrophy [33], which has to be avoided when studying vortical structures. To assess how well our simulations are resolved, we compare the Kolmogorov length scale η = (ν 3 / ε ) 1/4 , where ε = ν (∇ × u) 2 is the volume-averaged dissipation rate of the perturbation energy, to the grid cell sizes Since in the y-direction the boundary conditions are shear-periodic, standard FFT technique cannot be applied for calculating Fourier transforms along this direction during post-processing. Due to shear flow, harmonics drift along the k y -axis and, consequently, k y of an individual harmonic changes with time linearly, k y (t) = k y (0) − Sk x t (see e.g., Refs. [18,28]) and remains consistent with the shear-periodic boundary conditions. Following [34,35], we modified the standard FFT technique to take this effect into account as follows. As k y drifts, at each time t in the corresponding computational domain in (k x , k y )-plane (at fixed k z ) there are only Fourier harmonics with wavenumbers So, in this case Fourier transformation reduces to the decomposition into these wavenumbers. We calculate it by first transforming a physical quantity along k x and k z using the standard FFT, next by multiplying the result by exp(iSk x t ), that accounts for the timedependence of k y , and applying again FFT method along k y .
The time is normalized by the dynamical, or shear time S −1 , while an arbitrary length scale can be used for normalizing distances and other physical variables, since the problem is in fact scale-free. The Reynolds number is defined in terms of L z as Re z = SL 2 z /ν. Together with the Reynolds number, the problem has two dimensionless free parameters -aspect ratios of the domain, In this paper, we consider the model with the cubic box A xz = A yz = 1, resolution N x ×N y ×N z = 128×128×128 and the Reynolds number Re z = 4930. In this case, we found that Δx ≈ 1.66η, which implies that larger and intermediate scales, where, as we will show, the self-sustaining dynamics of turbulence is concentrated, are well resolved in this box. We normalize the wavenumbers k x , k y , k z , respectively, by Δk Since Δk x , Δk y and Δk z are the grid cell sizes in Fourier space, the normalized streamwise and spanwise wavenumbers become integers k x , k z = 0, ±1, ±2, ..., while k y , although changes with time due to drift, is integer at discrete moments t m = mA xz /(|k x |A yz ), where m is also a positive integer.
Starting from initial random perturbations, the flow settles into a self-sustained turbulent state. Our main goal is to investigate underlying self-sustaining process of the turbulence. We show that the turbulence is self-sustained by the interplay of linear and nonlinear processes and is characterized by quasi-periodic bursting events, which are related to the dominant mode. This scenario of a self-sustained state, based on a subtle cooperation between linear and nonlinear processes, is a keystone of the bypass concept of subcritical turbulence in spectrally stable shear flows [9,11,12,13,36,14]. Finally, it should be stressed that we base our study on DNS results of turbulence and in this respect is general, not relying on simplifying assumptions usually made in low-order models of self-sustaining processes (e.g., [24,25]).
It is well-known from the linear nonmodal theory of shear flows that optimal perturbations which are characterized by maximum transient growth during the dynamical -eddy turnover time, which can be assumed of the order of shear time, t d S −1 (or, t d 1 in normalized units) are responsible for most of the energy extraction from the mean flow [5,37]. In Fig.  1 we present the transient growth factor of energy for optimal harmonics during t d , which is calculated from the linearized inviscid equations, i.e., using only the first linear terms on the rhs of Eqs. (13)- (15). These terms, and hence the dynamics of harmonics in the inviscid case, depends only on the ratio of wavenumbers, k y /k x and k z /k x . In this figure, we show the growth factor as a function of these ratios. As one can see, the optimal growth factor is largest at 1 ≤ k z /k x < 4 and 0 < k y /k x < 1. So, one can cover as many of these optimal perturbations as possible and thus better account for their major role in the energy exchange processes between turbulence and mean flow if the computational box aspect ratios satisfy the conditions A xz < A yz and A xz ≥ 1. In agreement with the previous study of homogeneous shear turbulence by Pumir [26], in our simulations we also find that in addition to the nonsymmetric (k x = 0) optimal harmonics, a harmonic with a large scale spanwise variation, which   is uniform in the streamwise and shearwise directions, k x = 0, k y = 0, k z = ±1, plays a central in the dynamical processes. Due to its importance in turbulence dynamics, we call this harmonic as the basic mode. Other harmonics, whose energy grows more than 10% of the maximum value of the spectral energy at least once during the evolution, are considered to be actively participating in and forming the self-sustaining dynamics of turbulence. We show that quasiperiodic bursts of energy and Reynolds stress characteristic of homogeneous shear turbulence [26] are determined by a competition between the basic mode and other dynamically "active" harmonics. The number of these harmonics generally depends on the box aspect ratio and is quite large (in the considered here cubic box, it is 36, see Fig. 5, but can reach more than 200 for steamwise elongated boxes), indicating fairly complex dynamics, not amenable to reduced low-order description.  Below we analyze in detail the specific energy spectra and self-sustaining dynamics of turbulence in wavenumber/spectral space by calculating the Fourier transforms of linear and nonlinear terms in governing equations using the simulation data.

Nonlinear evolution
At the early stage of evolution, each harmonic contained in the initial perturbations grows transiently due to the flow nonnormality, leading to the growth of the total energy and Reynolds stress. Then after about t s = 30S −1 , the amplitudes of the perturbations become sufficiently large in the nonlinear regime and the flow settles into a self-sustained turbulence. As was found in the paper by Pumir [26], in the case of cubic box, the basic harmonic/mode with wavenumber k b = (0, 0, ±1), dominates over other ones, i.e., its spectral energy, E k (k b ), and stress, H k (k b ), are, respectively, maxima of E k and H k during an entire run. Figure 2 shows the time-evolution of the volume-averaged total energy and stress as well as the energy and stress of the basic mode in the saturated state, which appear to display similar patterns of quasi-periodic bursts in time. The volume-averaged Reynolds stress, − u x u y , and the stress of the basic mode, H k (k b ) are positive at all times and non-decaying as required by the self-sustaining process according to Eq. (7). The bursts in the total energy and stress closely follow amplifications (peaks) of E k (k b ) and stress H k (k b ) in the basic mode. This indicates that these bursts are related to the dynamics of the basic mode, which thus appears to be mainly responsible for the energy pumping from the mean flow into turbulence. During the bursts, the contribution of H k (k b ) to the total Reynolds stress is as much as 60% -70%, while other modes become important (although still remain smaller than the basic mode) in quiescent intervals between the bursts, when the total energy and stresses as well as the energy and stress of the basic mode are low. One can conclude, that most of the turbulent energy is primarily produced through the interaction of the basic mode with the mean shear flow.
The structure of the velocity components in the fully developed turbulence in physical space is presented in Fig. 3. The evolution of the volume-averaged perturbation kinetic energy and quadratic forms of streamwise, u 2 x /2, shearwise/vertical, u 2 y /2, and spanwise, u 2 z /2, velocities as well as the Reynolds stress are shown in Fig. 4. These plots show the dominance of the streamwise velocity over the shearwise and spanwise ones. In the (x, z) slices of u y and especially of the u x velocity, we can clearly see the signatures of the dominant basic harmonic with small scale fluctuations superimposed on it due to other harmonics.
To understand the role of other modes in the turbulence dynamics and energy exchange process with the mean flow, it is important to define the main, energy-containing area in kspace, which we also refer to as the vital area of turbulence. For this purpose we identify two sets of harmonics whose energies grow from 5% to 10% and larger than 10% of the maximum Black dots represent the harmonics whose energies grow more than 10% of the maximum spectral energy at least once during evolution. Bigger black dot corresponds to the basic mode having the maximum spectral energy at all times. Circles are the modes whose energies grow from 5% to 10% of the maximum spectral energy at least once during evolution. spectral energy at each time at least once during the evolution. These dominant energy-carrying, or dynamically important harmonics are presented in Fig. 5. They have small wavenumbers |k x | ≤ 2, |k y | ≤ 3, |k z | ≤ 3 (at |k z | = 3, not shown in this figure, mode energies always remain less than 10% of the maximum energy) which basically form the vital area of turbulence. Harmonics with larger wavenumbers |k x | > 2, |k y | > 3, |k z | > 3 lie outside the vital area and always have energy and stress less than 5% of the maximum value, therefore, not playing a major role in the energy exchange process between the mean flow and turbulence. Essentially, only these large scale modes in the vital area take part in the self-sustaining dynamics of turbulence, but the basic mode/harmonic still remains dominant -the basic harmonic is most active among the active harmonics. We emphasize that the total number of large scale "active modes", which we define as the harmonics whose energies grow more than 10% of the maximum value at least once during the evolution (black dots in Fig. 5) is equal to 36, implying that the dynamics of homogeneous shear turbulence cannot by any means be reduced to low-order models of selfsustaining processes. Figure 6 shows the time-averaged spectra of the kinetic energy in (k x , k y )-plane at different k z . These time averages are done over an entire time of the saturated turbulence in the simulation. The spectra are anisotropic due to shear -the isolines are of elliptical shape inclined opposite to the k y -axis. This fact indicates that on average modes with k y /k x < 0 have more energy than those with k y /k x > 0 at fixed k x . The plots show that the time-averaged spectral kinetic energy at fixed k z is larger for small k x and k y and decays at least by two order of magnitude at k 2 x + k 2 y ≈ 2. So, the active modes are situated just within this range of small wavenumbers, as also seen in Fig. 5. As mentioned above, the maximum of the spectral energy is reached at k x = k y = 0, k z = ±1 corresponding to the basic mode. The spectral energy is much reduced already for the next spanwise modes, k z = 2, and further decreases with increasing k z .
The integrated in (k x , k y )-plane the time-averaged spectral energy,Ê(k z ) = E k dk x dk y , and stress,Ĥ(k z ) = H k dk x dk y , as a function of k z are shown in Fig. 7. They both reach a maximum at k z = ±1 and rapidly decrease with k z .
The anisotropic nature of the kinetic energy spectrum, which clearly distinguishes it from typical isotropic turbulent spectra in the classical shearless case of Kolmogorov phenomenology, arises as a result of the specific interplay of linear and nonlinear processes in k-space. We show below that these processes are anisotropic over wavenumbers due to shear, resulting in a new phenomenon -the NTC of power in spectral space.

Energy injection, nonlinear transfers and their interplay in k-space
For the better understanding of the character of the anisotropic kinetic energy spectrum and nonlinear transfers, in Fig. 8 we present the distribution of the time-averaged energy injection, i.e., the Reynolds stress spectrum, H k , and the nonlinear transfer functions, N x , N y and N z , in k-space. Since these quantities are symmetric with respect to a change k → −k, without loss of generality, hereafter we concentrate on the upper part (k z > 0) of k-space. From these figures it is clear that, like the spectral energies, all these quantities exhibit anisotropy over wavenumbers, that is, depend on wavevector angle, or orientation. H k is primarily concentrated in the central part of spectral space, |k x | ≤ 2, |k y | ≤ 3, |k z | ≤ 2 and k y /k x > 0 -in the vital area, where it supplies the turbulence with energy via nonnormality-induced transient amplification of the dynamically active harmonics (see below). At k z = 0, H k is larger and mainly operates around |k y |, |k x | ∼ 1 (Fig. 8a), it is positive at k x /k y > 0, supplying harmonics with energy there and negative at k x /k y < 0, returning energy to the flow. On the other hand, at k z ≥ 1, the distribution of H k is changed -it is positive and appreciable around k x = k y = 0 (Figs. 8b and  8c). The Reynolds stress spectrum is highest, and therefore the energy extraction from the mean flow into turbulence occurs most intensively, at k z = 1 with the maximum at k b corresponding to the basic mode (Figs. 7b and 8b). Outside the vital area, at |k x | > 2, |k y | > 3, |k z | > 2, the stress and hence the energy input rapidly decrease and vanish at large wavenumbers (Fig. 7b). The total energy input is negligible at k z = 0,Ĥ k ≈ 0, as seen from Fig. 7(b), since H k changes sign at this k z in (k x , k y )-plane (Fig. 8a). Thus, the main energy supply for turbulence comes from large-scale modes, with wavelengths comparable to the box size. These are the modes H k is appreciable, and hence energy injection into turbulence mainly occurs, in the vital area at small wavenumbers, |k x | ≤ 2, |k y | ≤ 3, |k z | ≤ 2, on the k y /k x > 0 side. At k z = 1, H k reaches largest values and hence energy injection is most intensive at this spanwise wavenumber. N x , N y and N z terms transfer, respectively, the streamwise, shearwise and spanwise components of the spectral kinetic energy anisotropically, or transversely over wavevector angles in Fourier space, away from regions where they are negative N x < 0, N y < 0, N z < 0 (blue) to regions where they are positive N x > 0, N y > 0, N z > 0 (yellow and red). which most effectively interact with the mean flow and gain energy, whereas smaller scale modes only receive energy from the former due to nonlinear transfers. Below we describe the character of nonlinear transfers in the presence of shear and their interplay with the energy-injecting stress term, leading to a self-sustained turbulent state. Note also that the energy injection into turbulence mediated by the stress spans a range wavenumbers and in this respect differs from classical cases of forced turbulence in flows without mean shear, where forcing is applied at a specific, or a narrow band of wavenumbers (see e.g., [38]).
The nonlinear terms N x , N y and N z describe, respectively, transfers of the energies in the streamwise, shearwise/vertical and spanwise velocity components. As noted above, they do not represent a source of the total energy for the turbulence, but only redistribute energy of harmonics, which is extracted from the mean flow, over wavenumbers and together with the energy-injecting stress term, determine the spectra. It is seen from Figs. 8(d)-(l) that different components of the nonlinear term operate differently in k-space and hence play different roles in the self-sustaining process. For this reason, we plot them separately instead of their sum, which describes nonlinear transfer of energy in Eq. (19). Since we are interested in the interplay of linear energy-injecting stress terms and nonlinear cascades, in this figure we show the nonlinear transfer terms also at k z = 0, 1, 2 for which the stress is significant. First of all, these plots clearly demonstrate the anisotropic nature of the nonlinear processes in k-space. The transfer functions depend on wavevector angle and, consequently, redistribute the energies generally not only along wavevector, corresponding to direct/inverse cascades that are well-known in the classical turbulence theory, but also transverse to wavevector. This angular redistribution is the essence of NTC and occurs from the regions in k-space where the respective transfer functions are negative N x , N y , N z < 0 (blue areas in Figs. 8d-l) to regions where they are positive N x , N y , N z > 0 (yellow and red areas); at all other wavenumbers (light green areas) these functions are small. It is seen from Fig. 8(d)-(l) that the NTC primarily operates in the central vital area of the wavenumber space, |k x | ≤ 2, |k y | ≤ 3, |k z | ≤ 2, where, as described above, most of the energy supply process by H k is concentrated. Outside the vital area at large wavenumbers, where the nonnormality-induced energy exchange between the mean flow and harmonics vanishes, the transfer functions become nearly independent of the wavevector angle and depend mostly on its magnitude, i.e., the cascade is no longer transverse but direct -the situation is similar to the classical homogeneous isotropic turbulence described by Kolmogorov phenomenology. Below we describe the interplay of the NTC with this energy-injecting linear phenomena that can realize either positive or negative feedback. In the case of positive feedback, the NTC repopulates harmonics having the potential of transient growth.
To characterize transfers also along k z due to the NTC, we integrate N i over (k x , k y )-plane,  , y, z, k), and represent them as functions of k z in Fig. 9. Figures  8(d)-(f) and 9(a) show that N x is negative and large by absolute value in the vital area, where the stress, H k , is positive and appreciable, but it is small and positive at larger wavenumbers outside this area, |k x | > 2, |k y | > 3, |k z | > 3, where the stress is vanishing. In other words, N x provides a sink, or negative feedback for H k by taking out energy in the streamwise velocity component gained from the mean flow and transferring it to large wavenumbers and part of it to k z = 0, as evident from the dependence ofN x on k z . By contrast, N y is positive in those regions of the vital area where H k is also positive (Figs. 8g-8i and 9b). As a result, at these wavenumbers, it provides positive feedback by generating the energy in the shearwise velocity, |ū y | 2 /2, which, in turn, leads to the growth of the energy in the streamwise velocity, |ū x | 2 /2, due to the nonnormality/shearinduced linear growth mechanism (the first term of linear origin on the rhs of Eqs. 13 and 16, while N x < 0 at these wavenumbers and opposes this growth). This streamwise energy, which in fact appears to be largest among the other two components and hence giving the main contribution in the total energy, is then transferred to larger wavenumbers by N x . Note that this positive feedback by N y -crucial to the self-sustenance of turbulence -is realized owing to the NTC, which ensures repopulation with the shearwise energy of those harmonics with the wavevector orientations (k y /k x > 0) having H k > 0 and hence the capability of nonmodal growth. Figures 8(j)-(l) and 9(c) show N z , which is responsible for the transfer of energy of the spanwise velocity. However, since the Reynolds stress, providing energy supply for the turbulence, is produced by the product (correlation) of u x and u y , below we concentrate only on N x and N y functions governing nonlinear redistribution of energies in these velocity components. Figure 9(d) shows the nonlinear transfer function for the energy,N k =N x +N y +N z . It is dominated byN x and hence the nonlinear transfer of the total energy is similar to that of the streamwise energy, which in fact gives the dominant contribution to the former. Thus, the transverse cascade of energy appears to be a generic feature of nonlinear dynamics of perturbations in spectrally stable shear flows. The conventional description of shear flow turbulence solely in terms of direct and inverse cascades, which leaves such nonlinear transverse cascade out of consideration, might be incomplete and misleading. It should be emphasized that revealing of the complete picture of these nonlinear cascade processes has become largely possible due to carrying out the analysis in Fourier space. Because of the shear-induced anisotropy of cascade directions, often used angle-integrated energy spectra and transfer functions alone, clearly, are not fully representative of the actual, more general, over wavevector angles, nonlinear redistribution of the spectral energies in k-space.
We have described above how the interplay of linear and nonlinear process happens in a time average sense in the fully developed saturated turbulence. To gain a better insight in its self-sustenance and appearance of the bursts in the total energy evolution seen in Fig. 2, below we analyze the dynamics of basic harmonic. Namely, we show that the alternating growth and decay intervals in its evolution defines the behavior of the bursts.

The dynamics of the basic mode
We consider the dynamics of the dominant basic mode with k b = (0, 0, ±1), which during the whole evolution corresponds to the maximum energy and stress in spectral space. For this mode, the dynamical Eqs. (16)-(18) with S = 1 are reduced to the following system: ∂ ∂t ∂ ∂t  Figure 10. Evolution of (a) the streamwise, |ū x | 2 /2, and (b) shearwise, |ū y | 2 /2, spectral energies as well as the nonlinear transfer functions (c) N x and (d) N y for the basic mode with k b = (0, 0, ±1). The streamwise energy is much larger than the shearwise one due to nonnormality-induced linear growth. N x is negative during the whole evolution and acts as a sink for the streamwise energy, while N y oscillates and, when positive, causes increase in |ū y |; its time-averaged value over the evolution is positive 4.78 · 10 −4 .  Figure 11. Evolution of (a) the stress production rate, dH k (k b )/dt, and nonlinear transfer term N y (k b ) as well as (b) the spectral Reynolds stress, H k (k b ), and shearwise velocity, |ū y (k b )| of the basic mode. The stress production rate closely follows N y (k b ) (we scaled it by a factor of 2 to make comparison easier) and, consequently, the amplifications (bursts) of both shearwise velocity and stress occur when N y is positive and decrease when N y is negative.
Note that N z (k b , t) = 0 and Eq. (23) can be solved directly to get exponential decay of spanwise velocity of the basic mode,ū z (k b , t) =ū z (k b , 0) exp(−2νt). So, after a sufficiently long period of time the spanwise velocity vanishes and the velocity field of the basic harmonic becomes two-dimensional. Figure 10 shows the time development of the squared velocity amplitudes, |ū x (k b )| 2 /2 and |ū y (k b )| 2 /2 and the nonlinear transfer functions N x (k b ) and N y (k b ) for the basic mode. Because it is of large scale and the Reynolds number is high, the dissipation terms in these equations appear to be much smaller than the transfer functions. The velocities exhibit a sequence of amplifications (bursts) and decays, with the streamwise velocity being always much larger than the shearwise one. N x (k b ) is always negative and therefore acts as a sink in the streamwise velocity evolution Eq. (21), while N y (k b ) oscillates irregularly with alternating sign. When N y (k b ) is positive, it amplifies the shearwise velocity, while when negative causes its decrease, but the time-average of N y (k b ) is positive, though small (see also Fig. 8h), indicating the growth trend of the shearwise velocity on average. The shearwise velocity, in turn, leads to the amplification of the streamwise velocity due to the shear-related linear first term on the rhs of dynamical Eq. (13). This is essentially the nonnormality induced linear nonmodal amplification process, resulting in the growth of streamwise velocity much larger than that of the shearwise one. Note that this is the only cause of amplification for |ū x (k b )|, since the corresponding nonlinear term, N x (k b ), as mentioned above, is always negative and cannot amplify it. Then,ū y (k b ) and the amplifiedū x (k b ) jointly give rise to positive stress, H k (k b ) > 0 (Figs. 2b and 11b). As a result, the stress production rate, dH k (k b )/dt (calculated by combining Eqs. 13 and 14), closely follows, or correlates with N y (k b ) (Fig. 11a, the correlation coefficient is R(dH k /dt, N y ) = 0.93 1 and, consequently, the stress itself turns out to follow |ū y (k b )| (Fig. 11b, R(H k , |ū y |) = 0.97), i.e., the peaks and dips in the time evolution of these quantities coincide. So, it can be said that transfer function N y (k b ) plays a central role in the dynamical behavior of the basic mode by providing a positive feedback to the nonmodal growth. When N y (k b ) > 0, it causes the amplifications ofū y , which via nonmodal growth mechanism leads to amplification ofū x and of the corresponding stress. This bursts alternate with the decaying intervals of these quantities corresponding to negative N y (k b ) < 0. Because for the aspect ratio (1, 1), the basic mode is the dominant one among all the active modes, these bursts also manifest themselves in the total stress evolution, as evident in Fig. 2b. The same cooperative action of linear transient growth and nonlinear cascade by N y underlies the dynamics of other symmetric (k x = 0) and non-symmetric (k x = 0) active modes in the vital area, with wavenumbers in the vicinity of k b (Fig. 5), however, its quantitative character depends on mode wavenumbers. From the time-averaged plots in Fig. 8, we see that on average the Reynolds stress of modes with k y /k x > 0 is larger than that of the modes with k y /k x < 0 and hence the former modes play a main role in supplying turbulence with energy. From these plots we also see that on average N y > 0 at k y /k x > 0 in the vital area, i.e., the NTC ensures continual repopulation of just these quadrants of Fourier space with new harmonics that contribute most to positive stress and energetically maintain the turbulent state.

Summary and discussion
The considered here background flow with constant mean shear is linearly stable (with asymptotically decaying linear perturbations) according to classical/modal stability theory. Therefore, the only energy source for turbulence's self-sustenance in such a flow comes from linear transient amplification of perturbations due to the nonnormality associated with the flow shear. To understand the mechanism of homogeneous shear turbulence's self-sustenance, we performed DNS and then analyzed the turbulence dynamics in Fourier/wavenumber space based on the DNS data. We Fourier transformed Navier-Stokes equations and derived evolution equations for the perturbed streamwise, shearwise and spanwise velocities in wavenumber space. In these spectral equations, using the simulation results, we calculated individual terms, which are classified into two types -terms of linear and nonlinear origin. The term of linear origin, specifically the Reynolds stress, is responsible for energy exchange between the turbulence and the mean flow through transient, or nonmodal amplification of perturbation harmonics due to the flow shear. The nonlinear terms, which do not directly draw the mean flow energy, act to transversely, over wavevector angles, redistribute this energy in Fourier space, continually repopulating those perturbation harmonics which can undergo transient growth. This new type of transversal redistribution of harmonics referred to as the nonlinear transverse cascade appears to be a generic feature of nonlinear dynamics of perturbations in spectrally stable shear flows. It was also demonstrated to take place in 2D HD and MHD shear flows [19,20] and fundamentally differs from classical concepts of energy cascade processes (direct/inverse) in turbulence theory. So, the conventional description of shear flow turbulence solely in terms of direct and inverse cascades, which leaves such nonlinear transverse cascade out of consideration, might be incomplete and misleading. Thus, we have demonstrated that in the considered constant shear flow, turbulent state is sustained by the subtle interplay of linear and nonlinear processes -the first supplies energy for turbulence via shear-induced transient growth process (characterized by the Reynolds stresses) and the second plays an important role of providing a positive feedback that makes this transient growth process persist over long times and be self-sustaining against viscose dissipation.
We found that harmonics that undergo significant transient amplification ("active modes"), and thus are responsible for the majority of energy supply for turbulence, have larger length scales, comparable to the box size. According to the presented analysis, the transverse cascade also operates in this range of length scales. So, this central, small wavenumber area of wavenumber space is the arena of activity of the energy supply linear and redistributing nonlinear processes and their interplay. Consequently, this area is the crucial in the self-sustenance of turbulence and we call it the vital area.
The number of these active harmonics (e.g., the harmonics whose energies grow more than 10% of the maximum spectral energy at least once during evolution) in the vital area generally depends on the box aspect ratio and is quite large (in the considered here box size, it is 36, Fig.  5), indicating fairly complex dynamics, which evidently cannot be described by low-order models. Thus, carrying out a full DNS, we demonstrated the realization of the transverse cascade and its role in the self-sustenance of homogeneous shear turbulence, without truncating the number of modes as usually done in low-order models of self-sustaining processes (e.g., Refs. [24,25]).
We found that for the considered cubic box, the basic mode with the wavenumber k b = (0, 0, ±2π/L z ) plays a special role. It yields the maximum Reynolds stress during the entire evolution and undergoes the largest transient amplification among the nearby active modes. The amplification of the basic mode occurs in bursts, alternating with decaying intervals. The basic mode contributes up to 70% of the total Reynolds stress during the bursts and thus largely determines the behavior of the turbulence's total energy - Fig. 2 shows that the timehistory of the total energy and the stress closely follow those of the basic mode. However, at other aspect ratios (A xz , A yz ) = (3, 2), (1, 2) that we also examined, other nearby active modes grow comparably with the basic mode and hence their contribution to the total energy increases. Nevertheless, the largest bursts in the total energy and stress correspond to the growth (upsurge) of the basic mode for these aspect ratios too.