Evolution of channel flow and Darcy’s law beyond the critical Reynolds number

For incompressible channel flow, there is a critical state, characterized by a critical Reynolds number Rec and a critical wavevector mc along the channel direction, beyond which the channel flow becomes unstable in the linear regime. In this work, we investigate the channel flow beyond the critical state and find the existence of a new fluctuating, quasi-stationary flow that comprises the laminar Poiseuille flow superposed with a counter-flow component, accompanied by vortices and anti-vortices. The net flow rate is reduced by  ~ 15% from the linear, laminar regime. Our study is facilitated by the analytical solution of the linearized, incompressible, three-dimensional (3D) Navier–Stokes (NS) equation in the channel geometry, with the Navier boundary condition, alternatively denoted as the hydrodynamic modes (HMs). By using the HMs as the complete mathematical basis for expanding the velocity in the NS equation, the Rec is evaluated to 5-digit accuracy when compared to the well-known Orszag result, without invoking the standard Orr-Sommerfeld equation. Beyond Rec, the analytical solution is indispensable in offering physical insight to those features of the counter-flow component that differs from any of the pressure-driven channel flows. In particular, the counter flow is found to comprise multiple HMs, some with opposite flow direction, that can lead to a net boundary reaction force along the counter-flow direction. The latter is analyzed to be necessary for satisfying Newton’s law. Experimental verification of the predictions is discussed. Graphical Abstract Supplementary Information The online version contains supplementary material available at 10.1140/epje/s10189-023-00289-4

Consider the incompressible channel flow between two parallel solid boundaries separated by a distance 2H along the z direction, driven by a constant pressure gradient along the flow direction.We divide the velocity v into two components: P + v = u u , where P u is the Poiseuille flow velocity under an externally applied pressure, and u being the perturbation component, both with the Navier slip boundary condition at the fluid-solid interfaces.If we express length in units of H, mass density in units of the fluid mass density  , time in units of H/ 0 U , where 0 U denotes the maximum velocity of the Poiseuille flow [2,3] at the symmetry plane z=0 under the non-slip boundary condition, then Re= 0 / HU  , where  denotes the shear viscosity.The Navier boundary condition [2,4,14-16] has the form: The associated flow rate P Q per unit area in the yz plane and the associated pressure gradient are given by 0 0 0 12 4 , , 0 3 Re where 0 P denotes the applied pressure that gives rise to the Poiseuille flow.Equations (3a, 3b) represent the Darcy's law [5][6][7] in the simplest channel geometry.It expresses a steady channel flow resulting from the externally applied pressure drop, force-balanced by the dissipative viscous force inherent in the flow profile.It is noted that ( ) 0 PP = uu , hence theoretically the Poiseuille flow can persist to any Reynolds number.However, this cannot be the case because it has been shown that the laminar flow becomes unstable and can blow up with an infinitesimal perturbation, such as those caused by thermal fluctuations, beyond a critical Reynolds number Rec [1,2,[8][9][10].Since physically the channel flow is constrained by a finite energy input, hence what happens beyond Rec is an intriguing question that is the focus of this study.
By formulating the linearized NS equation into an eigenvalue problem as shown below, we have obtained the analytical solution of the complete set of eigenfunctions for the 3D channel flow with the Navier boundary condition.By expanding u in terms of the HMs [17][18][19][20], with P u being the zeroth order term (input) in the expansion, the timedependent expansion coefficients become the variables to be solved.That is, Eq. (1a), with the constraint of Eq. (1b), is reduced to an infinitely coupled autonomous system of ODEs for the expansion coefficients, in which time t is the only independent variable, with Re and s l being the two input parameters.
The use of analytical HMs as the basis function has the advantage of satisfying both the Navier boundary condition and the incompressibility condition, in contrast to the traditional basis functions such as the Fourier basis or the finite-element basis.Hence from the numerical point of view a small number of HMs can already attain high accuracy in the evaluation of the critical Reynolds number, without invoking the Orr-Sommerfeld (OR) equation [4,8].This point is shown below.
In the system of coupled ODEs, the nonlinear term, ( )  uu , takes the form of a third-rank tensor that couples pairs of the expansion coefficients to effect the time variation on the third.Without the nonlinear term, the right hand side of the equation represents a matrix acting on the coefficient vector.By solving for the eigenvalues of the matrix for a finite set of HMs, we found the critical Rec, accompanied by a critical wavevector   , as the lowest Re value at which the real part of one of the eigenvalues first changes sign, signaling the onset of instability.Our Rec is accurate to five significant digits when compared with the standard result of Orszag [8] for the non-slip boundary condition, without invoking the OR equation.The effect of a finite slip length on Rec was also obtained and shown to agree with the prior results [21] on its variation with increasing s l .
Since each HM represents an independent degree of freedom, in a thermal bath each can be excited with /2 B kT of energy, where B k denotes Boltzmann's constant and T is temperature [22][23][24].We examined the Poiseuille flow attendant with the excitation of the HMs by time-evolving the coupled autonomous set of ODEs.It was found that below Rec the excited HMs decay rapidly [2,8,10], and the Poiseuille flow is stable.However, above Re c the picture changes dramatically.Under the problem setup where we vary Re by changing  while simultaneously keeping 0 U constant (by stipulating the applied pressure gradient to vary as 1/Re), it is found that as we numerically time-evolve the autonomous ODE set, the net flow rate initially decreases but eventually settles at a fluctuating equilibrium state that is statistically stationary.
In what follows, we first present the analytical solution of the HMs for the 3D channel flow, in conjunction with a discussion on their properties.This is followed by the evaluation of the critical state and its parameter values, as well as the effect of a finite slip length.By expanding the perturbation velocity field u in terms of the HMs, we reduce linearized NS equation to a system of autonomous ODEs for the expansion coefficients in which t is the only independent variable.Results of the numerical time evolution are presented, with the initial states statistically sampled from a Gaussian distribution for simulating thermal excitations.A new fluctuating equilibrium state is revealed that exhibits a reduced net flow rate.In particular, interesting features of the flow profiles in the equilibrium state are presented through the lens of force balance.We give a brief discussion on the potential experimental verification of the predictions at the end of the manuscript, and conclude with a brief recapitulation.Details on the derivation of the analytical solution are given in the Supplementary Materials (SM), so as not to detract from the main line of exposition.

