Optimal energy growth in pulsatile channel and pipe flows

Abstract Pulsatile channel and pipe flows constitute a fundamental flow configuration with significant bearing on many applications in the engineering and medical sciences. Rotating machinery, hydraulic pumps or cardiovascular systems are dominated by time-periodic flows, and their stability characteristics play an important role in their efficient and proper operation. While previous work has mainly concentrated on the modal, harmonic response to an oscillatory or pulsatile base flow, this study employs a direct–adjoint optimisation technique to assess short-term instabilities, identify transient energy-amplification mechanisms and determine their prevalence within a wide parameter space. At low pulsation amplitudes, the transient dynamics is found to be similar to that resulting from the equivalent steady parabolic flow profile, and the oscillating flow component appears to have only a weak effect. After a critical pulsation amplitude is surpassed, linear transient growth is shown to increase exponentially with the pulsation amplitude and to occur mainly during the slow part of the pulsation cycle. In this latter regime, a detailed analysis of the energy transfer mechanisms demonstrates that the huge linear transient growth factors are the result of an optimal combination of Orr mechanism and intracyclic normal-mode growth during half a pulsation cycle. Two-dimensional sinuous perturbations are favoured in channel flow, while pipe flow is dominated by helical perturbations. An extensive parameter study is presented that tracks these flow features across variations in the pulsation amplitude, Reynolds and Womersley numbers, perturbation wavenumbers and imposed time horizon.


Introduction
Pulsatile flows are a common phenomenon in a variety of engineering flows, and they are ubiquitous in physiological configurations. The pulsatile flow through tubular geometries plays a key role in the haemodynamic system of many species as it is responsible for the transport of oxygenated blood to the organs and muscular tissue (Ku 1997;Pedley 2000). While in many of these configurations inertial effects are too weak to cause and sustain turbulent fluid motion, a variety of cardiovascular diseases can be linked to flow instabilities in the arteries (Chiu & Chien 2011). In addition, geometric modifications of the standard fluid-carrying vessels, such as stenoses, aneurysms or other pathologies, further amplify adverse flow effects and aggravate physiological consequences. For these reasons, a better understanding of pulsatile flows, and the perturbation dynamics they support, would be beneficial, if not mandatory, for improved diagnostics as well as the design of advanced medical devices.
Despite their importance in medical and engineering applications, pulsatile flowsand in particular their stability characteristics -have received far less attention than their steady analogues. Pulsatile flows comprise a steady as well as a time-periodic component. This is in contrast to oscillatory flows which consist of a harmonic part, but lack a steady background flow. The periodic time dependence precludes a standard modal approach, based on temporal Fourier normal modes, and instead calls for more complex methods, such as Floquet analysis. Furthermore, pulsatile flows are governed by a far larger suite of parameters than steady flows: besides the common Reynolds number Re and the wavenumbers of the perturbations, pulsatile flows depend on the pulsation amplitudes and the non-dimensional frequency (the Womersley number Wo). For a non-modal analysis, the time horizon over which growth or decay is measured and the phase shift of the perturbation within a base-flow cycle have to be accounted for as well. Within this high-dimensional parameter space, a rich and varied perturbation dynamics can be observed, with important transitions between distinct flow behaviours.
The stability of pulsatile flow has been addressed by a few key studies that laid the foundation for our current understanding of its perturbation dynamics. An account of the pertinent body of literature has been presented in Pier & Schmid (2017) with emphasis on the modal treatment via Floquet analysis. A resume of earlier work on general time-periodic flows has been presented in Davis (1976). Further notable work by von Kerczek (1982) has built on this foundation and established a framework for the analysis of flows with a harmonic base flow. Generic configurations such as a Stokes layer (Blennerhassett & Bassom 2002) or channel and pipe flow with time-periodic pressure gradients (Thomas et al. 2011), have been investigated with modal techniques and have been mapped out as to their stability characteristics across a range of governing parameters. The influence of wall modifications, such as stenoses or aneurysms, on the overall stability behaviour has been addressed via numerical simulations (see, e.g. Blackburn, Sherwin & Barkley 2008;Gopalakrishnan, Pier & Biesheuvel 2014).
The role of pulsation in the transition from laminar to turbulent pipe flow has been recently investigated by Xu et al. (2017) and Xu & Avila (2018). These studies in particular concentrated on the emergence and life cycle of localised 'puffs', together with their role in triggering transition in the presence of a pulsating flow component, since the occurrence of turbulent bursts in each cycle has been found to be sensitive on flow parameters and configuration details. A strong influence of the Womersley number has been reported, and a distinct regime-switching across three proposed parameter regions has been observed (Xu et al. 2017). These experimental findings have been further corroborated by direct numerical simulations initiated by a localised perturbation (Xu & Avila 2018). The earlier numerical study by Tuzi & Blondeaux (2008) concluded that at moderate but subcritical Reynolds numbers only parts of the harmonic cycle (around flow reversal) support turbulent flow via an instability and an associated break of the flow's symmetry.
While the early body of literature on time-periodic flows has concentrated on a modal (Floquet) approach, more recent studies have employed an initial-value perspective on the analysis of perturbation dynamics and energy growth. Biau (2016) has analysed the generic oscillatory Stokes layer as to its potential to support transiently growing perturbations over a forcing cycle. This study isolated the Orr mechanism as the dominant process by which energy amplification could be achieved efficiently for sufficiently high unsteady amplitudes. In particular the decelerating part of the forcing cycle has been identified as prone to strong non-modal growth. Complementary nonlinear simulations further verified that triggering by these mechanisms can yield subcritical transition to turbulence. A similar technique has been applied in a recent study by Xu, Song & Avila (2021) for oscillatory and pulsating pipe flow. Among others, they have reported that pulsating pipe flows are generally dominated by helical perturbations. In accordance with Biau (2016), a strong Orr-type mechanism has been found to dominate, once a threshold pulsation amplitude has been exceeded. Again, only half of the forcing cycle supported growth of the kinetic perturbation energy; disturbances have been observed to rapidly reach energy levels that facilitate a transition to turbulent fluid motion, often via localised disturbances.
These latter studies advocate the treatment of pulsatile flow as a generally time-dependent flow, distinct from a periodic Floquet ansatz. Over the past decades, the application of these non-modal techniques to hydrodynamic stability calculations has resulted in a more complete understanding of shear-driven instability phenomena. The generality of this approach (Schmid 2007) is well-suited for assessing pulsatile flow over a range of time scales, thus mapping out the optimal perturbation dynamics over partial and multiple pulsation cycles. This non-modal approach for time-dependent flows is based on a variational principle arising from a partial-differential-equation-constrained optimisation problem. It results in a direct-adjoint system of equations (Luchini & Bottaro 2014) that produce the maximum energy growth of perturbations over a prescribed time horizon. Time-dependent base flows are treated naturally within this formalism, and short-term energy amplification mechanisms, for example over a partial pulsation cycle, can be detected and extracted effectively. Over the past years, this computational framework has been successfully brought to bear on a variety of complex flow configurations (see, e.g. Magri (2019) and Qadri et al. (2021) for applications in reactive flows), and has furnished quantitative stability measures beyond the time-asymptotic limit and without the need for simplifying assumptions.
This article follows up on and extends earlier work (Pier & Schmid 2017) that demonstrated the influence of a pulsating flow component on the stability of channel flow via a linear (Floquet) and nonlinear analysis. In this present study, we focus on non-modal effects and the occurrence of transient energy amplification mechanisms under conditions that are asymptotically stable, both for rectangular channel and cylindrical pipe flows. The unsteady nature of the base flow lends itself to a formulation as a partial-differential-equation-constrained optimisation problem for the maximum energy gain which is subsequently solved by a variational approach based on direct-adjoint looping.
The main finding, and significance, of our investigation consists of the quantification of extremely large transient growth, brought on by the unsteady nature of the base flow. By considering both channel and pipe flows and carefully studying energy transfer mechanisms, we identify the fundamental mechanisms responsible for this huge growth, common to both geometries. This amplification potential translates directly into a strong sensitivity for the rise of coherent structures over one or many pulsatile cycles. While this feature of pulsatile flows has been observed and reported in previous studies, an encompassing treatment of this phenomenon, including its presence in parameter space and its manifestation in dominant spatial structures, is still missing in the literature on unsteady flows. Our findings also have a direct connection to classifying transition scenarios in wall-bounded flows under the influence of cyclic base flow variations, thus extending the classical scenarios for steady flows and potential routes for the transition to turbulence occurring during part of the cycle.
Despite our attempt to analyse pulsating channel and pipe flows comprehensively, judicious choices had to be made to arrive at an emerging picture for the perturbation dynamics prevailing in these configurations. The ensuing parameter ranges have been selected to capture the most compelling and representative flow phenomena, while limiting our focus to flows encountered in physiological and medical situations. Haemodynamic applications, across a range of blood vessel geometries, are well covered by our choice of parameters. Nonetheless, configurations outside this parameter range are touched upon as well, to establish continuity or bifurcations in flow behaviour and to connect to other studies that investigate such parameter regimes in more detail, e.g. Xu et al. (2021).
The present paper represents the culmination of several years of work; a preliminary version of the main results has been presented at the 12th European Fluid Mechanics Conference in Vienna (Pier & Schmid 2018).

