Identifying eigenmodes of averaged small-amplitude perturbations to turbulent channel flow

Eigenmodes of averaged small-amplitude perturbations to a turbulent channel flow – which is one of the most fundamental canonical flows – are identified for the first time via an extensive set of high-fidelity graphics processing unit accelerated direct numerical simulations. While the system governing averaged small-amplitude perturbations to turbulent channel flow remains unknown, the fact such eigenmodes can be identified constitutes direct evidence that it is linear. Moreover, while the eigenvalue associated with the slowest-decaying anti-symmetric eigenmode mode is found to be real, the eigenvalue associated with the slowest-decaying symmetric eigenmode mode is found to be complex. This indicates that the unknown linear system governing the evolution of averaged small-amplitude perturbations cannot be self-adjoint, even for the case of a uni-directional flow. In addition to elucidating aspects of the flow physics, the findings provide guidance for development of new unsteady Reynolds-averaged Navier–Stokes turbulence models, and constitute a new and accessible benchmark problem for assessing the performance of existing models, which are used widely throughout industry.

where u = (u, v, w) is the velocity, p is the pressure and Re is the Reynolds number. If Re is low enough such that flow is laminar, a simple unsteady solution of (1.1)-(1.2) can be obtained by imposing an initial condition u(x, y, z, 0) = (u IN (y), 0, 0), (1.4) that is independent of the wall-parallel coordinates x and z. In this case it is easy to verify that u = (u(y, t), 0, 0) at all times, the pressure gradient is independent of space and (1.1)-(1.2) reduces to the heat equation where d = dp dx . (1.6) To uniquely determine the solution, either the pressure gradient d or the mass flow rate must be prescribed. Assuming the latter, separation of variables leads to u(y, t) = U(y) + u (y, t), d(t) = D + d (t), (1.7a,b) where D is the constant pressure gradient associated with the steady state solution U(y), and u (y, t) = ∞ n=1 C n e λ n t u n (y), d (t) = (1.9a−c) All eigenvalues of this problem are negative, and hence u(y, t) → U(y) as t → ∞. Moreover the eigenvalues can be ordered such that λ n+1 > λ n for all n. The slowestdecaying eigenmode is given by u 1 (y) = sin(πy), d 1 = 0, λ 1 = −π 2 /Re, (1.10a−c) which is anti-symmetric. The second slowest-decaying eigenmode is given by u 2 (y) = cos(αy) − cos(α), d 2 = λ 2 cos α, λ 2 = −α 2 /Re, (1.11a−c) which is symmetric, where α = 4.4934 is the smallest solution of sin(α) = α cos(α). (1.12) Consequently, the long-time behaviour of the solution for a generic initial condition is typically u(y, t) ≈ U(y) + C 1 e λ 1 t u 1 (y) (assuming u 1 (y) is excited by the initial condition), whereas the long-time behaviour for symmetric initial conditions, which would not excite u 1 (y), is typically u(y, t) ≈ U(y) + C 2 e λ 2 t u 2 (y) (assuming u 2 (y) is excited by the initial condition).
1.3. Turbulent channel flow case If Re is high enough such that flow is turbulent, then the widely accepted triple decomposition (Hussain & Reynolds 1970) can be applied, and the velocity field can be written as where u = U + u is the averaged velocity, p is the averaged pressure and u u is called the Reynolds stress tensor. In the case of a plane channel flow the averaging can be done in either of the wall-parallel directions, or over an ensemble of channels. According to the widely accepted ergodic hypothesis these averaging methods and their combinations are equivalent, provided that the flow is statistically homogeneous in wall-parallel directions.
Analogous to the laminar case, one can attempt to obtain a simple solution by imposing an initial condition u(x, y, z, 0) = (u IN (y), 0, 0), (1.17) that is independent of the wall-parallel coordinates x and z. In this case it is easy to verify that u = (u(y, t), 0, 0) at all times, the averaged pressure gradient is independent of space and (1.14)-(1.15) reduces to Since (1.18) contains an unknown Reynolds stress term, it is not closed. We anticipate that the correct, yet unknown, closure exists. The objective of the current study is to ascertain whether in the limit of small-amplitude perturbations (1.18), with the correct closure, has solutions of a form analogous to (1.8), and in particular to identify the slowest-decaying anti-symmetric and symmetric eigenmodes that are analogous to (1.10) and (1.11). However, since the closure for (1.18) is unknown, the system cannot be solved directly. Therefore, the problem will be tackled via an extensive set of high-fidelity graphics processing unit (GPU) accelerated direct numerical simulations (DNS) of channel flow, from which averaged solutions will be obtained and analysed directly.
1.4. Related work The notion of a small-amplitude perturbation implies two flows are considered, one of which is designated a base flow, and the other as a sum of the base flow and the small-amplitude perturbation. Within this context, the base flow is often considered to be either an evolving time-dependent turbulent flow field, or an averaged flow field.
Analysis of small-amplitude perturbations to time-dependent turbulent base flows usually involves calculation of Lyapunov exponents, characterising the growth rate of the difference between the base flow and a small variation thereof. Recent studies in this area include e.g. Nikitin (2018). When flow is turbulent, small-amplitude perturbations defined in this way usually grow very quickly and then saturate at a finite amplitude, becoming in effect the difference between two unrelated realisations of the same turbulent state. Such small-amplitude perturbations are not considered in the present study.
Analysis of small-amplitude perturbations to averaged base flows began with the seminal work of Malkus (1956), who attempted to build a theory of turbulence based on an assumption that average velocity profiles are marginally stable with respect to infinitesimal perturbations as described by the linearised Navier-Stokes equations, whilst at the same time maximising energy dissipation rate. It was later shown that this assumption is not true (Reynolds & Tiederman 1967). However, the approach inspired further lines of research, including those motivated by the 762 A. S. Iyer, F. D. Witherden, S. I. Chernyshenko and P. E. Vincent idea that patterns, also called organised structures, observed in developed turbulent flows, could be explained by studying instabilities of the averaged base flow. Indeed, patterns in a finite-amplitude solution that evolved from a small-amplitude random perturbation are sometimes similar to patterns of the fastest growing infinitesimal disturbance. This is often the case for free-shear turbulent flows, which are strongly affected by the inviscid instability of the inflectional average flow profile (Gaster, Kit & Wygnanski 1985;Roshko 1993), but not for near-wall turbulent flows where viscosity becomes important, since it is often impossible to find the fastest growing infinitesimal disturbance, as all disturbances eventually decay if their evolution is calculated using the molecular viscosity. It was observed, however, that transiently growing solutions of the linearised Navier-Stokes equations exhibit patterns similar to those seen in developed turbulent flows. To investigate this further, the full Navier-Stokes equations can be rewritten as linearised equations forced by nonlinear terms, and then these nonlinear terms can be modelled as statistically stationary stochastic forcing (Farrell & Ioannou 1993), so that the amplitude of the fluctuations remains time independent. The patterns of the fluctuating motion are then a result of the filtering properties of the linearised Navier-Stokes operator i.e. its selectivity in amplifying or suppressing particular patterns. These filtering, or selectivity, properties can be studied directly using techniques known as input-output analysis, as well as indirectly by considering the pseudospectra of the operator (Trefethen et al. 1993) or optimal perturbations with maximum transient growth (Butler & Farrell 1993). The development of this philosophy proved fruitful, with results including elucidation of how streak spacing depends on wall distance (Chernyshenko & Baig 2005), how turbulent drag reduction due to spanwise wall oscillations depends on the frequency of the oscillations (Moarref & Jovanović 2012), how the angle of near-wall streaks depends on time in flows with drag reduction caused by wall oscillations (Blesbois et al. 2013) and how drag reduction by in-plane wavy motion of the wall depends qualitatively (although not quantitatively) on the frequency and streamwise wavenumber of this motion (Duque-Daza et al. 2012). Further progress, during which the above-described paradigm became known as a resolvent framework, was made when nonlinear effects in the form of triadic interactions were considered (see the review by McKeon (2017)). Another way of including nonlinear effects was proposed via a statistical state dynamics approach, in which they are restricted to the interaction of the streamwise average flow and perturbations (see the review by Farrell, Gayme & Ioannou (2017)). The above references serve as an entry point into this thriving area of research. However, we note that it is only indirectly related to the problem considered in the present study.
An alternative line of research stemming from the seminal work of Malkus (1956) has focused on the behaviour of averaged small-amplitude perturbations to averaged base flows. This is directly related to the problem considered in the present study. Following the observation that average turbulent flow profiles are stable with a laminar viscosity, Hussain & Reynolds (1970) introduced the now widely accepted triple decomposition, leading to the unclosed system (1.18). Importantly, it remains unclosed even after linearisation under the assumption of a small coherent component. Davis (1974) investigated deriving turbulent closures for averaged small-amplitude perturbations directly from the Navier-Stokes equations, and Gritsun & Branstator (2007) utilised the general fluctuation-dissipation theorem to gain additional insight. However, a full theoretical solution for the evolution of averaged small-amplitude perturbations is yet to be found. This leaves only three routes for further progress: the use of semi-empirical closures, experiments and DNS. Reynolds & Hussain (1972) themselves used a simple eddy-viscosity closure, and achieved good agreement with associated experimental data. However, more recent studies by Sen & Veeravalli (1998, Sen et al. (2007), which employed more complicated closures with anisotropic Reynolds stresses, produced very different results. In particular they observed growth as opposed to decay of averaged small-amplitude perturbations, although these results, too, were found to agree with their own experimental data. Such a contradiction highlights the challenges associated with making progress using semi-empirical closures and experiments. Finally, in terms of relevant DNS, Luchini, Quadrio & Zuccher (2006) studied the behaviour of a time-dependent channel flow subject to an impulse force, and more recently Russo & Luchini (2016) investigated a time-independent channel flow, where small perturbations were maintained by a constant body force. In the later study DNS results were found to disagree with those obtained using a turbulent closure, highlighting the fundamental difficulty of identifying the correct closure.

Overview
Our approach was based on performing ensembles of turbulent channel flow DNS in a finite channel. Specifically, ensembles of statistically independent snapshots of a fully developed turbulent channel flow were modified with either anti-symmetric or symmetric net-zero mass flux perturbations to their streamwise velocity, that depended solely on the wall-normal direction. The evolution of the ensemble-and wall-parallel averages of these perturbations were then analysed, and as the averaged perturbations became small their slowest-decaying anti-symmetric and symmetric eigenmodes were identified and extracted.
It was assumed that the averaged perturbations exhibit an initial nonlinear transient, followed by a linear decay phase as their amplitude becomes small. In this linear phase, it was further assumed that averaged perturbations can be represented as a sum of eigenmodes, each decaying at different rates, such that eventually after some time t s the slowest-decaying eigenmode dominates. It is clear that in order to identify this slowest-decaying eigenmode one must study the long-time behaviour of the averaged perturbations when t > t s . However, this poses a challenge. Specifically, the averaged perturbations will inherently include averaging errors, which fluctuate irregularly as a function of time, but remain finite. Hence after some time t e they are of the same order as the averaged perturbation. If the fluctuations are large enough such that t e < t s i.e. the averaging errors dominate the solution before only a single slowest-decaying eigenmode dominates, then it will be impossible to identify such a slowest-decaying eigenmode, even if it exists. The averaging errors can be reduced by averaging over more data. Hence, to overcome this challenge the ensemble of turbulent channel flows must be made large; such that t e t s . As will be shown, this requires a significant number of DNS to be undertaken, and indeed the approach has only recently become feasible with the emergence of high-fidelity GPU accelerated solvers.
The solver used to undertake the DNS here is described in § 2.2, the set-up for a single channel is described in § 2.3, the ensemble and averaging approach is described in § 2.4 and finally the procedure for identifying the slowest-decaying eigenmode is described in § 2.5.

Solver
All DNS were undertaken using PyFR ( reconstruction method of Huynh (2007), and can efficiently leverage the capabilities of modern massively parallel GPU platforms (Vincent et al. 2016). The compressible Navier-Stokes equations were solved at a Mach number of 0.1 to approximate incompressible conditions. A Rusanov Riemann solver was employed to calculate the inter-element fluxes, an explicit RK45[2R+] scheme (Kennedy, Carpenter & Lewis 2000) was used to advance the solution in time and all runs were performed using double precision arithmetic. Specifically, a patched version of commit ebf386bf4f30e9caaf1cbe8546169f6f200875a1 to the PyFR Git repository (e.g. https://github.com/vincentlab/PyFR) was employed. The patch, and the full version of the code, are provided as electronic supplementary material, available at https://doi.org/10.1017/jfm.2019.520.

Set-up
A schematic of the computational domain is shown in figure 1. It has a length of 8π in the x-direction, a height of 2 in the y-direction and a width of 4π in the z-direction. The origin is located at the centre of the domain. Periodic boundary conditions were applied in the xand z-directions, and an adiabatic no-slip condition was applied at y = ±1. Flow was driven through the domain via a volumetric momentum source, coupled with a feedback controller to ensure a constant bulk velocity of 1 in the x-direction. The viscosity was set to achieve a target friction Reynolds number of Re τ = 180, which corresponds to a bulk Reynolds number based on the bulk velocity and channel half-height of Re = 2767. Time units were based on the bulk velocity and channel half-height.
The domain was meshed with a structured grid of 62 × 19 × 60 hexahedral elements along the x-, yand z-directions respectively. The mesh was uniform in the xand z-directions, but non-uniform in the y direction, with resolution increased towards the channel walls. The flow solution within each element was represented using a fourth-order polynomial. Mesh independence, and comparisons with published direct numerical simulation data, are shown in appendix A and appendix B, respectively. The mesh is provided as electronic supplementary material.

Ensembles and averages
Consider an ensemble Π of statistically independent snapshots of a fully developed channel flow solution, obtained in the domain defined in § 2.3. The x-velocities u(x, y, z) e from each instance e of Π can be averaged as follows to obtain (2.1) which in the limit N → ∞ recovers the exact averaged x-velocity of the fully developed base flow. In the present study U(y) was calculated using N = 768 snapshots. These were obtained by running a series of independent statistically stationary channel flow simulations, from which solutions were extracted and stored every 4 time units. Now consider that each snapshot in Π is perturbed with either an anti-symmetric functionū to create an ensemble Π a of initial conditions with a purely anti-symmetric initial perturbation, or a symmetric function to create an ensemble Π s of initial conditions with a purely symmetric initial perturbation, where we note that bothū aIN (y) andū sIN (y) have net-zero mass flux, depend solely on the wall-normal direction, are zero at the walls, but are otherwise arbitrary. Finally consider restarting new simulations at t = 0 from the initial conditions in each ensemble Π a/s . The resulting x-velocities u a/s (x, y, z, t) e obtained from each instance e of Π a/s can be recorded, and their perturbations from the fully developed base flow U(y) can be averaged as follows to obtain (2.4) which in the limit N → ∞ recovers the exact averaged x-velocity perturbation after a purely anti-symmetric/symmetric initial perturbation. In the present studyū a/s (y, t) were each calculated using simulation results from N = 768 initial conditions. It will be demonstrated a posteriori that use of N = 768 snapshots to calculate U(y), and use of simulation results from N = 768 initial conditions to calculateū a/s (y, t), is sufficient to suppress averaging errors to the extent that single and hence slowestdecaying anti-symmetric/symmetric eigenmodes can be identified.

Identification of the slowest-decaying eigenmodes
We sought to identify the eigenvalue λ a and eigenfunction u a (y) associated with the slowest-decaying anti-symmetric eigenmode by studyingū a (y, t), since it was generated from a purely anti-symmetric initial perturbation. Similarly, we sought to identify the eigenvalue λ s and eigenfunction u s (y) associated with the slowest-decaying symmetric eigenmode by studyingū s (y, t), since it was generated from a purely 766 A. S. Iyer, F. D. Witherden, S. I. Chernyshenko and P. E. Vincent symmetric initial perturbation. To begin, the anti-symmetric component ofū a (y, t), denotedū aa (y, t), was extracted. This filtered out any symmetric components possibly generated during the nonlinear transient, and thus facilitated identification of the slowest-decaying anti-symmetric eigenmode. Similarly the symmetric component ofū s (y, t), denotedū ss (y, t), was extracted. This filtered out any anti-symmetric components possibly generated during the nonlinear transient, and thus facilitated identification of the slowest-decaying symmetric eigenmode. Subsequently, a single anti-symmetric/symmetric eigenmode was fitted to a series of time-windowed sections ofū aa/ss (y, t), with a view to identifying if sections ofū aa/ss (y, t) existed after the nonlinear transient phase, and after all but one of the eigenmodes had decayed, but before averaging errors became significant, where it could be well approximated by a single, and hence slowest-decaying anti-symmetric/symmetric eigenmode. Specifically, fitting was attempted in a series of windows t ∈ [t * , t * + t * ], for an equispaced distribution of window start times t * in the range 0 to 120, with candidate eigenvalues λ * a/s (t * ) and 'shifted' candidate eigenfunctionsû * a/s (y, t * ) associated with a single slowest-decaying anti-symmetric/symmetric eigenmode identified as solutions to the optimisation problem and associated relative errors R * a/s (t * ) obtained as where the 'shifted' candidate eigenfunctions were represented by 95 degrees of freedom, using the same piecewise polynomial discretisation as for the DNS simulations, and minimisation was performed using the Levenberg-Marquardt algorithm via the LMFIT package (Newville et al. 2014). Note that for technical reasons the time origin of the eigenmode in (2.6) was shifted by t * prior to fitting, hence whyû * a/s (y, t * ) are referred to as 'shifted' candidate eigenfunctions. Once obtained, these can be 'unshifted' to obtain candidate eigenfunctions u * a/s (y, t * ) = e −λ a/s t * û * a/s (y, t * ), (2.8) that can be directly compared. If a single slowest-decaying anti-symmetric/symmetric eigenmode is identifiable then there should exist a range of t * where: (i) The candidate eigenvalues λ * a/s (t * ) are approximately independent of t * , with a value that can be ascribed to the true eigenvalue λ a/s . (ii) The relative errors R * a/s (t * ) are small. (iii) The candidate eigenfunctions u * a/s (y, t * ) are approximately independent of t * , with a form that can be ascribed to the true eigenfunction u a/s (y). We note that, strictly, the above approach will identify slowest-decaying antisymmetric/symmetric eigenmodes excited by the initial perturbations (2.2) and (2.3), respectively. It is highly likely that these are also indeed the slowest-decaying anti-symmetric/symmetric eigenmodes per se. However, to confirm, it would be necessary to repeat the procedure using data from additional ensembles with various initial perturbations, and check that the eigenmodes all agree. Given the cost associated with DNS, such a check was not feasible here.
Finally, we comment that by applying a control volume approach in the normal way to the full Navier-Stokes (1.1) and continuity equations (1.2), averaging over the xand z-directions, and using the fact that the bulk velocity is constant, one can obtain an expression for the averaged pressure gradient perturbation associated with each eigenmode dp a/s dx (t) = Re e λ a/s t 2Re du a/s (y) dy +1 −1 .
(2.9) 3. Results which was obtained by nonlinear least squares fitting of a rational polynomial to U(y). . Plots of candidate eigenvalues λ * a (t * ) and relative fitting errors R * a (t * ) as a function of window start time t * , obtained using the procedure outlined in § 2.5 with a window size of t * = 30.

Slowest-decaying anti-symmetric eigenmode
eigenmode is real valued. Therefore, we assume that the corresponding eigenvalue λ a and eigenfunction u a (y) are real valued before proceeding with the analysis. Figure 4 plots candidate eigenvalues λ * a (t * ) and relative fitting errors R * a (t * ) as a function of window start time t * , obtained using the procedure outlined in § 2.5 with a window size of t * = 30. Within 30 t * 75 the candidate eigenvalues λ * a (t * ) are approximately constant and the relative fitting errors R * a (t * ) are small. the average of four normalised candidate eigenfunctions within 30 t * 75, obtained using the procedure outlined in § 2.5, and a window size of t * = 30. Figure 5 also plots their associated standard deviation, which is seen to be small.
We therefore conclude that within the interval 30 t * 75 the anti-symmetric component of the averaged x-velocity perturbation after a purely anti-symmetric initial perturbation can be described by a single, and hence slowest-decaying, anti-symmetric eigenmode. The associated eigenvalue can be approximated as . Plots ofũ * a (y), the average of four normalised candidate eigenfunctions within 30 t * 75, obtained using the procedure outlined in § 2.5, and a window size of t * = 30 (black line), along with their standard deviation (grey shaded region), and u a (y), obtained by least squares fitting of a rational polynomial toũ * a (y) (dashed white line). Note that u * a (y) (black line) and u a (y) (dashed white line) are almost exactly superposed.
obtained from the average of λ * a (t * ) within the interval 30 t * 75, which has a standard deviation of 0.001. The associated eigenfunction u a (y) can be approximated as u a (y) = (1 − y)(1 + y) −5.52y 7 + 11.81y 5 − 9.77y 3 + 3.52y 0.50y 4 − 1.47y 2 + 1 , which is also plotted in figure 5, and was obtained by nonlinear least squares fitting of a rational polynomial toũ * a (y). As a check, figure 6 comparesū aa (y, t), the anti-symmetric component of the averaged x-velocity perturbation after a purely anti-symmetric initial perturbation with the eigenmode C a u a (y)e λ a t within 30 t 75, where C a is real and depends on the initial condition, and is chosen here via a nonlinear least squares fit. The contour plots are similar, providing further confidence that a single, and hence slowest-decaying, anti-symmetric eigenmode has been identified.
We note that the shape of the eigenfunction is similar to that of a sine curve, which is the corresponding anti-symmetric laminar solution defined by (1.10). However, the eigenvalue λ a is found to be 11.8 times bigger than the analogous laminar value obtained from (1.10) with Re = 2767. This is in line with expectations since eddy viscosity is typically greater than molecular viscosity. Figure 7 plotsū as (0, t), the symmetric component of the averaged x-velocity perturbation, after a purely anti-symmetric initial perturbation, at y = 0. We note that it is non-zero, with a maximum value that is of the same order as the square of the amplitude ofū aa (0.76, t) shown in figure 3, i.e. the initial purely anti-symmetric perturbation generates a symmetric component. This is in contrast to the laminar case, where an anti-symmetric perturbation will always remain purely anti-symmetric, since the governing equation for a unidirectional laminar flow is linear. It can be shown (see appendix C) that if the evolution of the averaged perturbation is governed by a symmetric linear operator, and the initial perturbation is anti-symmetric, then the averaged perturbation must remain purely anti-symmetric for all time. Our observation of the counter indicates that the full governing equations for the averaged perturbation FIGURE 6. Contours ofū aa (y, t), the anti-symmetric component of the averaged x-velocity perturbation, after a purely anti-symmetric initial perturbation (a) and the eigenmode C a u a (y)e λat (b). must be nonlinear, even for a uni-directional channel flow case. Consequently, any turbulence models should remain nonlinear even for a uni-directional flow, a corollary of which is the requirement that for models employing eddy viscosity, such viscosity cannot be independent of the flow. The well-known Prandtl model, in which the eddy viscosity where l(y) is the mixing length, satisfies this requirement. However, the assumption that ν t = ν t (y), often used in studies of small averaged perturbations, does not. This finding is in line with results from previous studies into the response of turbulent flow to constant forcing (Russo & Luchini 2016;Luchini 2018). Finally, we note that in accordance with (2.9) the averaged pressure gradient perturbation associated with the anti-symmetric eigenmode is zero, since the eigenfunction u a (y) is anti-symmetric.
3.3. Slowest-decaying symmetric eigenmode Figure 8 plotsū ss (0, t), the symmetric component of the averaged x-velocity perturbation, after a purely symmetric initial perturbation, at y = 0. The evolution is oscillatory, which suggests that the slowest-decaying symmetric eigenmode is complex valued. Therefore, we assume that the corresponding eigenvalue λ s and eigenfunction u s (y) are complex valued before proceeding with the analysis. Figure 9 plots the real and imaginary parts of candidate eigenvalues λ * s (t * ) and relative fitting errors R * s (t * ) as a function of window start time t * , obtained using the procedure outlined in § 2.5 with a window size of t * = 30. Within 30 t * 75 the real and imaginary parts of candidate eigenvalues λ * s (t * ) are approximately constant and the relative fitting errors R * a (t * ) are small. Figure 10 plots the real and imaginary parts ofũ * the average of four normalised candidate eigenfunctions within 30 t * 75, obtained using the procedure outlined in § 2.5, and a window size of t * = 30. Figure 10 also plots their associated standard deviations, which is seen to be small.
We therefore conclude that within the interval 30 t * 75 the symmetric component of the averaged x-velocity perturbation after a purely symmetric initial perturbation can be described by a single, and hence slowest-decaying, symmetric eigenmode. The associated eigenvalue can be approximated as obtained from the average of λ * s (t * ) within the interval 30 t * 75, which has a standard deviation of 0.003 + i0.004. The real and imaginary parts of the associated . Plots the real and imaginary parts of candidate eigenvalues λ * s (t * ) and relative fitting errors R * s (t * ) as a function of window start time t * , obtained using the procedure outlined in § 2.5 with a window size of t * = 30. . Plots of the real and imaginary parts ofũ * s (y), the average of four normalised candidate eigenfunctions within 30 t * 75, obtained using the procedure outlined in § 2.5, and a window size of t * = 30 (black lines), along with their standard deviations (grey shaded regions), and the real and imaginary parts of u s (y), obtained by least squares fitting of rational polynomials to the real and imaginary parts ofũ * s (y), respectively (dashed white lines). Note that the real and imaginary parts ofũ * s (y) (black lines) and the real and imaginary parts of u s (y) (dashed white lines) are almost exactly superposed. eigenfunction u s (y) can be approximated as −0.19y 6 + 6.20y 4 − 6.56y 2 + 0.96 2.75y 6 − 4.23y 4 + 0.60y 2 + 1 , (3.8) Im{u s (y)} = (1 − y)(1 + y) 3.66y 8 − 8.72y 6 + 7.44y 4 − 1.41y 2 − 0.25 1.64y 6 − 2.11y 4 − 0.48y 2 + 1 , (3.9) which are also plotted in figure 10, and were obtained by nonlinear least squares fitting of rational polynomials to Re{ũ * s (y)} and Im{ũ * s (y)}, respectively. We note that since the problem is real valued, the complex conjugate of λ s and u s (y) will also form a valid eigenmode.
As a check, figure 11 comparesū ss (y, t), the symmetric component of the averaged x-velocity perturbation after a purely symmetric initial perturbation with the real part of the eigenmode Re{C s u s (y)e λ s t } within 30 t 75, where C s is complex and depends on the initial condition, and is chosen here via a nonlinear least squares fit. The contour plots are similar, providing further confidence that a single, and hence slowestdecaying, symmetric eigenmode has been identified.
We note that in contrast to the corresponding symmetric laminar solution defined by (1.11), the eigenvalue and eigenfunction are complex. This implies that the corresponding linear operator cannot be self-adjoint, even in the case of a uni-directional flow. Moreover, unlike for the anti-symmetric perturbation, the shape of the eigenfunctions differ substantially between the laminar and turbulent cases. In particular, the symmetric eigenfunction has a large magnitude approximately 12 wall units from the wall, in a region where turbulent fluctuations are known to attain their maximum value. We can hypothesise that the complex eigenvalue is a result of an interplay between two physical processes, one occurring in the bulk of the flow, and another associated with turbulent fluctuations concentrated near the walls. Also, the real part of the eigenvalue Re{λ s } is found to be 6.72 times bigger than the analogous laminar value obtained from (1.11) with Re = 2767. This is in line with expectations since eddy viscosity is typically greater than molecular viscosity. Figure 12 plotsū sa (0.76, t), the anti-symmetric component of the averaged x-velocity perturbation, after a purely symmetric initial perturbation, at y = 0.76. We note that it is zero i.e. the initial purely symmetric perturbation does not generate an anti-symmetric component. This is expected since the mean flow, the initial perturbation, and the boundary conditions are all symmetric.

Conclusion
Eigenmodes of averaged small-amplitude perturbations to a turbulent channel flow -which is one of the most fundamental canonical flows -have been identified for the first time via an extensive set of high-fidelity GPU accelerated DNS. Specifically, for a turbulent channel flow with Re τ = 180 and a constant bulk velocity, the slowest-decaying anti-symmetric eigenmode was found to have an eigenvalue λ a = −0.042 and an eigenfunction u a (y) approximated by the analytical expression (3.4), and the slowest-decaying symmetric eigenmode was found to have an eigenvalue λ s = −0.049 − i0.069 and an eigenfunction u s (y) approximated by the analytical expressions (3.8) and (3.9). While the system governing averaged small-amplitude perturbations to turbulent channel flow remains unknown, the fact such eigenmodes could be extracted constitutes direct evidence that it is linear. Moreover, given λ s and u s (y) were found to be complex this unknown linear system cannot be self-adjoint, even for the case of a uni-directional flow. In addition to elucidating aspects of the flow physics, the findings provide guidance for development of new unsteady Reynolds-averaged Navier-Stokes turbulence models, and constitute a new and accessible benchmark problem for assessing the performance of existing models, which are used widely throughout industry.
Future studies can employ the approach adopted here to identify further eigenmodes of averaged small-amplitude perturbations to turbulent channel flow. These could Configuration n x × n y × n z p include e.g. those associated with higher Reynolds numbers, and those associated with perturbations that vary in wall-parallel directions; with the objective of obtaining a more extensive set of eigenmodes analogous to those given by the Orr-Sommerfeld equation for laminar channel flow. has been achieved. Figure 14 shows variation with p for fixed n x , n z and n z . The two highest polynomial orders (62-19-60-p4, 62-19-60-p5) also produce visually identical results, also indicating that mesh independence has been achieved. One can conclude that the 62-19-60-p4 configuration is sufficient to achieve mesh independence of the velocity.
A.3. Velocity fluctuations Figures 15 and 16 show plots of x-velocity fluctuations u rms , y-velocity fluctuations v rms and z-velocity fluctuations w rms , as a function of y. Figure 15 shows variation with n x , n y and n z for fixed p. The three finest meshes (52-16-50-p4, 62-19-60-p4, 82-25-80-p4) produce visually identical results, indicating mesh independence has been achieved. Figure 16 shows variation with p for fixed n x , n z and n z . The two highest polynomial orders (62-19-60-p4, 62-19-60-p5) also produce visually identical results, also indicating that mesh independence has been achieved. One can conclude that the 62-19-60-p4 configuration is sufficient to achieve mesh independence of the velocity fluctuations.

Appendix B. Comparison with published DNS data
Results obtained using the 62-19-60-p4 configuration were compared with DNS results at the same Re τ obtained by Kim, Moin & Moser (1987). Figure 17 compares plots of u + as a function of y + , and u rms , v rms and w rms as a function of y. There is good agreement in the inner, transition and outer layers. This provides further evidence that the 62-19-60-p4 configuration effectively resolves the flow physics.

(C 6)
Therefore, if E t is a symmetric linear operator, and the initial perturbation is antisymmetric, then the ensemble-and x-z-averaged perturbation must remain purely antisymmetric for all time.