II. Formulation of the eigenvalue problem and its analytical solution
The HMs are the eigenfunctions of the linearized version of the incompressible NS equation [17,18] with the Navier boundary condition at the solid/liquid interface.They are given by [2,17,18,25]: where   is the component of the velocity field normal to the boundary, and   is its tangential component.In Eq. ( 4) we have used the condition that the time variation of u is given by  − , i.e., the eigenvalue denotes the exponential decay rate of the eigenfunction.The HMs are transient in nature; they cannot be sustained without the input of external energy.Without the loss of generality, periodic boundary condition is imposed on u and p in the  plane.Depending on the spatial dependence of the eigenfunctions, we classify the HMs into 3 disjoint groups: 1D HMs, 2D HMs and 3D HMs.In each group, according to their spatial symmetry property of the velocity components   ,   with respect to the  = 0 plane, they are further classified into the symmetric HMs group and the antisymmetric HMs group.Among these, the 2D and 1D HMs can be further grouped into two different branches, for flow in either the  plane (denoted the  branch) or the  plane (denoted the  branch).The 1D and 2D HMs have already been introduced in previous publications [17,18].Below we focus on the 3D HMs.We emphasize that compared to the NS equation constrained in the 2D spatial domain where vortex stretching is absent and only diffusive exchange of angular momentum exists, in 3D the vortex stretching plays an important role [26,27].Hence the 3D HMs are crucial for evaluating all possible manifestations of the velocity and vortices patterns.
For the 3D HMs, u is periodic along both the ,  directions, characterized by wavevectors  and , respectively.Here  ≠ 0,  ≠ 0 are real, continuous variables, implying infinite periodicity in the absence of a length scale defined either by the computational domain, or by a naturally-occurring length scale.In anticipation of later developments, we note that the emergence of a critical Reynolds number is accompanied by the appearance of a length scale in the  plane, which can be used as the inverse unit for the discretization of m, k in the numerical calculations beyond the critical state.