Flow configurations and governing equations
This investigation considers viscous incompressible flow through infinite channels and pipes of constant diameter. In this context, a flow is characterised by a velocity vector field u(x, t) and a scalar pressure field p(x, t) that depend on position x and time t and are governed by the Navier-Stokes equations ∂u ∂t where ν is the kinematic viscosity of the fluid, and the pressure has been redefined to eliminate the constant fluid density. The channel-flow configuration calls for a formulation using Cartesian coordinates, while cylindrical coordinates are appropriate for pipe flows. In order to address both configurations with similar mathematical and numerical tools, we adopt a general formalism using three spatial coordinates x 0 , x 1 , x 2 and associated velocity components u 0 , u 1 , u 2 . When analysing channel flow with respect to a Cartesian reference frame, the variables x 0 , x 1 and x 2 denote wall-normal, streamwise and spanwise coordinates, respectively, while they stand for radial, streamwise and azimuthal coordinates when studying pipe flow in a cylindrical setting. Whatever the configuration, the flow domain corresponds to |x 0 | < D/2 where D is the channel or the pipe diameter, and no-slip boundary conditions prevail along the solid walls at |x 0 | = D/2.
A formulation of the incompressible Navier-Stokes equations ((2.1) and (2.2)) in cylindrical coordinates comprises more terms than one in Cartesian coordinates. Nevertheless, the resulting equations have a very similar structure, and the above notations allow us to cast the governing equations into a single general system of partial differential equations, pertaining to both channel and pipe configurations, the details of which are given in Appendix A.

Base flows and non-dimensional control parameters
Pulsatile base flows driven by a spatially uniform and temporally periodic streamwise pressure gradient are obtained as exact solutions of the Navier-Stokes equations and consist of a velocity field in the streamwise x 1 -direction with profiles that only depend on time t and on the wall-normal/radial coordinate x 0 . Denoting by Ω the pulsation frequency, the base velocity profiles may be expanded as temporal Fourier series, and are associated with a periodic flow rate In the above expressions, the conditions ] ensure that all flow quantities are real (with denoting a complex conjugate).
By invariance of these base flows in the streamwise x 1 -direction, the different harmonics in the expansion (3.1) are not coupled through the nonlinear terms of the Navier-Stokes equations and the velocity components U (n) 1 (x 0 ) are analytically obtained by solving simple differential equations derived for each harmonic component.
following Womersley (1955), the profiles U (n) 1 (x 0 ) are obtained in terms of Bessel functions in cylindrical coordinates corresponding to pipe flows, while they are obtained in terms of exponential functions in Cartesian coordinates corresponding to channel flows.
Pulsatile channel or pipe flows are characterised by the Womersley number which is a non-dimensional measure of the pulsation frequency, and may be interpreted as the ratio of the pipe radius (or the channel half-diameter) to the thickness δ = √ ν/Ω of the oscillating boundary layers developing near the walls. A pulsatile base flow is then completely specified by the Fourier components Q (n) of its flow rate (3.2), and the velocity profiles of the different harmonics (3.1) are obtained as In the above expression, A denotes the relevant measure of the cross-section (A = D for channels and A = πD 2 /4 for pipes) and the function W is the normalised velocity profile pertaining to each harmonic component. The analytic expressions of W for channel and pipe flows are given in Appendix B.
In this investigation, we only consider pulsatile flows with a non-vanishing mean flow rate Q (0) . Thus, the definition of the Reynolds number may be based on mean velocity Q (0) /A, diameter D and viscosity ν, leading to Re ≡ Q (0) ν for channels and Re ≡ Q (0) ν 4 πD for pipes. (3.5) Moreover, using Q (0) as reference, the flow rate waveform is completely determined by the non-dimensional ratiosQ corresponding to the amplitude (and phase) of the oscillating flow rate components (n > 0) relative to the mean flow.
In order to reduce the dimensionality of the control-parameter space for the rest of this paper, we will only consider base flow rates with a single oscillating component where the pulsation amplitudeQ ≡ 2Q (1) may be assumed real without loss of generality. Note that the theoretical and numerical methods developed for the present investigation are also suitable for studying the dynamics of pulsating base flows with higher harmonic content.

