Variational Nonlinear Optimization in Fluid Dynamics: The Case of a Channel Flow with Superhydrophobic Walls

: Variational optimization has been recently applied to nonlinear systems with many degrees of freedom such as shear ﬂows undergoing transition to turbulence. This technique has unveiled powerful energy growth mechanisms able to produce typical coherent structures currently observed in transition and turbulence. However, it is still not clear the extent to which these nonlinear optimal energy growth mechanisms are robust with respect to external disturbances or wall imperfections. Within this framework, this work aims at investigating how nano-roughnesses such as those of superhydrophobic surfaces affect optimal energy growth mechanisms relying on nonlinearity. Nonlinear optimizations have been carried out in a channel ﬂow with no-slip and slippery boundaries, mimicking the presence of superhydrophobic surfaces. For increasing slip length, the energy threshold for obtaining hairpin-like nonlinear optimal perturbations slightly rises, scaling approximately with Re − 2.36 no matter the slip length. The corresponding energy gain increases with Re with a slope that reduces with the slip length, being almost halved for the largest slip and Reynolds number considered. This suggests a strong effect of boundary slip on the energy growth of these perturbations. While energy is considerably decreased, the shape of the optimal perturbation barely changes, indicating the robustness of optimal perturbations with respect to wall slip.


Introduction
The original view of turbulent transition in shear flows, consisting in the linear amplification of unstable eigenmodes followed by secondary instability and nonlinear mixing, has been recently replaced by new scenarios focusing on receptivity [1], transient growth [2], and coherent flow structures [3]. In fact, even when all eigenmodes are damped, disturbances may be strongly amplified by means a non-normal growth mechanism arising from the constructive interference of non-orthogonal eigenfunctions [4]. As this growth is linear for infinitesimal disturbances, the concept of optimal perturbations was introduced to identify the maximum possible growth that a perturbation to the laminar "base" flow can experience. Optimal disturbances are defined as initial perturbations which yield the largest amplification of their energy over a time/space interval [5,6]. When the Reynolds number is smaller than its critical value (i.e., the value for which the linearized Navier-Stokes operator has unstable eigenmodes), these perturbations may rapidly lead the flow to turbulent transition. Subcritical transition in shear flows has been studied using linearized approaches for computing the perturbations to the laminar solution able to induce the largest possible growth at a given target time, namely, the linear optimal perturbations [7]. In shear flows, such as channels, boundary layers, and plane Couette flow, linear optimal perturbations consist of pairs of counter-rotating streamwise vortices, capable to elicit streamwise streaks by the lift-up effect [8]. A common objection to this transition scenario is that optimal streamwise-invariant initial disturbances can be rarely observed in real flows. In fact, in most practical cases, the flow undergoes transition by receptively selecting and amplifying freestream turbulence perturbations [9,10], or localized disturbances, such as those arising from the presence of roughness elements or gaps on the wall. As an example, Figure 1 provides sketches of three different scenarios of transition triggered by freestream turbulence of increasing intensity: while in low-noise configurations modal instabilities may arise, increasing environmental disturbances trigger bypass transition which relies on the generation of coherent structures such as streamwise streaks or hairpin vortices. Linear optimal perturbations are found to reproduce some of these flow structures very well, as suggested by the comparison of the optimization results with experimental data [6]. This surprising result may be easily explained by considering that, when the largest singular values of the problem are well separated, perturbations having a large projection onto the optimal one would provide a large contribution to the energy amplification [6], leading to a flow dominated by optimal and near-optimal structures even when the flow is excited by environmental, random disturbances. However, many studies have demonstrated that linear optimization is inefficient at estimating thresholds for transition and, for example, the so-called oblique transition mechanism, as well as suboptimal disturbances obtained in a linearized framework [11] succeed in triggering transition at a much lower initial disturbance level than linear optimal ones [12,13]. In fact, linear optimizations do not provide any relevant information about the minimal energy for which a given perturbation may actually experience transition, nor give hints about the mutual interaction between linear mechanisms such as, for instance, transient growth of streaks and their secondary instability. Indeed, in the linearized system, although the energy can transiently grow, transition to turbulence is not possible. Thus, one can conclude that non-normality of the linearized operator is not the only responsible for subcritical transition, nonlinearity being an essential ingredient that one cannot neglect. Some efforts on including both non-normal and nonlinear effects into simple flow models [14] or in a simplified set of equations [15][16][17] have been done before this decade. However, only in recent years, the increase of available computational resources has pushed the development of numerical strategies allowing to optimize nonlinear problems constituted by a large number of degrees of freedom such as the Navier-Stokes equations for fluid flows. The first nonlinear optimization in shear flows have been conducted, almost at the same time, by two different research groups, for pipe [18] and boundary-layer [19] flows. Both studies found nonlinear optimal perturbations very different from linear ones [20,21]. These optimal perturbations were found able to optimally connect different linear energy growth mechanisms by means of nonlinearity. Similar optimal perturbations have been found for different objective functions [22] and different flows [23][24][25][26][27], or regimes [28,29]. The relevance of these optimal perturbations for the determination of minimal energy thresholds for transition has been also assessed for different shear flows [30][31][32][33]. Noticeably, the nonlinear variational procedure is well adaptable to passive or active flow control [23,[34][35][36], to the computation of invariant flow states [37] and orbits connecting them [38]. A thorough review of nonlinear optimization of perturbations in shear flows can be found in [39,40].
Despite all this body of work, it is still not clear to what extent these nonlinear optimal energy growth mechanisms are robust with respect to external disturbances or imperfections such as the presence of environmental noise, freestream turbulence, isolated, or distributed wall roughness, which are potentially present in all practical applications of these shear flows. In particular, distributed surface roughness is of strong interest in the context of drag reduction by passive control means, since the discovery of biomimetic surfaces allowing to reduce friction at the wall. An example of these bioinspired surfaces are Superhydrophobic Surfaces (SHS) [41], which are composed by a hierarchical nanostructure able to trap air within the roughness elements (see sketch in Figure 2). When a liquid flows over such a nanostructured surface, the presence of gas-liquid interfaces reduces the wetting by limiting direct contact of the liquid with the solid substrate. Thus, when the nonwetted [42] state is maintained, the overlying liquid flow remains only partially in contact with the solid substrate, being partially sustained by the underlying mattress of trapped gas bubbles (plastron). The overall effect corresponds to a lubrication the overlying liquid flow. In recent years, a number of experimental works have studied the effect of SHS surfaces on wall-bounded flows, at first in laminar microchannels [43][44][45] and more recently in fully turbulent shear flows [46][47][48][49]. For laminar flows, Ou et al. [43] experimentally demonstrated the potential of these surfaces for drag reduction, finding resulting slip velocities greater than 60% of the average velocity measured at the wall. Further studies on wetting stability [50] and drag reduction [51][52][53][54] on laminar flows have reported a drag reduction which increases approximately linearly with the texture size, up to the point where wetting transition occurs [55], reaching a maximum of 40%. Numerical studies have been pushed by the work of Ybert et al. [56], who have shown that these surfaces can be accurately modeled as a partial wall slip. In the following years, researchers have reported even larger drag reduction [46,57,58], reaching up to 90% of the initial drag [49]. Even in turbulent flow conditions [59][60][61], superhydrophobic surfaces have demonstrated their potential in reducing turbulent drag for a number of flow configurations [62]. On the other hand, very few studies have by now focused on the effect of SHS on the transitional dynamics. The first study tackling this issue [63] considered the K-type transition scenario on a channel flow with slippery boundary. Transition to turbulence has been found to be considerably delayed for sufficiently large (but still practically viable) slip lengths. This pioneer study has been followed, some years later, by two studies in which different transition scenarios have been analyzed in the presence of slippery boundaries [64] and feature-resolved SHS [65]. These studies reported very different behaviors of the transitional channel flow in the presence of SHS, depending on the particular type of perturbations used for triggering turbulence. This suggests that the effect of wall slip on different energy growth mechanisms is not trivial at all. Some recent studies have assessed the effect of boundary slip on the optimal growth of perturbations under linearized conditions [66], finding little effect on the linear energy amplification, while linear stability [63,66,67] appears to be much more affected. At the same time, a recent study [68] has assessed the influence of slip at the wall on nonlinear traveling waves solutions of the Navier-Stokes equations. Depending on their structure, the transition to turbulence of these nonlinear solutions can be delayed or advanced. Thus, it appears that these engineered surfaces may considerably affect strong energy growth mechanisms relying mostly on nonlinearity. However, their potential effect on the optimal growth of nonlinear perturbations has never been investigated. This work aims at filling this gap, trying to unveil the influence of SHS on the shape and amplification of nonlinear optimal perturbations in a channel flow. This analysis bears an enormous potential interest for passive flow control, as superhydrophobic surfaces that succeed at hindering the laminar-turbulent transition process may potentially lead to a reduction of friction losses up to 50%, ensuring a huge energy saving in the flow transport process. The present work aims at presenting the variational mathematical framework used for nonlinear optimization and its application to the study of subcritical transition in channel flows with no-slip or superhydrophobic walls. Section 2 provides the mathematical and computational framework, with Section 2.1 describing the optimization procedure and the related gradient algorithms (Sections 2.1.1 and 2.1.2), whereas a description of the flow setting is provided in Section 2.2. Section 3 provides results of the applications of the presented algorithm to the case of the channel flow with no-slip or superhydrophobic walls. In particular, Section 3.1 presents in detail the no-slip case, Section 3.2 reviews recent results about transition to turbulence in channel flow with superhydrophobic walls, while Section 3.3 describes the optimization results for the channel with superhydrophobic walls. Section 4 concludes the paper.

Mathematical and Computational Framework
The dynamics of an incompressible Newtonian fluid flowing in a channel is governed by the Navier-Stokes equations: where U = (U(x, t), V(x, t), W(x, t)) T is the velocity field and P(x, t) is the pressure field. The Reynolds number is defined as Re = U c h/ν, the reference velocity U c being the centerline one, h being half the height of the channel and ν the kinematic viscosity of the fluid. For a given value of Re, the flow rate is kept constant no matter the wall slip introduced. The reference frame x = (x, y, z) T is chosen so that x is the streamwise, y the wall-normal and z the spanwise direction.

Variational Optimization
Let us define a base flow Q = (U, P) T which is a steady solution of the Navier-Stokes equations. A perturbation q = (u, p) T superposed at t = 0 to the base flow can increase or decrease its energy, depending on linear stability, non-normality, and nonlinearity of the equations. The aim of the nonlinear optimization used here is to maximize, at a certain time T (the target time), the energy of the perturbation defined as E(T) = u(T), u(T) = ||u(T)|| 2 2 , given the scalar product V being the volume of the considered computational domain, and u, v two generic velocity fields in that domain. The objective function of the optimization procedure, E(T), is subject to partial differential constraints such as the initial energy of the perturbations, E 0 , and the incompressible Navier-Stokes (NS) equations which should be now written in a perturbative formulation with respect to the base flow: The mentioned constraints are imposed within the optimization procedure via scalar product with the Lagrange multipliers (or adjoint variables) (u † , p † , E † ), yielding the augmented functional L u, p, u † , p † , u(T), u(0), E † : Integrating by parts and setting to zero the first variation of L with respect to (u(t ∈ (0, T)), p) provides the adjoint equations Whereas setting to zero the first variation of L with respect to u(T) yields the compatibility conditions δL δu(T) = u(T) − u † (T) = 0; which relate the direct to the adjoint variables at target time. Furthermore, the first variation of L with respect to u(0) yields the gradient which should be iteratively nullified in order to reach the maximum of the constrained objective function. This is attained using an iterative procedure which is described in the following.

2.
The direct NS Equations (4) are numerically integrated from t = 0 to t = T.

4.
Starting from these values of the adjoint variables at t = T, the adjoint NS Equations (6) are integrated backward in time.

5.
At t = 0, the initial direct state is updated in a direction determined by the gradient (8) using different gradient-based methods, which are detailed below.

6.
The objective function E(T) and its gradient are evaluated to determine the convergence of the algorithm. If convergence is achieved, the loop is stopped; otherwise, the procedure is continued from step 2.
The whole direct-adjoint looping procedure is schematized in Figure 3. For checking the convergence of the algorithm, the residual = (E(T) n − E(T) n−1 )/E(T) n is measured at each iteration. When its value is lower than the given threshold 10 −7 , convergence is achieved. The main computational drawback of this direct-adjoint looping procedure is the requirement of storing all direct variables at each time step during the direct NS integration, since they are necessary in the backwards time integration (see the advective terms of Equation (6)). This can lead to storage issues, which may be avoided by using lighter storage alternatives such as the receding-horizon algorithm [35,69]. As mentioned above, at each iteration, once the gradient has been computed based on the time integration of the direct and adjoint problems, different algorithms based on the gradient can be used for updating the initial perturbation allowing to initialize the direct problem in the next iteration. In this work, most of the optimizations have been carried out using the gradient-rotation method [25,70], which avoids direct determination of E † by projecting the perturbation on a hypersphere of energy E 0 . However, simpler methods, such as the steepest ascent [18] or the conjugate gradient [11] methods, have been used with success, although they computationally perform worst than the gradientrotation previously discussed and very sensitive to the choice of E † . Details on two possible algorithms that can be used for ensuring convergence of the direct-adjoint looping procedure, are given below.

The Gradient Rotation Algorithm
Once the gradient δL/δu(0) is computed making use of the results of the direct and adjoint time integrations, the initial solution for the next iteration should be computed. In particular, when the gradient rotation algorithm is chosen, the initial condition at the (n + 1) th iteration is updated using the gradient ortho-normalized with respect to u(0). This is done in order to allow a robust exploration of the state space and avoid convergence in local maxima. Thus, the update of the initial condition at the (n + 1) th iteration is performed in the following way, where the angle α can have any value in the range 0 < α < 2π and the second term on the right is the gradient in Equation (8) after ortho-normalization with respect to u(0). In the following, we will explain in detail how the ortho-normalization of the gradient is carried out, in order to fulfill the imposed constraints. As the initial condition at iteration (n + 1) must have initial energy equal to E 0 , we have to ensure that Injecting the expression in Equation (9) in the constraint above, we obtain Notice that the first term of this equation is equal to E(0)cos 2 α, whereas the second term of this equation is zero, as the ortho-normal gradient is, for definition, perpendicular to u(0) n . Thus, to ensure that the sum of the first and third terms of Equation (11) are indeed equal to E 0 , the energy of the orthonormalized gradient should be equal to E 0 as well. Thus, we first orthogonalize the gradient with respect to u(0) obtaining the orthogonal gradient (δL/δu(0)) ⊥ and then we normalize its energy, simply scaling it for a constant such that where G 0 represents the energy density of the perpendicular gradient. A sketch of this method is shown in Figure 4. Thus, when gradient rotation is chosen, the step 6 of the direct-adjoint looping is split in the following substeps.
(a) Orthogonalize the gradient computed at the present iteration with respect to u(0) n (b) Normalize the orthogonal gradient as in Equation (12) (c) Update the initial perturbation for the iteration (n + 1) using Equation (9) (d) Evaluate convergence evaluating the residual . If < 10 −7 , the algorithm is stopped, u(0) (n+1) being the desired optimal solution. Otherwise, the whole procedure continues from step (2).

The Conjugate Gradient Rotation Method
Although the gradient rotation method is very efficient in leading to the optimal solution, simpler methods such as the steepest ascent [18] or the conjugate gradient [11] methods have been used with success, although they computationally perform worst than the gradient rotation previously discussed and are very sensitive to the choice of E † . However, the performance of these algorithms can be considerably improved, when they are combined with the previously described rotation method. For instance, concerning the conjugate gradient rotation, the direction of the update is computed using the gradient at iterations n and n − 1, providing the following update, where is obtained computing the conjugate gradient and orthonormalizing it with respect to u(0) n . In the expressions above, the superscript C indicates the conjugate gradient and, using the Polack-Ribière formula, Notice that, as the update (13) has almost the same form as in (9), the orthogonal conjugate gradient can be easily normalized to the energy E 0 using an equation equivalent to (12), obtained by injecting the expression for the updated initial perturbation (13) in the constraint (10). Thus, for the present algorithm, the step 6 of the direct-adjoint looping is split in the following substeps.
(a) Compute the conjugate gradient using Equation (14) (b) Orthogonalize the conjugate gradient computed at the present iteration with respect to u(0) n (c) Normalize the orthogonal gradient as in Equation (12) (d) Update the initial perturbation for the iteration (n + 1) using Equation (13) (e) Evaluate convergence evaluating the residual . If < 10 −7 , the algorithm is stopped, u(0) (n+1) being the desired optimal solution. Otherwise, the whole procedure continues from step (2).
Notice that the first update must be carried out with the gradient rotation method, which is based on the gradient at the iteration n only.
Both gradient rotation and conjugate-gradient rotation methods have been used for optimizations in the channel flow, leading to virtually the same results with comparable efficiency. For most of the optimizations shown in this paper the gradient rotation method has been used, since it reduces the storage (the gradient at n − 1 does not need to be stored or orthonormalized), with a comparable numerical cost with respect to the conjugategradient one. It is also worth to remark that, for ensuring convergence of these methods to the maximum of the objective function, it is essential to carefully choose the value of α used for the update of the solution, and even change it during the optimization when the residual is not decreasing enough. In particular, in the presented optimizations we start with an arbitrary value of α and check the value of the objective function at the n th -iteration: if this is smaller than that achieved at the n − 1 th iteration, we reduce the value of α and repeat the update of the solution at the n − 1 th iteration.

Modeling of the Superhydrophobic Surfaces
In the last decades, a number of theoretical and numerical studies in laminar [43,71,72] and turbulent flows [46,48,73] have established that microtextured superhydrophobic patterns in wetting-stable conditions can be effectively modeled via an equivalent spatially homogeneous Robin (or Navier) boundary condition: where [U, W] and ∂[U, W]/∂y are the slip velocity and the shear rate in the wall normal direction y at the boundary (y = ±1), respectively, and L s denotes the slip length, as sketched in Figure 2. Here, we use of a single slip length L s since we will restrain our analysis to the case of isotropic surfaces, homogeneous in the wall-parallel directions. This homogenized approach [74] allows for a huge reduction of the computational cost of the overlying flow as a the computational grid can be up to 12 times rougher than that required in an equivalent computation using heterogeneous slip/no-slip boundary conditions over the micro-roughnesses, as discussed by Seo and Mani [73]. Within this homogenized framework, the base flow assumes the form as found by Philip [75]. Notice that the centerline velocity, classically used as a reference velocity in channel flows, changes with the slip length. In this work, following Min and Kim [76], slip length values of L s = 0.005, 0.01, 0.02, 0.05 have been considered. Many studies have shown that, in most flow conditions, L s is linearly dependent on the texture size L, up to a threshold value for which wettingtransition [77] is observed, and the homogenized approximation is no longer valid [47,49,78]. Thus, to verify that these slip length values can be realistically obtained by wetting-stable superhydrophobic surfaces, we follow Seo and Mani [73] evaluating the quantity where L + = LRe τ and L + s = L s Re τ , respectively, Re τ = u τ H/ν being the friction Reynolds number, and φ s is the solid fraction, quantifying the fraction of solid contact within the SHS roughness pattern, typically spanning from 0.1 to 0.3 [47,79]. Setting φ s = 0.25. and Re τ ≈ 200 (which is the value of the friction Reynolds number obtained in the considered flow cases after transition occurs), we find L + ≈ 19, which is within the wettingstable conditions reported by Seo et al. [78]. Moreover, an L + of O(10) corresponds to a micro-roughness size comparable with those successfully employed in experimental campaigns [46,47,49], and is at the same time sufficiently small to avoid any direct interaction with the O(100) coherent structures characteristic of a turbulent flow [80].

Results
In this section, we present the results of the nonlinear optimization procedure described above, with particular focus to the plane Poiseuille flow. At first, some results of nonlinear optimizations in channel flow with no-slip walls are reviewed, which are then used for selecting meaningful parameters for the successive optimizations in the case of the channel flow with superhydrophobic surfaces. In both cases, the optimization has been performed on perturbations of the laminar base flow, as discussed below. In all the presented cases, at the wall, no-slip boundary conditions are imposed for the wall-normal component of the velocity, whereas for the streamwise and spanwise components, Robin boundary conditions are used, with L s = 0 in the no-slip cases. Periodicity is prescribed in the streamwise and spanwise directions. Computations are performed using the spectralelement code NEK5000 [81], with Legendre polynomial reconstruction of degree 7 and second-order accurate Runge-Kutta time integration [82]. Depending on the selected target time, 40 to 80 direct-adjoint iterations are needed for reaching convergence for one set of parameters, each optimization needing 60, 000 to 800, 000 CPU hours on an IBM cluster Intel ES 4650 on a minimum of 128 processors. The numerical setting has been validated by comparing the optimal solutions and the energy gain with those obtained in [25] for the case with Re = 4000, E 0 = 2 × 10 −6 . The optimal perturbation is almost identical to the reference one, and the energy gain is retrieved with a deviation of ≈0.6% with respect to the value reported in [25].

Optimal Perturbations in Channel Flow with No-Slip Walls
We first consider the plane Poiseuille flow with no-slip walls at Re = U c h/ν = 4000 (U c being the centerline velocity and h the half-height of the channel) in a domain l x × l y × l z = 2π × 2 × π, where y is the wall-normal direction, x and z are the streamwise and spanwise (periodic) directions. Nonlinear optimizations have been performed for several initial energies in the range [10 −7 , 10 −5 ], for small target times, T = 10, 20, 30, 40, 50 (for the largest ones, the value of the initial energy has been limited to 10 −6 in order to avoid transition at the target time, allowing convergence of the optimization algorithm). One can notice from Figure 5 that, depending on the initial energy and on the target time, different types of optimal perturbations have been found. For very low (O(10 −7 )) values of the initial energy (triangles in Figure 5), the perturbation is constituted by x-independent streamwise vortices at t = 0, as shown in the left frame of Figure 6, showing the optimal solution obtained for T = 20 and E 0 = 10 −7 . These lead to streamwise streaks at target time, similarly to the typical optimal perturbation found in a linearized framework in most of the shear flows [5]. However, when the value of the initial energy is slightly increased above the long-dashed line in Figure 5 the optimal perturbation becomes more fragmented in the streamwise and spanwise direction, characterized by localized vortices inclined with respect to the streamwise direction, as provided in the middle frame of Figure 6, showing the optimal solution obtained for T = 10, E 0 = 10 −6 .
It has been discussed in [26] that the upstream tilting with respect to the wall-normal direction, which is observed also in some linear optimal perturbation in non-parallel flows [11] is due to the Orr mechanism. The conservation of circulation allows a transient energy growth when the mean flow tilts downstream a spanwise-vorticity region initially inclined opposed to the base flow. On the other hand, in the linear case, no inclination of the vortices in the spanwise direction is observed, meaning that it is a direct effect of nonlinearity. It is also worth noticing that these "weakly" optimal perturbations are constituted by a basic flow structure which is replicated in the three flow directions. This building block is composed by spanwise-inclined streamwise vortices alternating in the spanwise direction and flanking a region of high negative streamwise disturbance. Notably, this is very similar to the basic building block obtained for the Blasius and the asymptotic suction boundary layer flows [20,32]. Slightly increasing the initial energy (for instance, for T = 10 and E 0 = 1.25 × 10 −6 ) does not affect the main structure of the perturbation, but leads to a localization of the perturbation on one of the two walls, as already observed in other shear flows [20,83]. Moreover, for T ≥ 20 and for a slight increase of the initial energy above the top solid line of Figure 5, a further localization is observed in the spanwise and streamwise directions, leading to an optimal perturbation characterized by only one building block of the previous solutions. These "highly nonlinear" optimal is shown in the right frame of Figure 6, and is very similar to that previously found in Couette flow [21][22][23]. Is it worth to remark that spatial localization allows an increase of the perturbation initial amplitude for a given initial energy. This enhances nonlinear effects leading to an energy growth much stronger than the linear one. As a result, these inclined vortical perturbations rapidly create bent, large-amplitude streaks that directly lead the flow to turbulence, bypassing the long phase of linear growth and saturation of small amplitude, streamwiseindependent streaks, characterizing (quasi-)linear optimal perturbations. Notice that these initial perturbations already contain the three main ingredients of the self-sustained process theorized in [84], namely, streaks, vortices and sinuous modulations of the streaks, rapidly allowing the establishment of a cycle sustaining turbulence. Figure 5. Optimizations performed depending on T and E 0 : triangles for linear optimals, circles for weakly nonlinear optimals, squares for highly nonlinear optimals, and diamonds for hairpin-like nonlinear optimals (the black symbols indicate optimal solutions which are shown in the following figures). The long-dashed line roughly represents the boundary between linear and weakly nonlinear solutions, whereas the solid line indicates the approximate energy threshold above which hairpin-like and highly nonlinear optimal perturbations can be found.Results are obtained for a plane Poiseuille flow with Re = 4000 and domain size l x × l y × l z = 2π × 2 × π. Figure 6. Different types of nonlinear optimal perturbations obtained depending on T and E 0 : the linear (black triangle in Figure 5), weakly nonlinear (black circles in Figure 5), and highly nonlinear optimal perturbations (black square in Figure 5) are represented using, respectively, for the linear optimal, iso-surfaces of the streamwise vorticity (red/green for positive/negative); for the weakly nonlinear optimal, iso-surfaces of the Q-criterion coloured in yellow/ green for positive/negative values of the streamwise vorticity; and for the highly nonlinear optimals, iso-surfaces of the streamwise vorticity (white/black for positive/negative) and of the streamwise velocity (yellow/blue for positive/negative values).
On the other hand, keeping the target time as small as T = 10 and further increasing the initial energy within the range 1.25 × 10 −6 < E 0 < 10 −5 , the nonlinear optimal perturbation changes again its main structure, localizing in the streamwise and spanwise direction but keeping the spanwise symmetry, as shown in the left frames of Figure 7. Once again, we observe tubes of counter-rotating vorticity placed side by side and tilted against the base flow, confirming the presence of linear energy growth mechanisms such as the Orr [85] and the lift-up [8] ones. However, the downstream tilting now causes the merging of one extremity of the two counter-rotating vortices, creating a characteristic hairpin vortex structure at target time (see the right frames of Figure 7). For a further increase of the initial energy, this hairpin-like basic structure observed in the top frames of the figure slightly extends in space, although remaining confined on only one of the walls and keeping a very similar flow structure. Notably, the optimal disturbance providing the highest energy growth is the most localized one, obtained for E 0 = 2 × 10 −6 . These results suggest that, for the considered flow, hairpin structures are highly energetic flow structures able to provide the maximum possible energy growth, much larger than that achievable under linearized conditions. In the next subsection of the paper, we would rapidly analyze the effect of superhydrophobic surfaces on different scenarios of transition, in order to select relevant values of the initial energy and target time for the nonlinear optimization in a turbulent channel flow over superhydrophobic walls.

Transition Scenarios in a Superhydrophobic Channel: A Review
To motivate our choice of target time and initial energy for the nonlinear optimizations in the SHS cases, we shall rapidly review how different transition scenarios in channel flow are influenced by the presence of slip at the wall. While the influence of SHS as a mean of passive flow control has been extensively studied in both laminar and turbulent regimes, to the authors knowledge only few studies have by now focused on the impact of SHS onto the laminar-turbulent transition process. The first work reporting about the potential of SHS in delaying the onset of transition [86] made use of an experimental pipe flow facility to show that the transition Reynolds number increased when superhydrophobic coating were applied to the walls of the pipe. Some years later, the authors of [76] found that introducing homogeneous streamwise slip stabilizes the Tollmien-Schlichting waves. Similar results were successively found in the case of anisotropic SHS [66,87]. Using direct numerical simulations with homogeneous slip as boundary condition, Min and Kim [76] investigated the effect of the imposed slip length on a specific transition path, the K-type scenario [88]. They showed that this particular type of transition can be delayed by wall slip; however, they did not explain how and why this effect is achieved, and they did not extend the analysis to other transition scenarios. The recent work of Picella et al. [64] has analyzed the influence of slip at the wall on three different transition scenarios in channel flow, that we rapidly review here: (i) the K-type transition scenario, which is triggered setting as initial disturbances a sum of a twodimensional TS wave with streamwise and spanwise wavenumber (α, β) = (1.12, 0.0), and a sum of equal and opposite oblique three-dimensional fundamental Tollmien-Schlichting waves characterized by spatial wavenumbers (α, β) = (1.12, ±2.10) with amplitude 0.033 and 0.0011, respectively; (ii) an initial linear optimal perturbation computed in the chosen flow conditions and scaled with a sufficiently high amplitude for triggering nonlinear effects and consequently leading to transition; and (iii) an ad hoc volume forcing, constructed by a superposition of optimal forcing functions, able to trigger a large-amplitude response in the flow as a consequence of receptivity, which is referred to as F-type transition [89]. In all these cases, slippery walls with slip length L s = 0.00, 0.01, 0.02 have been considered, with the highest value resulting from a SHS having a roughness texture period of L + ≈ 15 in the turbulent regime [73]. For detecting when (and whether) transition occurs, in Figure 8 we have tracked in time the evolution of the friction Reynolds number, defined as where · represents the spatial average computed on the wall-parallel planes x − z at a given time t and the laminar Reynolds number is set to the constant value Re = 5000 [88]. Figure 8, showing the friction Reynolds number evolution in the case of K-type transition, shows that in both laminar (t < 100) and turbulent (t > 500) regimes, the mean value of Re τ changes with L s . In fact, combining Equation (17) with Equation (19) one obtains Looking at the transitional phase of the curves in Figure 8, we observe that while Re τ increases monotonically for L s = [0.00, 0.01], this is not the case for L s = 0.02, where a transient growth of the friction Reynolds number is observed, followed by a slow decrease, corresponding to saturation, and a subsequent rapid increase towards the turbulent value. In the top right frame of Figure 8, a snapshot of the no-slip case is reported for comparison purposes. In this case, the two-dimensional TS wave first saturates, breaking its spanwise symmetry due to secondary instability [90]. This promotes the formation of λ-vortices at t ≈ 103, which then connect downstream forming the heads of characteristic hairpin vortices, shown in the snapshot of Figure 8 (top right) [91], which then lead the flow to turbulence [92]. In the slippery case with L s = 0.01, shown in the right middle frame of Figure 8, transition appears to be qualitatively similar to the no-slip case in its early stages (t < 100), whereas some notable differences arise in the late stages, where the λ− vortices weaken and stretch in the streamwise direction, creating an elongated and alternated streamwise velocity distribution at the wall. Then, at t ≈ 130 (see middle frame of Figure 8), hairpin vortices are created, but they show smaller heads with respect to the classical ones and straight long legs directly connected to the wall, which leave a clear footprint in the slip velocity at the wall. Breakdown to turbulence follows. This process is radically modified for the largest slip length considered here, L s = 0.02, provided in the bottom frame of Figure 8. Even if the initial condition is designed to promote the formation of λ vortices [91], the reduced wall shear induced by the presence of slip at the wall restrains their formation, inhibiting the consequent creation of hairpin vortices, leading to the transient decrease of the friction Reynolds number. However, at t ≈ 200, streaky structures are created as a response of the flow to the residual perturbations present in the computational domain [6]. For 180 < t < 200, these streaks increase their amplitude due to a transient growth mechanism triggered by flow receptivity. Once the streaks reach a threshold amplitude, nonlinear effects set in leading to a saturation of Re τ . Finally, at t > 400, secondary instability of these nonlinearly saturated flow occurs leading to breakdown to turbulence.
A radically different behavior is observed when the transition in the channel is initiated by optimal disturbances or by random noise [89]. In both cases, if the initial perturbation energy is set to a sufficiently high value (E(0) > 10 −6 ), after an initial linear growth phase (0 < t < 200), the streaks saturate nonlinearly. As soon as this amplitude overcomes a critical threshold [93], secondary instabilities of the nonlinearly saturated optimal streaks set in (t ≈ 400, see the left column of Figure 9 for the no-slip case), immediately followed by breakdown to turbulence. Introducing a slippery boundary does not modify this scenario, as depicted in the right column of Figure 9. Comparing these snapshots, one can see that the behavior for the cases L s = 0.00, 0.02 is virtually the same, except a time lag of about 20 time units. The same behavior is found for the F-type scenario, which is practically unaffected by the presence of slip at the wall. In fact, in both these cases the perturbations are localized in the bulk of the channel, as opposed to TS waves which develop close to the wall. As the superhydrobhobic surfaces act by modifying the flow shear near the wall, it appears clear why optimal streaks and centerline noise are both unaffected by this boundary modification even when nonlinearities kick in, as depicted in Figure 9, whereas structures that originate right from the wall, such as the λ vortices observed in the K-type scenario, are strongly influenced by the slip at the wall.

Nonlinear Optimal Perturbations in a Superhydrophobic Channel
The considerations made above explain the reason why, in order to find a relevant influence of the superhydrophobic walls on the nonlinear optimal perturbations, we choose to focus on high initial energies and short target times leading to hairpin-like optimal perturbations, which are created towards the rapid formation of λ vortices.
Thus, nonlinear optimizations of the flow in a channel with superhydrophobic walls have been carried out at target time T = 10, for four different values of the Reynolds number, and different initial energies. In particular, we have started at Re = 4000 with E 0 ≈ 2 × 10 −6 and then we have increased/decreased the initial energy until, quite suddenly, the initial optimal perturbation starts to localize in space leading to the previously shown hairpin-like solution, with a consequent increase of the energy gain of more than 10%. The value of E 0 for which this happens has been successively bisected for determining, within a two digits accuracy, the threshold energy E(0) th for which these hairpin-like optimal perturbations are found. The values of this energy threshold are shown in the left frame of Figure 10 for all of the considered values of Re and L s . Increasing the slip length, the energy threshold for obtaining the optimal hairpin-like solution slightly rises, confirming the stabilizing effect of a slippery wall on perturbation growth. However, for all the considered slip lengths the threshold initial energy scales with the Reynolds number with the approximated scaling law E(0) th ≈ Re −2.36 . This scaling law for the energy transition threshold, corresponding to about Re −1.2 for the threshold amplitude, is not far from the value provided by the experimental analysis of [94], who found a scaling law for the transition threshold amplitude of A ≈ Re −3/2 using wall-normal flow injection for inducing hairpin vortices. This may suggest that the mechanism of formation of hairpin vortices from λ− vortices remains robust even in the presence of slip, but a slightly higher threshold energy is needed for allowing the λ− vortices to grow. The corresponding energy gain at target time is provided in the right frame of Figure 10. For all of the considered Reynolds numbers and slip lengths, the energy increases of two orders of magnitude in only 10 time units. In all cases, the energy increases with Re but larger slip lengths are associated to lower gains. In particular, for Re = 5000 and L s = 0.05, the energy gain reached by the hairpin-like optimal is almost halved with respect to the no-slip case, indicating a strong effect of the boundary slip on the energy growth of these type of perturbations. As a reference, we recall that the energy growth of the linear optimal perturbation is affected by less than 2%, as reported in several references [63] and verified here using local nonmodal stability analysis. One can also observe that the energy gain of the hairpin-like optimal scales with the Reynolds number with a positive exponent which strongly depends from the imposed slip length, suggesting an even stronger effect of slip on nonlinear optimal growth at larger Reynolds numbers.

2000
3000 4000 5000 Re Despite the fact that the initial and final energies of these nonlinear optimal perturbations are noticeably affected by the presence of slip at the wall, the overall shape of the optimal perturbations remains robust with respect to this parameter. Figure 11 shows the isosurfaces of the streamwise perturbation for the optimal perturbations computed for slip lengths L s = 0, 0.02, Re = 4000 and for two initial energies close to the threshold E(0) th , namely E(0) = 1.5 × 10 −6 , 2.5 × 10 −6 . One can see that, for the lower energy a nonlocalized perturbation corresponding to the weakly nonlinear solution is found, whereas, increasing E 0 above the threshold, the localized hairpin-like solution is found. Most importantly, the spatial structure of the perturbation appears to be not affected at all by the presence of slip at the wall, indicating that the optimal mechanism leading to these two different types of solutions is indeed robust, despite the differences found in the corresponding energy gain.
Whereas, placing ourselves at the intermediate initial energy E(0) = 2 × 10 −6 leads to important differences in the perturbation shape, as a consequence of the slight change of the energy threshold E(0) th with the slip length. In fact, in the no-slip case this value of E(0) is slightly larger than E(0) th , thus we observe a hairpin-like solution, virtually identical to that obtained for E(0) = 2.5 × 10 −6 and shown in the bottom frames of Figure 11. While, in the slippery case, as this initial energy is slightly lower than E(0) th , a weakly nonlinear solution is obtained, virtually identical to these obtained for E(0) = 1.5 × 10 −6 and shown in the top frames of Figure 11. The time evolution of these two optimal perturbations is shown in Figure 12, providing the no-slip case on the left column, and the slippery case with L s = 0.02 on the right one. The hairpin-like optimal first generates one single pair of λ-vortices (Figure 12a,b, left) with a subsequent hairpin vortex that creates other ones (Figure 12c,d, left), leading to break down towards turbulence (Figure 12e, left). Whereas, the weakly nonlinear optimal perturbation creates from the first instants an array of staggered hairpin heads (Figure 12a, right), which are then connected by legs creating hairpin vortices (Figure 12b, right) that rapidly become turbulent (Figure 12c-e, right). Notice that this transition scenario in the presence of slip is very rapid, even faster than that observed in no-slip conditions [26], leading to transition in a few instants of time. Interestingly, transition is found to occur in both cases, although in the no-slip case a perturbation of much larger amplitude (equal energy, but much small extent) is needed to trigger it. Figure 11. Isosurfaces of the streamwise velocity component (yellow for positive and purple for negative values) of the nonlinear optimal perturbation for T = 10, Re = 4000, two initial energies indicated within the figure and for the no-slip (left column, (a,c)) and slippery case with L s = 0.02 (right column, (b,d)).

Discussion
The results provided above indicate that boundary slip may have a strong effect on the energy amplification of nonlinear coherent flow structures such as λ− and hairpin vortices, despite the fact that their overall shape barely changes due to the presence of the slip. These findings confirm the recent conclusions drawn by Davis and Park [68], who investigated the effect of slip surface on the transition dynamics of nonlinear flow structures. In particular, they found that the shape and position of critical-layer vortical structures similar to those typical of K-or H-type transition remain mostly unchanged no matter the slip lengths considered, despite their swirling amplitude may be reduced with slip length. On the other hand, in the case of nonlinear "sinuous" coherent structures placed in the core of the channel, transition is promoted by the presence of wall slip. In fact, slip appears to promote the propagation of these structures towards the wall inducing a consequent break up to turbulence through bypass-like transition. The authors conclude that slip surfaces promote the prevalence of strong wall-toward motions (sweeps events), that may destabilize core vortices, but stabilize those structures characterized by strong ejection events, allowing them to remain intact and propagate downstream for a longer time. Our results appear to confirm this scenario. In fact, the hairpin-like optimal perturbation studied in the present work has a critical-layer structure characterized by a strong ejection event (see in [25]). Although the presence of slip does not influence its shape, its growth in energy is considerably reduced by the presence of the slip, which is consistent with a stabilizing effect due to the a prevalent motion towards the wall. Whereas, transition triggered by the weakly nonlinear optimal with slip conditions appears to be even faster than that observed in no-slip conditions. This can be easily explained noticing that these vortical perturbations rapidly create sinuous streaks which bear many resemblances with the core structures mentioned by Davis and Park [68], whose transition is promoted by the presence of slip. This may appear in contrast with what found by Picella et al. [64], where slip appeared to have no effect at all on the transition behavior of core perturbations such as the linear optimal disturbances. However, this behavior could stem from the linear character of these last perturbations, which have been computed in the total absence of nonlinear effects.
Despite these results appear to suggest that the optimal energy growth mechanisms are indeed robust with respect to boundary slip, it should be considered that these investigations have been carried out using an homogenized approach for the SHS modeling. Thus, the effect of the gas-liquid interface motion on the optimal growth and on the consequent transition scenarios, remains to be investigated. Using a moving-boundary technique for modeling the SHS surfaces, Picella et al. [65] found that the strong ejection events induced by the interface motion can contrast the sweep-like effect of the boundary slip, leading to a transition scenario very similar to that observed in the no-slip case. This study suggests that the optimal solutions may remain robust also in the presence of interface motion. Nevertheless, this point remains for the moment open and it might be worth investigating in future studies.

Conclusions
Avoiding or delaying transition to turbulence is of a huge practical interest for applications requiring the transport of flows in channels or pipes, where turbulence is the primary cause of friction losses. As it is known, transition to turbulence induces a consistent drag increase in flows over bounding surfaces. Superhydrophobic surfaces (SHS) are currently used to reduce turbulence intensities in fully turbulent flows, but they have yet not been tested for delaying or avoiding turbulence. Hindering the laminar-turbulent transition process may potentially lead to a reduction of friction losses up to 50%, ensuring a huge energy saving in the flow transport process. Within this framework, this work aims at investigating how superhydrophobic surfaces affect optimal energy growth mechanisms relying on nonlinearity. Nonlinear optimizations based on a variational direct-adjoint procedure coupled with a gradient rotation algorithm have been carried out in a channel flow with no-slip or slippery boundaries, mimicking the presence of superhydrophobic surfaces. At first, optimizations have been carried out in the no-slip case, in order to identify relevant flow cases of potential interest. Changing the initial energy and the target time, four different families of nonlinear optimal perturbations have been found, characterized by inclined vortical structures with different localizations and symmetries. Among those, we have selected the hairpin-like optimal structure for the computations in the presence of boundary slip, due to recent results showing that hairpin vortices are strongly affected by wall slip. Increasing the slip length, the energy threshold for obtaining hairpin-like nonlinear optimal perturbations slightly rises, scaling approximately with Re −2.36 no matter the slip length. More importantly, the corresponding energy gain increases with Re with a slope that reduces with the slip length, being almost halved for the largest slip and Reynolds number (Re = 5000) considered. Moreover, the energy gain of the hairpin-like optimal scales with the Reynolds number with a positive exponent which strongly depends from the imposed slip length, suggesting an even stronger effect of slip on nonlinear optimal growth at Re > 5000. This suggests a strong effect of boundary slip on the energy growth of nonlinear perturbations, while small influence were found on nonmodal growth in linearized conditions. While the energy is considerably decreased, the shape of the optimal perturbation barely changes, indicating the robustness of optimal energy growth mechanisms with respect to wall slip. However, the amount of energy growth and consequent time to transition appears to be considerably influenced by wall slip, leading to damping effect in the case of hairpin-like structures, while an even more rapid transition is observed for core sinuous flow structures. These results appear to confirm the scenario suggested by [68], in which wall slip stabilizes hairpin-like coherent structures constituted by ejections, while destabilizes sinuous ones located in the core of the channel. Thus, the present work shows that wall slip has a strong potential in stabilizing optimal ejection-like events, although it can be not effective, or even detrimental, in the case of core perturbations. However, it has been recently shown [65] that, for sufficiently large (but wetting-stable) micro-roughnesses, the gas-liquid interface motion may partially hinder the transition delay effect of the wall slip. Thus, it would be crucial to investigate the robustness of these conclusions about the effect of SHS on optimal energy growth mechanisms with respect to the presence of gas-liquid interface motion. Until this computationally demanding task will be tackled, the validity of these results remains limited to the case of very small micro-roughnesses, for which the gas-liquid interface motion is negligible.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.