The general expressions of eigenfunctions
u can be formulated as: (5a) They are to be considered in conjunction with the incompressiblity constraints as well as the Navier slip boundary conditions.In component form, the boundary conditions for u are .
( 1) 0 As HMs must display the translational symmetry in the xy plane,  1 ,  2 ,  3 ,  4 are the constants for expressing the four independent degrees of freedom for adjusting solution's phase in the lateral plane.
For the  1 branch, ( ) Here , ,  are the wavevectors characterizing the eigenfunctions along the , ,  axis, respectively.For each given pair of wavevector (, ), we can obtain a countably infinite number of {  ,  ∈  + } from the following corresponding dispersion relation: The symmetric HMs can be expressed similarly as For the  1 branch, The corresponding dispersion relation is given by In Eqs.(6,7),  represents the wavevector in the xy plane, whereas  is the wavevector along the z direction, which is also the eigenvector of the HMs, to be solved from the dispersion relation.It is to be noted that the eigenvalue , i.e., the decay time of the HM, is always bundled together with Re.Therefore the change in the value of Re can be completely compensated by the decay time of the HMs, without affecting the spatial configurations of the HMs.
The complete HM expression has the form ( ) , where  is a real positive number ensured by the self-adjoint property of the NS equation [29,30].Eigenfunction's time variation is important because the force balance of HMs is between the viscous force and the inertial acceleration, as well as the pressure field associated with each HM's velocity pattern.In contrast, for the Poiseuille flow the force balance is between the viscous force and the externally applied pressure gradient.All the 3D HMs satisfy both the incompressibility condition and the boundary conditions, as well as the linearized NS equation.These are explicitly verified in SM Part I [28].As the basis set for expanding the channel flow velocity, the HMs are consistent with the channel geometry and boundary condition, in contrast to the Fourier basis [31].The latter is more suitable for the infinite system with no boundaries.
For completeness, below we also write down the 1D and 2D HMs.For details, see SM Parts I and II [28].The solution method for  is detailed in SM Part III [28].
The 1D antisymmetric HMs have the form: The 1D symmetric HMs have the form: The 2D antisymmetric HMs have the form: The 2D symmetric HMs have the form  In Fig. 1, we graphically illustrate the 3D, 2D, and 1D HMs with Re=10 4 .Among all the 1D, 2D and 3D HMs, the 1D symmetric HMs are special because they display center of mass motion, which is crucial for giving rise to the change in the net channel flow rate beyond the critical state (Rec, mc).In contrast to the isotropic infinite systems in which the center of mass motion cannot be allowed due to the absence of external force, there can be two external forces in the channel flow: the frictional reaction force from the boundaries, and the force from the externally applied pressure gradient.According to Newton's law, the external forces can alter the center of mass momentum in the channel flow.Beyond the critical state, nonlinear couplings between the HMs with zero net momentum, and the 1D symmetric modes with net flow, are shown to generate boundary reaction force in a direction opposite to the pressure force, leading to a decrease in the net flow.An account of the energy conservation and force balance under this scenario is given in the last section.
The HMs completely separate the NS equations into the spatial and temporal components.They can describe all possible fluid particles' motions in the spatial domain [17,18,24] from the continuum perspective.Each HM satisfies both the linearized NS equation and the incompressibility condition, and the relevant pressure field is coupled to the velocity field through the pressure Laplace (Poisson) equation [32], with the boundary condition given by It is shown below that the pressure gradient term can be eliminated when the NS equation is reduced to a coupled set of time-varying autonomous ODEs, by expanding the velocity in terms of the HMs.Hence the pressure gradient term is not involved either in time-evolving the autonomous ODEs (SM Part IV [28]), or in analyzing the critical state (Rec, mc).
In order to simplify the classification scheme, below we introduce a unified expression for the HMs that includes all the cases stated above.
Here  is the Kronecker delta function and , , ,   ,   are the 5 hyper-indices used to denote different branches of the HMs.Values of , , ,   ,   ∈ {0 or 1}.Here  = 0 denotes the antisymmetric HMs and  = 1 denotes the symmetric HMs;  = 0 denotes the  branch and  = 1 denotes the  branch;  = 0 denotes the 1D HMs and  = 1 denotes the 2D or 3D HMs depending on the values of (, ) tuple;   = 0 denotes the  () phase term of   and   = 1 denoting the  () phase term of   ; and   = 0 denotes the () phase term of   , and   = 1 denotes the  () phase term of   .
As shown in Fig. 1, except for the 1D symmetric HMs, the remaining HMs have zero net momentum and zero net angular momentum within each periodic domain, mostly comprising clockwise and counterclockwise pairs of vortices.When used in calculations, every HM, denoted below as  u , is normalized so that its norm is 1.Normalization is over the volume specified by , where A denotes the area in the xy plane covered by the volume integration.12) is a first order system of coupled ODEs.They can be written in the matrix form as