Mathematical formulation
This entire study considers the dynamics of small-amplitude perturbations developing in the basic pulsatile channel and pipe flows specified in the previous section. The incompressible Navier-Stokes equations are, therefore, linearised about these base flows.
Considering that the base flows do not depend on the streamwise coordinate x 1 nor on the spanwise/azimuthal coordinate x 2 , infinitesimally small velocity and pressure perturbations may thus be written as spatial normal modes of the form where α 1 and α 2 are streamwise and spanwise/azimuthal wavenumbers, respectively. Separation of total flow fields into basic and perturbation quantities and substitution of the expansions (4.1) and (4.2) into the governing equations (2.1) and (2.2) linearised about the relevant time-periodic base flow then yields a system of coupled linear partial differential governing equations of the form Here, the superscript d refers to components of the direct problem, to be distinguished from the adjoint variables below (4.6). The spatial differential operator L(x 0 , t) in (4.3) is a 4-by-4 matrix and its coefficients involve ∂ 0 -differentiation, depend on the wavenumbers α 1 and α 2 as well as on the base velocity profiles U 1 (x 0 , t); see Appendix C for explicit expressions of all these terms. When studying transient growth effects and searching for optimal initial perturbations that are maximally amplified over a finite-time horizon, it is necessary to choose an appropriate measure of disturbance size (Schmid 2007). Using a classical energy-based inner product, the adjoint governing equations associated with the direct problem (4.3) are routinely obtained as (4.5) and the adjoint differential operator L † (x 0 , t) is explicitly given in Appendix D. In contrast with the direct (4.3), the adjoint equations (4.5) have negative diffusion coefficients and the adjoint fields are integrated backwards in time.
Denoting by {q(x 0 , t i ), |x 0 | < D/2} an initial perturbation at time t i , the evolution of this disturbance at subsequent times t > t i and the associated perturbation energy E(t) are then obtained by solving the initial-value problem corresponding to (4.3) with q(x 0 , t i ) specified for |x 0 | < D/2. The temporal evolution of the perturbation amplitude is then characterised by the ratio The maximum possible amplification of a disturbance over the interval t i < t < t f is obtained as by optimising over all possible initial conditions at t = t i . Note that, since the base flow is time-periodic, the amplification factor depends not only on the duration t f − t i of the temporal evolution but also on the phase of its starting point t i within the pulsation cycle. The particular initial condition at t = t i that achieves the largest amplification at t = t f is referred to as the optimal perturbation and the resulting flow fields at t = t f as the optimal response. In practice, the amplification factors G(t i , t f ) and associated optimal perturbations and responses are iteratively computed by successive direct-adjoint loops, consisting of temporal integration of the direct (4.3) from t i to t f and of the adjoint equations (4.5) from t f to t i , using the numerical methods described in the next section.
In linearly stable configurations, all perturbations eventually decay and the maximal transient growth for given wavenumbers α 1 and α 2 , is well defined and takes finite values. Obviously, G max (α 1 , α 2 ) also depends on the base flow configuration and its control parameters. For a given pulsating base flow, the largest possible transient amplification that may be achieved is obtained as by considering all possible wavenumbers.

Numerical implementation
The direct and adjoint temporal evolution problems (4.3) and (4.5) are first order in time and involve spatial differential operators in the wall-normal x 0 -coordinate. For spatial discretisation we use a Chebyshev spectral method with collocation points spanning the whole diameter of the channel or the pipe. Whether considering channel or pipe flows, all computations are restricted to half of the domain, 0 ≤ x 0 ≤ D/2, by taking into account the symmetry or antisymmetry of the different flow fields and using the associated discretised differential operators of corresponding symmetry. For channel flow calculations carried out in Cartesian coordinates, the parity of the different flow fields depends on the sinuous or varicose nature of the perturbation under consideration. Note that for all the channel flow configurations considered in this paper, the dynamics is dominated by sinuous perturbations. For pipe flow calculations carried out in cylindrical coordinates, it is the value (even or odd) of the azimuthal mode number that determines the parity of each of the different flow fields. Note that the singularities in the differential operators at the pipe axis (x 0 = 0) are only 'apparent' (Boyd 2001): the exact solution is analytic at the axis even though the coefficients of the differential equations are not. Thus a consistent implementation of the symmetry/antisymmetry conditions at the axis removes any apparent singularities and guarantees that the spectral method yields smooth solutions.
Time-marching of the direct and adjoint incompressible Navier-Stokes equations uses a second-order accurate predictor-corrector fractional-step method, derived from Raspo et al. (2002). In classical fashion, the maximal gain G(t i , t f ), together with optimal initial perturbation and final response, is then obtained by direct-adjoint loops, maximising the energy growth from t = t i to t = t f . All subsequent quantities G max and G max max are derived from the gain G, by maximising over t i and t f , and over α 1 and α 2 .
Resorting to the general formulation of the governing equations detailed in the Appendix A and taking advantage of the relevant symmetry properties of the different flow fields thus leads to a numerical implementation capable of handling all configurations of the present investigation.
This entire numerical solution procedure is a generalisation of an approach already used in our previous investigation (Pier & Schmid 2017), and its implementation in C++ is based on the 'home-spun' PackstaB library (Pier 2015, § A.6). The interested reader is referred to these references for further details of the general method.

Pulsating channel flow
The objective of the present section, which is the core part of the paper, is to investigate how the well known transient-growth properties of steady channel flow are modified by the presence of a pulsating base flow component. Starting with a steady Poiseuille flow, the approach consists of studying the influence of pulsation as the amplitudeQ is increased from 0 for different values of the Womersley number Wo.
First, we consider the growth rates G of streamwise-invariant and spanwise-periodic streaks since they display the largest transient growth for Poiseuille flow. Then, the strikingly different behaviour observed for two-dimensional (spanwise-invariant) flows calls for a systematic computation of all possible three-dimensional perturbations. Having established the optimal amplification rates G max that prevail over the whole wavenumber plane, we are then in a position to derive the maximal achievable energy amplification G max max for a given pulsating base flow and to document its dependence on the pulsation amplitudeQ, the Womersley number Wo and the Reynolds number Re. Finally, a detailed discussion of the energy transfer mechanisms allows us to highlight the various growth mechanisms that come into play during the different stages of the evolution and to explain the huge growth factors that are observed for pulsating flows, already for moderate pulsation amplitudes. We recall that sinuous perturbations prevail for all the situations investigated here; thus, all the results presented in this section correspond to flow fields of sinuous symmetry.
The vast parameter space of the problem requires a systematic exploration of the flow physics and a concentration on essential characteristics by a progressive compression of the governing parameters. To this end, we successively investigate the growth of streaks, two-dimensional and three-dimensional disturbances, before focusing on transient energy growth and the energy transfer mechanisms that accompany the observed amplifications. We conclude by isolating the shape and dynamics of the two-and three-dimensional structures that optimally exploit the unsteady background flow and thus exhibit maximal energy growth. Along this analysis, we present a sequential reduction of the parameter space, starting from the effect of cycle length and cycle phase, via spatial scales to the time horizon for optimal growth. Within each step, the essential features of the transient instability will be presented, before reducing the parameter dependency for the subsequent analysis. This section then culminates in the detailed examination of the most amplified disturbances, for the two-and three-dimensional case, under the influence of a pulsatile background flow.

Growth of streaks
In steady channel Poiseuille flow, largest transient growth is known to occur for initial conditions which are spanwise periodic and consist of streamwise aligned vortices, thus corresponding to perturbations with α 1 = 0 and α 2 / = 0. Figure 1(a) shows the optimal transient amplification at Re = 1000, 2000, . . . , 5000 computed for α 2 = 4, which is near the most transiently amplified spanwise wavenumber. (Throughout this paper, length scales are non-dimensionalised with respect to the channel (or pipe) diameter D.) For a steady base flow, the energy growth factor G(t i , t f ) only depends on the duration t f − t i , here measured in mean-flow advection units τ Q ≡ D 2 /Q (0) . Replotting these data for G/Re 2 and measuring the duration t f − t i in diffusion units τ ν ≡ D 2 /ν = Reτ Q , the curves in figure 1(b) confirm the known scaling laws, leading to a maximum transient growth of G max 1.1 × 10 −4 Re 2 at t max /τ Q 1.9 × 10 −2 Re.
Adding to this steady base flow a pulsatile component of given amplitude and frequency, the transient growth properties are characterised by G(t i , t f ) which then depends both on the phase of the initial perturbation t i within the pulsation period T ≡ 2π/Ω and on the duration t f − t i of the temporal evolution. For Re = 2000 and Re = 5000, figure 2 shows plots of the growth factors G(t i , t f ) for pulsation amplitudesQ = 0.4 and 1.0 at Wo = 10. It is found that the amplitudeQ of the base flow modulation only weakly influences the streak growth. Even increasingQ to values larger than unity (corresponding to negative flow rates during part of the pulsation cycle), does not significantly alter the distribution of G(t i , t f ): the maximum amplification remains at the same level and the growth hardly depends on the phase t i /T. Thus streamwise-invariant (α 1 = 0) perturbations appear to be almost unaffected by the time-dependent component of the base flow and to display a dynamics predominantly dictated by the time-averaged base flow. The discussion of energy transfer mechanisms in § 6.5 below will shed further light on this observation. Comparing figure 2(a) with 2(c), and 2(b) with 2(d), the similarity observed between plots at different Re and sameQ also indicates that the scaling of G with Re 2 remains valid for the transient growth of streaks in pulsating base flows.   6.2. Growth of two-dimensional perturbations Two-dimensional spanwise invariant perturbations, corresponding to α 2 = 0 and α 1 / = 0, exhibit much weaker transient amplification than streaks for the same steady Poiseuille flow. Figure 3 plots the transient growth properties prevailing for Poiseuille flow at Re = 1000, 2000, . . . , 5000 for perturbations with α 1 = 2 and α 2 = 0, near the most unstable two-dimensional perturbation. Here, the maximal amplification G max scales linearly with the Reynolds number and reaches much lower values than those corresponding to streaks (see figure 1a); note that this maximal amplification is also reached for a much shorter time horizon. The evolution of two-dimensional transient growth properties as the amplitudeQ of the pulsating component is increased is given in figure 4. After increasingQ, a second maximum emerges in the plot of G, located around t i /T = 0.2 and (t f − t i )/T = 0.5. In contrast with the situation prevailing for streaks, this second maximum is seen to rapidly grow withQ and to become the dominant feature, here already forQ 0.1. While streaks display much larger transient growth for steady Poiseuille flow, these two-dimensional perturbations are found to become the most efficient optimal perturbations for pulsatile base flows, beyond some threshold value of the pulsating amplitudeQ. This overwhelming growth of two-dimensional perturbations for pulsatile conditions will be explained in § 6.5, below, by detailed monitoring of the amplification process in comparison with the dynamics of temporal Floquet eigenmodes.

Growth of three-dimensional perturbations
The very different transient growth behaviour observed for streaky and two-dimensional perturbations calls for a systematic investigation in the entire (α 1 , α 2 )-wavenumber plane. For a given pulsating base flow, the optimal energy amplification G max (α 1 , α 2 ) (4.8) is obtained by maximising the transient growth G(t i , t f ; α 1 , α 2 ) over all values of t i and t f at each prescribed wavenumber. We have systematically explored the control-parameter space spanning the ranges 1000 ≤ Re ≤ 5000, 5 ≤ Wo ≤ 15 and 0 ≤Q ≤ 1, and a few characteristic results are presented below.
The plot of G max for steady Poiseuille flow (Q = 0) at Re = 4000 (figure 5a) confirms that strongest transient growth occurs for streaks (with α 1 = 0) and that the largest value of G max 1763 is reached at α 2 4.09 (indicated by a black dot). Two-dimensional perturbations (with α 2 = 0) experience growth factors that are two orders of magnitude smaller, with G max 30 for α 1 = 3.1.
The distribution of maximal amplification factors G max in the (α 1 , α 2 )-plane evolves significantly as the amplitudeQ of the pulsating component is increased for a given pulsation frequency. Figure 5(b-g) reveal that, asQ is increased, the maximum energy growth (indicated by a black dot) rapidly switches over from streaks to two-dimensional perturbations that experience growth factors sharply increasing withQ while those    5h), the isolines of G max display a similar structure as for Re = 4000 (figure 5a) but with lower levels. After increasing the amplitudeQ of the pulsating flow component at Wo = 10, figures 5(i,j) show that two-dimensional perturbations again eventually dominate the response. However, at this lower Reynolds number, a larger value ofQ is required for the two-dimensional perturbations to emerge, and the increase of G max withQ also occurs at a lower rate. Thus G max is found to reach values of the order of 10 5 at Re = 2000 for Q = 0.5 and Wo = 10 (figure 5j), while at Re = 4000 values in excess of 10 11 are observed (figure 5e).