III. Velocity expansion and the reduction of the Navier-Stokes equation to an autonomous system of ODEs
Here represents the coefficients matrix of the linear terms and represents the 3 rd order tensor term.By utilizing the incompressibility property and periodicity of the HMs in the xy plane, the pressure gradient term can be eliminated in Eq. (12).For details, see SM Part IV [28].The NS equation is thereby transformed formally into a system of coupled ODEs: The two operators , are only functions of Re and slip length   .It should be noted that the Poiseuille flow enters only in the matrix elements of the linear term, , and is absent in .Theoretically we should be able to obtain all the channel fluid dynamics information by time-evolving the above set of coupled ODEs.Here we choose the explicit Runge-Kutta scheme of 4 th order accuracy [33,34] to numerically compute the time evolution trajectories of Eq. ( 13).Details of the numerical scheme are presented in SM Part V [28].By choosing an appropriate time step Δ, the numerical errors can be controlled.
For the linear operator , the pressure corresponding to each  u can be recovered by solving a pressure-Poisson equation with the boundary condition In anticipation of later developments, we note that the exponential decay of each HM as a function of time is implicitly incorporated in Eq. ( 13) through the matrix/tensor elements that contain the (viscous dissipation of its) spatial velocity pattern of the HMs.The sustained presence of the HMs in the flow pattern, as will be shown to be the case beyond the critical state, can only arise when there is sustaining work done by the externally supplied pressure.A detailed consideration of the energy conservation law is given below.

IV. Evaluation of the critical state parameters
The stability of the Poiseuille flow under the nonslip boundary condition is characterized by the existence of a critical value Rec [8][9][10], so that when Re<Rec the Poiseuille flow is stable under small perturbations; while the perturbed velocity field can increase exponentially for Re>Rec, accompanied by the emergence of a critical wavevector.
By considering infinitesimal perturbed velocity field ( ) ,t ur in the linear regime, Eq. (13b) may be expressed as Equation ( 14) represents a linear system whose dynamic stability is determined by the maximum real part of the linear operator spectrum.The critical (Rec, mc) corresponds to the value of (Re, m) when one of 's eigenvalues' real part first crosses 0. The perturbations decay exponentially below the critical state and increase exponentially otherwise.In the search for the critical state parameter values, we have fixed k=0, since from the literature [35] it is already that the most unstable mode can be found in the 2D case with k=0.This point is further elaborated below.We have numerically calculated the linear operator's spectrum from its matrix formulation  as given in Eq. (13b).Here we note that the  matrix has a block-diagonal structure, in which each 2N × 2 block corresponds to a unique value of m (k=0), with N being the number of 's used in the search.By fixing the number of 's at 400, with ∆ = 0.0001 as the step size in the m scan, we focused on the subspace spanned by the finite number of HMs: Bisection method was employed to locate the precise value of Rec.By denoting the maximum real part of the eigenvalues of  as Real() and the imaginary part as Imag(), we scanned the (Re, m) values to determine the first sign change of Real(), and denote the corresponding (Rec, mc) as the critical state.We compare our result with that of Orszag [8]-Re c = 5772.22,  = 1.02056.Our formalism yields Re c = 5772.2and   = 1.0203.Excellent agreement is seen.It should be noted that Imag(  ) = 0.269 also agrees well with the Orszag result, and represents an oscillation frequency of the critical state.A video illustration of this oscillating critical state is shown in SM Part IV [28].
For the 2-dimensional channel flow, Orr-Sommerfeld stability equation [8,10,35] utilizes the nonslip boundary conditions to fix the normal derivative   /=0.In this way, the perturbed NS equation can be compressed into a 1D, 4 th ODE.The aim is to locate the value of Re when Real() = 0. Furthermore, if the eigen wavevector  ≠ 0, then under the nonslip boundary condition the Squire transformation [35] shows that any 3D channel flow's perturbation is equivalent to a specific 2D channel flow perturbation.Hence the task of obtaining the 3D channel flow's critical state is transformed into solving for the (Re c ,   ) of its corresponding 2D flow.By specifying a given input  ≠ 0, the critical   as well as its corresponding Re c can be uniquely determined through the Squire transformation.In Fig. 2 we plot the Re c as a function of   , .The conventional definition of Re c corresponds to the minimum value of the curve at k=0,   = 1.02, which is the most unstable point and labeled as a red dot in Fig. 2. As anticipated previously,   =1.02 defines a natural (inverse) length scale, which will be used in our discretization of the wavevectors in the xy plane for numerical evaluation of the channel flow beyond the critical state.The Orr-Sommerfeld stability equation generally does not apply to the case when   ≠ 0 , but our formalism can easily deal with the more generalized case.No extra modifications are needed.In Fig. 3, the evaluated (Re c ,   ) are plotted as a function of slip length   ranging from 0 to 14× 10 −4 .We see that there exists a minimum   = 0.006 at which the flow has the smallest Re c .It describes an interesting behavior: Increase in the slip length will first lead the system to be less stable, and only after the minimal point is passed will the system become more stable with increasing slip.However, this behavior is conditioned on allowing the Poiseuille flow rate to increase with   as given by Eq. (3b).If we switch to the condition of constant flow rate independent of   , then the variation of Re c with a finite slip length would revert to a monotonic behavior, in agreement with that shown in [21].
An aspect not previously observed is the actual configuration and dynamics of the critical state.Since in our case the critical state configuration is a by-product of the critical eigenvalue of the linear operator , its configuration can be easily constructed from the set . It turns out that while Real(  ) = 0, its imaginary part, Imag(  ) = 0.269 (for the nonslip case), represents an oscillation frequency for the critical state.In SM Part IV [28], we give a video illustration of the oscillating critical state that comprises a vortex and antivortex pair that oscillates at a time periodicity of 2/Imag(  ) = 23.3.Obviously, this oscillation behavior arises from the participation of the Poiseuille flow as the energy source, whose velocity profile has explicitly participated in the matrix elements of the linear operator as shown in Eq. (13b).