Maximal transient growth
The maximal transient energy amplification achievable for a given base flow has been defined as G max max (4.9) and is derived by maximising G max (α 1 , α 2 ) over the entire wavenumber plane.  Figure 6 plots the evolution of G max max as the pulsation amplitudeQ is continuously increased for Womersley and Reynolds numbers in the range 5 ≤ Wo ≤ 15 and 1000 ≤ Re ≤ 5000, respectively. At low values ofQ, the pulsating flow component has a very weak influence and G max max remains near the value prevailing for steady Poiseuille flow at the same Reynolds number. For these low pulsation amplitudes, the optimal initial perturbation corresponds to streaks (with α 1 = 0) and the associated growth duration t f − t i remains very close to that prevailing for the equivalent mean Poiseuille flow (see also figure 8, below). Beyond some critical value ofQ, the amplification factor G max max starts to increase exponentially withQ, as illustrated by the nearly constant slopes in the logarithmic plots of figure 6 (note the different vertical scale used in figure 6a for Re = 1000). This critical valueQ c of the pulsation amplitude depends on Wo and Re as shown in figure 7(a): increasing the Reynolds number is found to promote the two-dimensional perturbations which become the dominant feature already forQ > 0.1 around Re = 5000. ForQ >Q c , the rate of the exponential growth of G max max withQ corresponds to the slopes seen in figure 6 and significantly increases as the Womersley number decreases. As a result, G max max rapidly reaches 'astronomical' values, several orders of magnitude beyond the amplification rates prevailing for the corresponding steady Poiseuille flows. These exponential rates have been computed as and their variation with Re and Wo is given in figure 7(b). In this plot the values of κ have been computed by taking the average overQ c <Q <Q c + 0.1, but note that the growth rate remains nearly constant over much larger intervals ofQ in the two-dimensional regime. Obviously, the growth rates are enhanced with the Reynolds number and they also significantly increase towards the lower Womersley numbers, corresponding to longer pulsation periods. The regime change in the transient growth behaviour occurring forQ =Q c is further illustrated in figure 8 at Re = 4000. The evolution with pulsation amplitudẽ Q of the streamwise α 1 and spanwise α 2 wavenumbers associated with the maximally amplified optimal perturbations display a sharp transition from streaky (α 1 = 0, α 2 / = 0) to two-dimensional (α 1 / = 0, α 2 = 0) perturbations. For Wo = 6 and 8, the spanwise wavenumber here directly switches from α 2 4 to 0. For higher values of Wo, however, a small range inQ is observed where the maximally amplified perturbations consist of oblique waves with small but finite values of α 2 . This corresponds to configurations where the amplification of two-dimensional perturbations is still in competition with streaks, so that the maximum of G max in the (α 1 , α 2 )-plane occurs slightly off the α 1 -axis, as illustrated by the black dot in figure 5( f ) and corresponding dots in figure 8(a,b). It is then only for higher values ofQ that purely two-dimensional (α 2 = 0) perturbations prevail. The transition from streaky to two-dimensional maximally amplified perturbations is also accompanied by a significant change in the duration of the growth phase t f − t i , shown in figure 8(c) in mean-flow advection units τ Q and in figure 8(d) in units of the pulsation period T. These two time scales are associated with different dynamical features and related as Wo 2 T = (π/2)Reτ Q . At weak pulsation amplitudesQ, the duration t f − t i remains very close to the value prevailing for streaks developing in the equivalent steady Poiseuille flow, here approximately 75τ Q at Re = 4000 (compare with figure 1). At higher pulsation amplitudes, when two-dimensional perturbations dominate, maximal amplification occurs over intervals t f − t i that approximately correspond to half a pulsation period, T/2. Thus, the transition from streaky to two-dimensional perturbations also coincides with a change in the dynamical time scale: from streak growth essentially dictated by the mean flow to two-dimensional perturbations strongly amplified over half a pulsation cycle.

Discussion of energy transfer mechanisms
In this final subsection on channel flows, we investigate the energy production and dissipation mechanisms in order to explain the different transient-growth scenarios that have been identified.
Following the notations introduced in § 4, we consider a perturbation of the form using a complex-valued three-dimensional velocity vector u(x 0 , t). Such a perturbation is associated with an instantaneous kinetic energy per unit volume of represents the local energy density. Thus, the temporal energy variation, follows from the dynamics of u(x 0 , t), governed by the Navier-Stokes equations linearised about the pulsating base flow (4.3). Separating terms due to interaction with the base flow from those involving viscous dissipation leads to where The term π(x 0 , t) accounts for energy transfer between the pulsating base flow and the perturbation: it essentially represents energy production due to base-flow shear, but negative values may occur and its profile across the channel crucially depends on the relative phases of u 0 (x 0 , t) and u 1 (x 0 , t).
Another quantity of interest is the instantaneous growth rate particularly relevant during phases of near-exponential amplification. Close monitoring of the spatiotemporal development of the base-flow interaction π(x 0 , t) and the dissipation θ(x 0 , t) terms will clarify the amplification mechanisms that govern the different stages of the dynamics.
We focus on two characteristic configurations that have already been discussed: pulsating base flows at Re = 4000 and Wo = 10 with two different pulsation amplitudes, Q = 0.1 andQ = 0.2, associated with streaky and two-dimensional maximally amplified perturbations, respectively.

Streaky maximally amplified optimal perturbation
For the lower pulsation amplitude ofQ = 0.1, a maximal amplification of G max max = 1.77 × 10 3 is achieved from t i = 0.166T to t f = 1.389T for streamwise invariant and spanwise periodic perturbations with α 1 = 0 and α 2 = 4.073. The associated temporal evolution of the perturbation energy E(t) is shown in figure 9(a), with the corresponding instantaneous growth rate σ (t) in figure 9(b). Here, the transient growth is seen to follow the classical pattern prevailing for steady Poiseuille flow: a strong and very short initial boost for t i < t < t = 0.175T (blue parts of the curves), followed by a phase of gradually weakening growth for t < t < t f (in red) towards the maximum response. And indeed, these curves in figure 9(a,b) are almost identical to the accompanying insets that correspond to the maximally amplified perturbations for steady Poiseuille flow at the same Reynolds number, characterised by α 1 = 0, α 2 = 4.088, G max max = 1.76 × 10 3 . This evolution is the result of energy production Π(t) and dissipation Θ(t), shown in figure 9(c). As can be seen by plotting these quantities relative to the instantaneous energy in figure 9(d), viscous dissipation plays here a minor part in the transient growth throughout the entire process from t i to t f .
The temporal evolution of the spatial structure of the maximally amplified streaky perturbation is illustrated in figure 10. Selected snapshots correspond to the thick black dots in figure 9: t i = 0.166T optimal initial perturbation (thick blue curves); t = 0.169T (thin blue curves); t = 0.175T at maximal instantaneous growth (thick green curves); t = 0.500T (thin red curves); t f = 1.389T optimal response (thick red curves). In order to enable comparison of these profiles throughout the temporal evolution, they have here all been normalised to unit total energy. As expected, the initial perturbation consists in streamwise aligned vortices, that fill the entire channel cross-section, with a vanishing streamwise velocity component: see thick blue curves in figure 10(a-c) and corresponding vector plot in figure 10(e). Transient amplification promotes streamwise velocity while reducing wall-normal and spanwise velocity components, leading to a final response that solely consists of streamwise velocity: see thick red curves in figure 10(a-c) and u 1 -isolines in figure 10( f ). The energy production profiles π shown in figure 10(d) result from the interaction of base flow shear with u 0 and u 1 , and are therefore significant only around t = 0.175T (green curve), while displaying vanishing levels near t i and t f . Dissipation profiles θ (not shown) remain at small values throughout the entire evolution.
Clearly, in this regime, the oscillating component of the base flow has a very weak influence, the amplification process operates as for the equivalent steady Poiseuille flow by redistributing streamwise momentum by streamwise vortices, and the dynamics is essentially dictated by the lift-up effect. This insensitivity to the pulsating base-flow component explains why the maximal growth factors G max max prevailing for streaky optimal perturbations remain at nearly constant level asQ is increased, as observed in figure 6. 6.5.2. Two-dimensional maximally amplified optimal perturbation A markedly different scenario prevails at higher pulsation amplitudes when the largest transient amplification is achieved for two-dimensional (streamwise periodic and spanwise invariant) perturbations.
As already shown in figure 5(d), for a pulsation amplitude ofQ = 0.2, the maximally amplified optimal initial perturbation at Re = 4000 and Wo = 10 occurs for α 1 = 2.619 Spatial profiles of flow fields (normalised to unit total energy) over half-channel 0 ≤ x 0 ≤ D/2 at different snapshots: t i = 0.166T optimal initial perturbation (thick blue lines); t = 0.169T (thin blue lines); t = 0.175T at maximum growth rate (thick green lines); t = 0.500T (thin red lines); t f = 1.389T optimal response (thick red lines). Envelope of (a) wall-normal |u 0 (x 0 , t)|, (b) streamwise |u 1 (x 0 , t)| and (c) spanwise |u 2 (x 0 , t)| velocity perturbations. (d) Energy production π(x 0 , t) terms. Snapshots of velocity fields in half-channel over two spanwise wavelengths (λ 2 = 2π/α 2 ): (e) vector plot of (u 0 , u 2 ) for initial perturbation at t i and ( f ) equispaced isolines of streamwise component u 1 of response at t = t f .  168T to t f = 0.698T, leading to G max max = 5.48 × 10 4 . (b) Corresponding instantaneous growth rate σ (t). (c) Energy production Π(t) (solid line) and dissipation Θ(t) (dashed line). (d) Production and dissipation terms relative to instantaneous energy. Growth occurs in two stages: phase I (in blue) for t i < t < t and phase II (in red) for t < t < t f , separated by stall at t = 0.251T (green vertical line). and α 2 = 0 and leads to an amplification of G max max = 5.48 × 10 4 from t i = 0.168T to t f = 0.698T. The temporal evolution of the perturbation energy is given in figure 11(a), with the associated instantaneous growth rate in figure 11(b). The transient growth that occurs over the interval t i < t < t f is here seen to develop in two stages: first a relatively short period (phase I, blue curves) of rapid growth followed by a longer interval (phase II, red curves) of weaker but almost constant growth. Between these two stages, the amplification stalls and the instantaneous growth displays a minimum value, which is here slightly negative around t = 0.251T. This two-stage evolution results from production and dissipation contributions, as illustrated in figure 11(c,d): a first peak in Π(t) during phase I is responsible for the rapid growth of the perturbation, followed by a sustained nearly exponential increase during phase II. The contribution of the relative dissipation Θ(t)/2E(t) is significant only at the very beginning, before rapidly dropping to low values.
The mechanisms responsible for the growth of the perturbation differ in both phases, as illustrated by the profiles in figure 12. These plots show the evolution of the spatial distribution of various fields by selected snapshots, corresponding to the thick black dots in figure 11. Note that the perturbation profiles have again been normalised to unit total energy. The envelopes of wall-normal |u 0 | and streamwise |u 1 | velocity components in figure 12(a,b) show that the initial perturbation at t i (thick blue curves) is localised toward the wall with maximum amplitude around x 0 = 0.4D, before spreading over the entire channel cross-section in the subsequent evolution. The spatial distribution of the base-flow π θ x 1 /λ 1 Figure 12. Evolution of spatial structure of maximally amplified (two-dimensional) optimal perturbation for Re = 4000, Wo = 10 andQ = 0.2. Spatial profiles of flow fields (normalised to unit total energy) over half-channel 0 ≤ x 0 ≤ D/2 at different snapshots: t i = 0.168T optimal initial perturbation (thick blue lines); t = 0.186T and t = 0.210T in phase I (thin blue lines); t = 0.251T at stall (thick green lines); t = 0.300T and t = 0.500T in phase II (thin red lines); t f = 0.698T optimal response (thick red lines). (a) Envelope of wall-normal |u 0 (x 0 , t)| and (b) streamwise |u 1 (x 0 , t)| velocity perturbations. (c) Energy production π(x 0 , t) and (d) dissipation θ(x 0 , t) terms. Snapshots of (u 0 , u 1 ) velocity fields in half-channel over two streamwise wavelengths (λ 1 = 2π/α 1 ): (e) initial perturbation at t = t i and ( f ) response at t = t f .
interaction terms π(x 0 , t) (shown in figure 12c) reveals that the driving mechanism is strong and concentrated around x 0 = 0.4D in phase I (blue curves) while weaker and evenly spread out in phase II (red curves). In contrast, plots of θ(x 0 , t) (figure 12d) show that dissipation is only significant in the initial stages for t t i and reduced to a very thin boundary layer near the wall throughout the rest of the evolution. The vector plot of the initial perturbation in a streamwise channel cross-section (figure 12e) highlights the flow structures concentrated near the wall and characteristically tilted upstream. In contrast, the final response (figure 12f ) fills the entire channel.  Figure 13. Temporal energy evolution of maximally amplified perturbation (blue and red curves) compared with least stable normal mode (thick black curve), at Re = 4000, Wo = 10 andQ = 0.2, for α 1 = 2.619. Normal-mode energy exponentially decays in the long term, according to a Floquet multiplier of μ F = 0.0766, but displays intracyclic amplification by a factor of G nm = 4.30 × 10 3 , approximately during half a pulsation period. Optimal perturbation is amplified by G max max = 5.48 × 10 4 . Phase II (in red) closely follows normal mode while phase I (in blue) is very similar to optimal growth prevailing for steady Poiseuille flow at same parameters (inset). Grey sinusoidal curve represents base flow rate Q(t) (not to scale).
Finally, we compare the dynamics of the present maximally amplified optimal perturbation with the development of the temporal normal mode prevailing for the same pulsating base flow at the same spatial wavenumbers. Such normal modes have been extensively computed and characterised in our previous investigation (Pier & Schmid 2017), using both Floquet eigenproblems and linearised direct numerical simulations.
All pulsatile base flows under consideration here are linearly stable so that temporal eigenmodes decay in the long term. The thick black curve in figure 13 shows the temporal evolution of perturbation energy for the least stable normal mode at Re = 4000, Wo = 10 andQ = 0.2, with α 1 = 2.619 and α 2 = 0. The negative mean slope in this logarithmic plot confirms the decay, governed by a Floquet multiplier of μ F = 0.0766. Thus, the perturbation energy of this normal mode is reduced by a factor of μ 2 F after each pulsation period. However, within each pulsation cycle, significant modulation occurs. This intracyclic growth and decay has been shown to approximately coincide with base-flow deceleration and acceleration phases (Pier & Schmid 2017), as indicated by the grey sinusoidal line representing Q(t). Here, the normal mode displays an intracyclic amplification of G nm = 4.30 × 10 3 .
Comparison of optimal-perturbation and normal-mode energy curves reveals that, during phase II (t < t < t f , red part of curve in figure 13), the optimal perturbation closely follows the normal-mode dynamics. And, indeed, optimal perturbation and normal mode also display very similar flow fields during that interval.
In the initial phase I (t i < t < t , blue part of curve in figure 13), however, the maximally amplified perturbation takes advantage of the optimal initial condition responsible for the initial boost in the response through the Orr mechanism. The amplification during phase I is almost identical to the maximal growth experienced by a two-dimensional optimal initial condition for steady Poiseuille flow at the same Reynolds number and same streamwise wavenumber, shown in the inset in figure 13.
These considerations reveal that the maximally amplified two-dimensional perturbation is an optimal combination of Orr mechanism (phase I) and intracyclic normal-mode growth over half a pulsation cycle (phase II). Growth during phase I is essentially determined by the equivalent steady Poiseuille base flow: the resulting amplification factor therefore scales approximately linearly on Re while being largely independent of Wo and Q. In contrast, growth during phase II closely follows the intracyclic amplification of the associated temporal eigenmodes, and the magnitude of this intracyclic growth has been shown to strongly depend on Wo andQ: whatever the Womersley number, it increases almost exponentially withQ, and the increase is fastest at the lower values of Wo (Pier & Schmid 2017). This exponential growth withQ explains why two-dimensional optimal perturbations always eventually prevail over streaky perturbations, as observed in figure 6.
The maximal growth factor G max max is obviously always larger than either contribution of phase I or of phase II to the total growth. But while the contribution of phase I remains at moderate levels (one or two orders of magnitude, as for steady Poiseuille flow), it is phase II that is responsible for the huge amplification factors prevailing as the modulation amplitudeQ increases. As a result, except for weak pulsation amplitudes, the Orr mechanism only contributes a small factor to the maximal growth G max max , while most of the amplification process is due to modal growth during base-flow deceleration.