V. Channel flow beyond the critical state A. Initial state sampling
Thermal fluctuations constitute an inherent part of the physical system at equilibrium.They are manifest as a source of noise arising from the thermal bath.Without thermal excitations, the Poiseuille flow can be absolutely stable to any Reynolds number.In previous works [17,18], it was shown that in terms of the complete set of HMs, the fluctuation-dissipation theorem [23,24,36,37] may be expressed in a simple expression involving the average of the inverse of the HMs' eigenvalues.By considering the fluid molecular motions as superpositions of the HMs, the thermal fluctuations in fluid can be represented from the continuum perspective as canonical ensembles of the projection coefficients ()with initial Gaussian distributions following the Heisenberg picture of phase space averages.Moreover, since each HM represents one degree of freedom [24,17,18], from the equipartition theorem each HM should have the same thermal energy [24].Hence the initial probability density function of the projection coefficients () should obey the Gaussian distribution where 2  may be regarded as the thermal energy.Given the initial state with coefficients distributed in accordance to (), the ensemble averaged perturbed net flow ( ) Qt can be computed as ( ) ( ) ( ) 1 symmetric HMs flow rate where  1 is the vector of 1D symmetric HMs' flow rates.
To simulate the averaged flow rate under thermal excitations, we first prepare a set of independent and normally distributed random variables  ≔ { 1 , … ,   } with variance  2 .Assuming that a total of  trajectories are used to calculate the ensemble average, we repeat the above process  times to obtain a set of random variable vector sets  0 = {  1 ,  2 , … ,   }.Then by performing the time evolution stipulated by Eq. ( 13) with each element in  0 set to be the initial state, we obtain  independent trajectories { c 1 (), c 2 (), … , c  ()} .The ensemble-averaged flow rate at any time  may be evaluated as In our numerical time evolution of Eq. (13b), we used  = 10 evolution trajectories, where the initial state of each trajectory comprises the 2D antisymmetric modes along the x and y branches with discretized values of m and k in multiples of   × 2  , where i=0,1, 2, 3, 4, 5. Multiples of 2 were chosen because they represent the maximum couplings mediated by the tensorial term [38].The amplitudes of the modes were generated by the Gaussian random process described above, with  2 = 0.04.Only 2D antisymmetric modes were chosen as the initial states since the symmetric modes have not been found to be unstable, a fact that has not yet been rigorously proved.In our numerical evaluations, Re was fixed at 10 4 , and the computational discretization in the x, y directions was set at ∆ =   , so that the length L and width W of our computational domain in the xy plane are fixed at It should be noted that even though the initial states were all 2D antisymmetric modes, as time evolved, the 3D modes appeared.Hence the 3D HMs are absolutely needed in the numerical experiment.In the time evolution as stipulated by Eq. ( 13), a total of 6596 HMs were involved.
Poiseuille flow below the critical state is always stable with all excited HMs decaying exponentially as a function of .Hence the ensemble-averaged perturbed net flow

( )
Qt rapidly decays to 0 below the critical state.Above the critical point, the situation changes due to the amplification effect of the linear operator .However, the amplification effect cannot be sustained since the total kinetic energy is continuously being dissipated, while the other perturbation components are also excited through the nonlinear coupling term ( )  uu .In this manner, the exponential increase is suppressed.
What happens when an extremely long time  passes?From the following simple argument we can see that the perturbed flow cannot be zero.Suppose at a time t all the perturbation components tend to 0 above the critical state, then the linear analysis tells us that the thermal excitations will again make those HMs with positive real part of the eigenvalues exponentially increase in magnitude, and the whole process can start over again.It follows that the final equilibrium state should have a nonzero ensemble-averaged perturbed net flow, accompanied by fluctuations.

B. Reduced net flow rate
Physically, this anticipated outcome originates from the excitation of the 1D symmetric modes mediated by the tensorial term, leading to counter flow in a direction opposite to the Poiseuille flow.This is indeed the result as shown in Fig. 4.

FIG. 4.
A plot of the normalized net flow rate, in unit of P Q , as a function of time .The time evolution is evaluated from Eq. ( 13) and Eq. ( 18), with the Reynolds number set at Re=10 4 .The ensemble-averaged net flow rate, with ten independently generated initial states, is seen to decrease as a function of time from the initial Poiseuille flow rate.The yellow curve, indicative of the case for the finite slip length, is seen to approach the equilibrium flow rate slower than that of the nonslip case.
We have adopted the 4 th order explicit Runge-Kutta numerical scheme to evaluate the time evolution of Eq. ( 13).Details of the numerical scheme can be found in SM Part V [28].In Fig. 4 we plot the time evolution of the flow rate variation as evaluated by Eq. ( 13) and (18), at Re=10 4 .The curves were generated by ensemble-averaging ten independent initial states, run in parallel on 4 Amazon servers for about 3 months.The net flow rate shows the deviation from P Q .We should be reminded that even though the Reynolds number used in the evaluation of Eq. ( 13) is above the critical value, the Poiseuille flow, in the absence of perturbed modes, should remain the same.This is because even though the pressure varies inversely as 1/Re, the decrease in viscosity exactly compensates the lowered pressure to maintain the same flow rate.Hence the deviation as shown in Fig. 4 can only be ascribed to the appearance of the perturbed HMs, especially the 1D symmetric modes.The net flow rate is seen to decrease until it reaches a plateau that is ~15% lower than the Poiseuille flow rate.However, not seen from Fig. 4 is the presence of the many HMs with vortices and anti-vortices that accompany this state.There are perceptible fluctuations around the plateau.
An interesting observation is that in accordance to Newton's law, the reduced net flow would require the boundary reaction force to be along the direction of the counter flow, i.e., against the Poiseuille flow.That would be opposite to what one would usually expect.This puzzle is resolved in the context of the counter flow profiles and force balance, presented below.

C. Energy conservation law
With the solution of the velocity expansion coefficient set {  }, we can also verify, at every time instant, the energy conservation for each independent trajectory as a validation of our numerical implementation.The energy conservation law states that   , must be equal to the net work.In Fig. 5 we plot the ensemble-averaged, time-integrated kinetic energy as a function of , while the inset shows the same for one selected evolution trajectory.We see excellent consistency between the two curves as evaluated from two sides of Eq. (19b), indicating energy conservation to be well obeyed by our numerical scheme for every evolution trajectory.
It is clear from Fig. 5 that the general trend in the kinetic energy evolution is along a decreasing path.This trend is due to the fact that most of perturbation modes' energy becomes concentrated in the 1D symmetric modes as t increases, with the rest of the HMs having at least an order of magnitude less energy.Since the effect of the 1D symmetric modes is to decrease the net flow, hence the kinetic energy of the whole system decreases.While the other zero-momentum modes (mostly in the form of coupled vortices and antivortices) have small amplitudes, they provide the energetic link between the 1D symmetric modes and the Poiseuille flow.Otherwise the 1D symmetric modes would decay and disappear.

FIG. 5.
The ensemble-averaged kinetic energy as evaluated from the energy conservation relation, Eq. ( 19), plotted as a function of time .Here the kinetic energy is numerically time-integrated from Eq. ( 19).The blue line indicates the left hand side of Eq. ( 19), and the red dashed line indicates the right hand side.The nonslip boundary condition was used in this case.The same quantity for a single trajectory is plotted in the inset.It is seen that the energy conservation law is very well obeyed in the time evolution, and the kinetic energy decreases as a function of time, until it reaches a plateau.The vertical arrow indicates the time, t=7000, at which the flow profile is plotted in Fig. 6.