Pulsating pipe flow
After having presented detailed results for channel flows, we now turn to the transient growth properties of pulsating flows through circular pipes. The organisation of this section is similar to the previous one. However, since most features are equivalent, many details may be omitted here.
By adopting the general formulation appropriate for both Cartesian and cylindrical coordinates, the analysis of pulsating pipe flows is carried out with the same numerical codes as previously used for pulsating channel flows. Due to periodicity in the azimuthal coordinate, the wavenumber α 2 only takes integer values, but otherwise the numerical implementation proceeds as for a Cartesian formulation. Recall that the apparent singularity at the pipe axis (x 0 = 0) resolves itself by taking advantage of the symmetry properties relevant for each flow component, since all flow fields are either symmetric or antisymmetric in the radial coordinate x 0 .

Transient growth of streaks and helical perturbations
Since for steady Hagen-Poiseuille flow, streamwise invariant streaks with α 2 = 1 and α 1 = 0 undergo the largest non-modal growth, we first consider the transient amplification features prevailing for the same type of perturbations developing in pulsating pipe flows. Figure 14 shows the amplification factors G(t i , t f ) obtained at Wo = 10, for Re = 2000 and 5000,Q = 0.4 and 1.0. The control parameters are the same as those used in figure 2 for pulsating channel flow, and it is observed that the transient growth properties are very similar.
For streamwise periodic (α 1 / = 0) perturbations, the least stable temporal modes correspond to helical perturbations with α 2 = 1. Investigation of transient growth characteristics for α 1 / = 0 also confirms that perturbations with α 2 = 1 dominate over axisymmetric perturbations (α 2 = 0) as well as over those of higher azimuthal order (α 2 ≥ 2). Figure 15 illustrates the transient growth properties for α 2 = 1 and α 1 = 2 at Wo = 10 and Re = 2000 and Re = 5000 as the amplitudeQ of the pulsating base flow component is increased. As for pulsating channel flow, a second maximum emerges that rapidly dominates the dynamics beyond some value of the base flow modulation amplitude. This maximum is again located near t i /T = 0.2 and (t f − t i )/T = 0.5 and corresponds thus to amplification over half a pulsation cycle.

7.2.
Optimal growth at given wavenumbers The maximal transient growth G max , computed by optimisation of G(t i , t f ) over all values of t i and t f for fixed wavenumbers α 1 and α 2 , is shown in figure 16. In each plot the evolution of G max curves for 0 < α 1 < 6 is illustrated asQ is increased fromQ = 0 tõ Q = 1 in steps of 0.1. Panels (a-c) compare the growth of axisymmetric perturbations, α 2 = 0 in panel (a), with that of helical perturbations, α 2 = 1 in panel (b) and α 2 = 2 in panel (c). Clearly, under pulsating flow conditions, axisymmetric initial conditions undergo transient amplification that is not much larger than for the equivalent steady Poiseuille flow, as demonstrated by the nearly overlapping curves in figure 16(a). Non-axisymmetric perturbations, however, experience transient amplification that rapidly grows withQ, and strongest growth occurs for α 2 = 1 (figure 16b). Computation of G max for all α 2 ≤ 6 (results not shown) reveals that the same scenario prevails at higher azimuthal order, but the rate of increase of G max withQ is significantly weaker for higher α 2 .
Evolution of the growth characteristics for Re = 4000 with different Womersley numbers, Wo = 8 in panel (d), Wo = 12 in panel (e) and Wo = 14 in panel ( f ), confirms again that largest amplification factors occur for lower pulsation frequencies, i.e., longer pulsation cycles.
Finally, values obtained at lower Re = 2000 for Wo = 8 in panel (g), 10 in panel (h) and 14 in panel (i) show that the general trend is similar but with lower values of G max , as expected for lower Re.

Maximal amplification
Finally, the maximal amplification G max max achievable for a given pulsating pipe flow is obtained by optimising G max (α 1 , α 2 ) over all streamwise wavenumbers α 1 and azimuthal mode numbers α 2 . Figure 17 shows the variation of G max max as the pulsation amplitudeQ is increased for Womersley numbers in the range 5 ≤ Wo ≤ 15 and Re = 2000, 3000, 4000 and 5000. The behaviour is again found to be similar to that prevailing for pulsating channel flows: at low pulsation amplitudes, G max max hardly departs from the value corresponding to the equivalent steady Poiseuille flow; beyond a critical valueQ c of the pulsation amplitudeQ, transition to approximately exponential growth of G max max with Q takes over. The results shown in figure 17(a) perfectly match those of figure 4(a) of Xu et al. (2021), for the subset of control parameter values that is common to both studies. This agreement further validates our methods.
The variation with Wo and Re of the critical valueQ c for crossover between the two regimes is shown in figure 18(a). In the exponential regime prevailing forQ ≥Q c , the growth rates κ corresponding to the slopes in figure 17 are given in figure 18(b), computed according to (6.1). For a given Reynolds number, the curves ofQ c in figure 18(a) are seen to display a minimum for moderate values of the Womersley number, while they increase both for large and small values of Wo. The increase ofQ c with Wo for Wo ≥ 10 is strongest at lower values of the Reynolds number. In contrast, for Wo ≤ 10 the values ofQ c depend much less on Re.  Comparison of the values ofQ c and κ for pipe flows (figure 18) with those prevailing for channel flows shown in figure 7, reveals that pipe flows require larger pulsation amplitudes to switch to the regime with exponentially growing amplification factors G max max . This is especially true for lower Womersley numbers (see also figure 20 below with additional data for Wo = 3). Also, while the growth rates κ display very similar trends for both channel (figure 7b) and pipe configurations (figure 18b), the values for pipe flows are approximately half those of channel flows.
The streamwise wavenumber α 1 associated with the most amplified perturbation as the pulsation amplitudeQ is varied for a range of Womersley numbers is monitored in figures 19(a) and 19(b) for Re = 2000 and Re = 4000, respectively. These plots demonstrate that the regime change occurring atQ c is indeed associated with a jump in streamwise wavenumber from α 1 = 0 forQ <Q c to finite α 1 -values forQ >Q c . In contrast with channel flows, however, for all pulsating pipe flow configurations investigated here, the optimal perturbations always occur with azimuthal mode number α 2 = 1. Thus the critical valueQ c always corresponds to a transition from streaky (α 1 = 0, α 2 = 1) to helical (α 1 / = 0, α 2 = 1) optimal perturbations, at the same azimuthal mode number. This regime change is also associated with a discontinuity in the duration of the growth phase t f − t i for the optimal amplification process, as illustrated in figure 19(c,d) for Re = 4000. The optimised duration t f − t i is given in mean-flow advection units τ Q in figure 19(c) and in units of the pulsation period T in figure 19(d). These plots illustrate that pulsating pipe flows display similar transient dynamics as channel flows: forQ <Q c , the optimal duration t f − t i remains close to the value prevailing for the average parabolic flow profile; forQ >Q c , when helical perturbations dominate, maximal amplification occurs over intervals corresponding approximately to half a pulsation period. Thus our results confirm the findings of Xu et al. (2021) that helical perturbations dominate the transient growth at large pulsation amplitudes. By our detailed comparison of channel and  pipe configurations at moderate pulsation frequencies, we have been able to highlight the fundamental growth mechanisms, which are common to both geometries.