D. Net flow profile and the counter flow
Beyond the critical Reynolds number the net flow comprises both the Poiseuille flow and the 1D symmetric modes.In Fig. 6(a) we plot the total velocity profile,   , at t=0 and t=7000 for a comparison.Here we consider only the nonslip boundary condition; the same holds true for what follows.The time instance t=7000 falls in what we denote as the "equilibrium state" time period.It is seen that the two velocity profiles are largely similar.However, the presence of the 1D symmetric modes can be made apparent by subtracting off the t=0 Poiseuille profile from the net flow profile.The result, denoted as the counter flow    , is plotted in Fig. 6

E. Force accounting
The NS equation expresses Newton's law in the form   ⁄  + ( • ∇) = ∇ •  ⃡, where  ⃡ denotes fluid's stress tensor.In view of the findings above, an interesting question arises: Is it possible to have a statistical stationary (or equilibrium) state at Re>Rec?A stationary state would imply the vanishing of the inertial force on the left hand side of the NS equation.At t=0, the Poiseuille flow is obviously stationary, but in the presence of counter flow at t=7000, the vanishing of the inertial force seems very unlikely, even if it is on the statistical basis of taking the time-averaged flow profile.We show below that a detailed force accounting yields a somewhat surprising answer.
The natural starting point is to volume-integrate both sides of the NS equation, and convert the right hand side of the NS equation to a surface integral over the boundaries of the computational domain so that the integrated equation expresses the equality of the total inertial acceleration to the total external forces exerted on the system.The latter comprises the force from the externally applied pressure gradient, plus the frictional reaction force from the solid/liquid interfaces.Since the applied pressure gradient is always given by −2/Re, hence we denote 4 × 2 ×  The fluctuating nature of the curve in Fig. 7(a) raises an obvious question: What if one takes a time average over this (near-equilibrium) time period?The result is shown as the flat dashed line in Fig. 7(a), which is slightly negative but very close to zero, i.e., the time-averaged boundary reaction force caused by the counter flow is nearly zero.In Fig.

F. Experimental possibilities
Our numerical experiment was limited by the available computational resources.Hence it is not possible to test whether the final result would be altered by considering a very large set of initial states, or by considering many more coupling possibilities via the tensorial term.However, experimental verification of the predictions should be possible.Viscosity can be varied by mixing glycerol and water in varying proportions, for example.Two experiments can be carried out by using the same channel but with different applied pressure across the sample, with compensating viscosity values so as to maintain the same 0 U associated with the Poiseuille profile.The value of Re can thus be varied, with one experiment below the critical state and another one above it.In this manner the flow rate of the two experiments can be measured and compared.In another experiment, perhaps more difficult, is to visualize the 2D critical oscillating state, if even in a transient condition.

V. Conclusion
By obtaining the analytical eigenfunctions of the linearized incompressible NS equations in the channel geometry, we reduced the full NS equation to a set of coupled autonomous ODEs with t being the only variable.The critical Reynolds number and the critical wavevector were obtained to five and four significant figures, respectively, together with an explicit configuration of the oscillating critical state.By considering only a limited number of tensorial couplings, and with a small set of initial perturbed HMs, we reveal that beyond the critical state, the net flow rate decreases due to the nonlinear coupling between the 1D symmetric HMs and other HMs with zero momentum.There is a plateau regime for the net flow rate, reduced by ~15% from that in the linear regime, almost independent of the slip length.The 1D symmetric modes are shown to co-exist with the Poiseuille flow and the vortices.However, it remains an open question as to the ultimate net flow rate in the limit of t→ ∞, or Re→ ∞, or both.The resolution of this question may require mathematical analysis of the problem, which is beyond the scope of this work.

where sl
is the slip length.The laminar Poiseuille profile is given by