Conclusion
Considering pulsating flows through both channels and pipes, we have investigated the non-modal transient energy amplification resulting from optimal initial conditions. Our study has systematically covered the pulsating base flows for 1000 ≤ Re ≤ 5000, 5 ≤ Wo ≤ 15 and 0 ≤Q ≤ 1. While channel and pipe flows display quite different linear modal stability characteristics, their non-modal transient growth features are found to be very similar in situations that are linearly stable. Optimal energy growth occurs according to two distinct scenarios. At weak pulsation amplitudesQ, the behaviour is similar to that resulting from the equivalent steady Poiseuille flow, and the oscillating flow component appears to have only a small effect. Beyond a critical value ofQ, however, transient growth increases exponentially withQ and reaches astronomical values, already for moderate pulsation amplitudes. In this latter regime, optimal growth mainly occurs over half a pulsation period, during the slow part of the pulsation cycle, and closely follows the intracyclic amplification of the associated Floquet eigenmodes. We have previously shown (Pier & Schmid 2017) that the intracyclic modulation amplitudes derived from temporal normal modes may be huge, even for linearly decaying eigenmodes. The maximal transient amplification factors G max max computed in the present investigation are even larger since they take advantage of both this normal-mode intracyclic growth and non-modal Orr-type amplification, which contributes in the early stage of the growth process.
These findings have been firmly established by a comprehensive investigation of pulsating flows over the range 5 ≤ Wo ≤ 15, deemed to be the most relevant for applications in the haemodynamic context. In order to explore the expected behaviour beyond that frequency range, figure 20 shows the maximal achievable transient growth G max max for channel and pipe flows at Re = 5000, extending the results of figures 6(d) and 17(d) by including data at lower and higher pulsation frequencies, Wo = 3 and Wo = 20, respectively. At high pulsation frequencies, it is observed that the pulsating component is rather inefficient in producing G max max factors beyond those prevailing for steady base flows, a result closely related to the fact that high-frequency pulsation also has a strong stabilising effect on modal temporal growth rates. In the low frequency regime, the curves for Wo = 3 indicate that strong growth is still possible but requires larger pulsation amplitudesQ. When lowering Wo, the critical valueQ c for onset of the exponential regime increases moderately for channel flows and significantly for pipe flows (Q c 0.98 for Wo = 3). These plots are in agreement with observations already made by Xu et al. (2021) for pulsating pipe flows. Concerning the spatial structure of the optimal perturbations, the results of figure 20 follow the same scenario as previously discussed: while streaky perturbations prevail forQ <Q c , at larger pulsation amplitudes two-dimensional and helical perturbations dominate, respectively, in channel and pipe flows.
It should be noted that a major difference between channel and pipe flows concerns their linear modal instability features. Indeed, for channel flows there exists a critical Reynolds number beyond which linear instability occurs. This is well known for steady Poiseuille channel flow, and the dependence of this critical Reynolds number with the pulsating flow parameters has been extensively discussed in our previous work (Pier & Schmid 2017). By contrast, steady pipe Poiseuille flow remains linearly stable, whatever the Reynolds number. For time-periodic base flows, linear instability has been found for purely oscillating pipe flows (Thomas et al. 2011). However, the presence of a non-vanishing mean flow rate appears to have a stabilising effect and all pulsating pipe flows considered in the present study are far from temporal instability.
Another difference is that two-dimensional (spanwise invariant and streamwise periodic) perturbations are the most unstable or the least stable for channel flows, whereas the leading linear instability for pipe flows occurs for helical modes with α 2 = 1 and α 1 / = 0, which dominate over perturbations of higher azimuthal order (α 2 ≥ 2) as well as over axisymmetric (α 2 = 0) ones. It is found that this remains true for pulsating pipe flows. While channel flows are rapidly dominated by two-dimensional sinuous perturbations, pipe flows are dominated by helical perturbations in similar pulsating flow regimes. For pipe flows, axisymmetric perturbations never prevail. But note that the Cartesian equivalent of axisymmetric perturbations are two-dimensional perturbations of varicose symmetry, which never prevail either. The closest equivalent to a two-dimensional sinuous perturbation in a cylindrical geometry is a helical perturbation (with α 2 = 1).
This study gives a detailed and comprehensive perspective on the perturbation dynamics in pulsatile channel and pipe flow, treating these configurations within a time-dependent, initial-value problem formalism and thus avoiding restrictive assumptions of a modal, time-asymptotic approach. This analysis identified a rich perturbation behaviour driven by parametric and transient excitation over one or multiple forcing cycles and the dominance of an Orr-type amplification mechanism at early times that acts efficiently and selectively across a significant parameter range, once a critical pulsation amplitude has been surpassed.
Our study lays the foundation for a future analysis of pulsating base flows with higher harmonic content, such as blood flow rates resulting from the cardiac pulse. The present approach could also be generalised to take into account compliant walls or to address nonlinear fluid effects.
∂ t u 2 + (u · ∇)u 2 + 1 x 0 u 0 u 2 = ν ⎛ ⎝ Δu 2 + 1 with the notations In these expressions, the terms enclosed in boxes are only present in the formulation using cylindrical coordinates and pertaining to the pipe flow configuration. Resorting to such a general formalism is particularly useful when developing numerical codes to solve both channel and pipe flows: the boxed terms may be switched on or off depending on the flow configuration.
All the profiles above are normalised so that their cross-sectional average equals unity. Thus, the pulsating base flow velocity components (3.4) are simply obtained by multiplying these profiles with the flow rate coefficients Q (n) .

Appendix C. Linear governing equations of direct problem
In the direct formulation of the incompressible Navier-Stokes equations (4.3), the spatial differential operator L(x 0 , t) is a 4-by-4 matrix of the form Its coefficients involve ∂ 0 -differentiation, depend on the spatial wavenumbers as well as on the base flow velocity profiles. Their explicit expressions are the following: L 20 (x 0 , t) = −L 02 (x 0 , t) = 2ν iα 2 with

Appendix D. Adjoint problem
In the adjoint formulation of the incompressible Navier-Stokes equations (4.5), the spatial differential operator L † (x 0 , t) is obtained as