(
of complete, normalized HMs with subscript  as the labeling index, denotes inner product, and all the spatially-dependent variables are converted into analytical scalar expressions through inner product projections.Equation (

FIG. 2 .
FIG. 2.Plot of Re c as a function of the eigen-wavevectors   , .The conventional definition of Re c =5772.2 is noted to be the minimum point of the discretized curve (i.e., the least stable mode), labeled by a red dot.Here k is an input parameter.

FIG. 3 .
FIG. 3.A plot of the critical wavevector   (in black) and the critical Reynolds number Re c (in blue), both as a function of the slip length   .Here the Poiseuille flow rate is allowed to vary with the slip length as stipulated by Eq. (3b).
the rate of energy input exerted by the external pressure, ( ) d Wt is the rate of viscous dissipation, and the left hand side of Eq. (19a) is the rate of change of the kinetic energy.In SM Part VI[28] we give the derivation of the expressions for of the HMs and the velocity expansion coefficients.This form of the energy conservation law can be expressed as coefficients of the total velocity field v , including the Poiseuille flow velocity.The latter is expanded in terms of the 1D symmetric HMs.Here N=6596 denotes the total number of HMs involved.Derivation of Eq. (19b) is given in SM Part VI[28].The first term on the right hand side represents the energy density input exerted by the external pressure, while the second term represents the energy dissipation rate caused by the viscous stress tensor.The energy conservation law tells us that the change in the total kinetic energy, (b) at three time instances.The    is composed of different 1D symmetric modes, not all of them flowing along the same direction, e.g., at t=7000,    = −0.158cos(1.571)+ 0.01 cos(4.712)+ 0.0075 cos(7.854)− 0.0074 cos(10.996)+ ⋯.Hence one can see that in the vicinity of the solid boundary, there can be opposite slopes of [ different directions of the boundary reaction force at difference time instances.In particular, the blue curve at t=7000 is seen to have a boundary reaction force along the same direction as the counter flow.That is in contrast to the Poiseuille flow, where the boundary reaction force is always opposite to the flow direction.For almost all    (t) prior to reaching the equilibrium state, the flow profile in the vicinity of the boundary is similar to the blue curve in Fig. 6(b).It explains how the initial counter flow can provide the necessary boundary reaction force for decreasing the net flow, in accordance to Newton's law.Below we examine in more detail the question of force balance and the necessary condition for a statistical equilibrium (or stationary) state.

FIG. 6 .
FIG. 6.(a) Flow profile at two time instances, where vx denotes the total velocity.At t=0, the profile is parabolic, corresponding to the Poiseuille flow.At t=7000, the net flow rate is decreased, even though the profile is still very similar.(b) The counter-flow profile    at three time instances, obtained by subtracting off the Poiseuille profile from the total flow profile.The inset shows an enlarged view of the velocity profile in the vicinity of the solid boundary.The counter flow profile can be expressed as a superposition of the 1D symmetric modes, not all flowing along the same direction.It is noted that the boundary reaction force is in the negative x direction at t=7000 (blue curve).At other two time instances, t=7560 (red) and t=7345 (black), the boundary reaction force is either along the positive x direction or nearly zero, respectively.The fluctuating boundary profile as a function of time has a direct bearing on the time-averaged frictional reaction force from the solid boundary arising from    .The blue curve is representative of almost all the counter flow profiles prior to t=7000.

2 Re=
0 as the unit of force (along the x direction) in this accounting exercise, i.e., total force from applied pressure gradient=1.By denoting the counter flow velocity profile as    , it is clear that the frictional boundary reaction force ∆   arising from    must be equal to the inertial force ∆   , since pressure force is completely compensated by the boundary reaction force from    , and the total boundary reaction force is the linear addition of that from    and    .In Fig. 7(a), we plot the inertial force ∆   as a function of time (blue dots), together with the boundary reaction force ∆   caused by    (red line).Complete agreement is seen, since the two represent the two sides of the same equation.

FIG. 7 .
FIG. 7. (a) A plot of the normalized inertial force and the boundary reaction force as a function of time.The time period is chosen to be in the near-equilibrium state, i.e., after t=7000.The blue dots denote the inertial force, and the red curve denotes the boundary reaction force.Complete agreement is seen, since the two represent the two sides of the same equation.The black dashed line represents the time average of the 7(b)  we plot a time-averaged profile of    over the same period.It is seen that at the upper and lower boundaries, [ indicating a close-to-zero frictional reaction force from the boundary!This is a somewhat surprising result from our limited data, but consistent with the finding of a statistical equilibrium state that requires the absence of inertial force in the time-averaged sense.This scenario is conceivable because    comprises different 1D symmetric modes with both +  − flow directions, thereby making it possible to have no net boundary reaction force via cancellations.This point has already been foreshadowed in Fig.6(b), which shows that there can be opposite slopes of [ instances in the vicinity of the boundary, including one that shows [ In summary, our analysis yields a physical picture of the equilibrium state beyond Rec with the following characteristics: (1) fluctuations; (2) co-existence between Poiseuille flow, 1D symmetric modes, and the vortices, with the transient nature of the latter underling the fluctuations; and (3) reduced flow rate from the Poiseuille flow.
the solution denotes the parity of   ,   with respect to the  = 0 plane.However,   must have the opposite parity to that of   ,   as dictated by the incompressibility condition.This can be easily seen by expanding the  dependence of