Floquet stability analysis of pulsatile flow in toroidal pipes

The linear temporal stability of the fully-developed pulsatile flow in a torus with high curvature is investigated using Floquet theory. The baseflow is computed via a Newton– Raphson iteration in frequency space to obtain basic states at supercritical Reynolds numbers in the steady state for two curvatures, 𝛿 = 0 . 1 and 𝛿 = 0 . 3, exhibiting structurally different linear instabilities in the steady case. The addition of a pulsatile component is found to be overall stabilising over a wide range of pulsation amplitudes, in particular for the higher values of the curvature. The pulsatile flows are found to be at most transiently stable with large intracyclic growth rate variations even at small pulsation amplitudes. While these growth rates are likely insufficient to trigger abrupt transition at the parameters in this work, the trends indicate that this is indeed likely for higher pulsation amplitudes, similar to pulsatile flow in straight pipes. At the edge of the considered parameter range, subharmonic eigenvalue orbits in the local spectrum of the time-periodic operator, recently found in pulsating channel flow, have been confirmed also for pulsatile flow in toroidal pipes underlining the generality of this phenomenon.


Introduction
Unsteady flows in curved pipes are not only ubiquitous in engineering applications such as heat exchangers (Vashisth et al. 2008), mixers (Jarrahi et al. 2011) and the exhaust manifolds of internal combustion engines (Kalpakli Vester et al. 2016).Also in biological processes the flow can be modelled as unsteady pipe flow, most prominently the pulsatile blood flow from the heart through the aortic arch in the cardio-vascular system (Pedley 2003) and the respiratory flow to the lungs (Tanaka et al. 1998).
While the steady flow in straight pipes is stable to infinitesimal perturbations even to very high Reynolds numbers (Meseguer & Trefethen 2003) and the transition occurring at Re ≈ 2040 (Avila et al. (2011), based on the pipe diameter) is due to transient amplification of finite amplitude disturbances (Schmid & Henningson 1994), the pulsatile flow in the same geometry exhibits both non-modal and modal growth of disturbances during the pulsation cycle which can completely change the optimal transient growth mechanism (Xu et al. 2017(Xu et al. , 2021)).It is well-known that time-dependent, and in particular time-periodic flows, can instantaneously exhibit both strong stabilisation and destabilisation over time, even if the flows are categorised as linearly asymptotically stable (Davis 1976).The rapid growth of perturbations, typically during the deceleration phase of the pulsation cycle, can be so large that non-linear effects are triggered leading to abrupt transition to turbulence which has been extensively documented experimentally for straight pipes (Merkli & Thomann 1975;Hino et al. 1976).This tableau of stability and transition in straight pipes contrasts with the situation in curved tubes, in which transition to turbulence is known since the late 1920s to occur at much higher Reynolds numbers (White & Appleton 1929;Taylor 1929).Recently, careful experiments by Kühnen et al. (2014Kühnen et al. ( , 2015) ) as well as extensive numerical simulations (Noorani et al. 2013;Canton et al. 2016Canton et al. , 2017) ) have been able to resolve some inconsistencies in the literature and establish that almost all curved pipes (  /  =  > 2 × 10 −3 , where   and   are the pipe radius and its radius of curvature at the centreline, respectively) exhibit linear instabilities and that, for  > 0.028 (which includes most pipes of practical interest), the road to turbulence begins with a supercritical Hopf bifurcation (Canton et al. 2020).
For unsteady flow in curved pipes, on the other hand, the results are considerably more sparse.To the authors' knowledge, the only relevant results are by Papageorgiou (1987) who studied the centrifugal linear stability of the oscillatory flow (zero mean) in curved pipes in the limit of high frequencies based on earlier work by Seminara & Hall (1976) as well as a later extension by Shortis & Hall (1999) including a weakly non-linear analysis.These results do not carry over to pulsatile flows with a non-zero mean flow, however, since the instability mechanisms in these two flow regimes differ fundamentally.In fact, oscillatory flows are governed solely by the imposed frequency whereas pulsatile flows also depend on the timescales of the steady flow component.The associated instabilities in pulsatile flow are therefore intimately linked to their steady counterpart.
In order to consider questions of stability, first a basic state must be established.While analytical profiles are readily available for a number of plane oscillatory and pulsatile shear flows (e.g.oscillating Stokes layers (Stokes 1851) and pulsating plane Poiseuille and pipe flows (Sexl 1930;Womersley 1955;Uchida 1956)), the geometry of curved pipes leads to intrinsically two-dimensional flows even in the fully-developed regime which considerably complicates analytical treatment.Although asymptotic solutions for steady flow in toroidal pipes in the limit of low curvatures were derived in the late 1920s (Dean 1927(Dean , 1928)), the first results for considerably more difficult problem of oscillatory flow were only found two decades later by Lyne (1971) via matched asymptotic expansions of the viscous boundary layer and the inviscid core region.Lyne (1971) discovered that for high enough pulsation frequencies, the dipole structure described by Dean (1928), characteristic of steady flow, is replaced by a quadrupole vortex structure.Since these seminal papers, there have been numerous works aimed at verifying the results as well as extending the solutions to more broader parameter ranges.While Zalosh & Nelson (1973) confirmed Lyne's results with an alternative method, Smith (1975) and later Simon et al. (1977) extended the approach to pulsatile flow and Siggers & Waters (2008) was able to relax the requirement for asymptotically small curvatures typical of earlier perturbation approaches.In tandem, a number of experimental campaigns have been able to confirm the theoretical results of Lyne (1971) and Zalosh & Nelson (1973), in particular Bertelsen (1975) and Munson (1975), and add a great deal of detail to the description by mapping out the parameter space (Sumida et al. 1984;Sudo et al. 1992).In spite of this progress, many practical parameter combinations are out of reach for perturbation approaches.The need for accurate flow predictions in both engineering and medicine combined with the increasing availability of computational resources and efficient simulation codes has led to the development of a number of numerical techniques to simulate these complicated flows, mainly based on finite-difference schemes (Lin & Tarbell 1980;Rabadi et al. 1980;Hamakiotes & Berger 1988;Dwyer et al. 2001).Simultaneously, the field has branched out to cover not only unsteady flow in classical toroidal pipes (e.g.Lin & Tarbell (1980); Sudo et al. (1992)) but also includes studies of pulsatile entry flow (Talbot & Gong 1983;Rindt et al. 1991), intermittent pulsations (Komai & Tanishita 1997), turbulent flows (see Kalpakli Vester et al. (2016) for a recent review) and increasingly complicated physiological geometries (e.g.Morris et al. (2005)).
In the context of time-periodic flows, the general temporal linear stability problem can be expressed in terms of Floquet theory where the evolution of perturbations over one pulsation cycle are related via the Floquet operator.The asymptotic stability is then governed by the eigenvalues of this operator called Floquet exponents.There are fundamentally two approaches to the solution of the resulting generalised eigenvalue problem.The first option involves the time-integration of the linear evolution equation over the pulsation period.While this is comparatively inexpensive when the baseflow is given analytically (e.g.Rozhdestvenskii et al. (1989)), the procedure is considerably more involved of the basic state is itself solved numerically.This time-stepper approach has successfully been applied in a number of flows, for example in the cylinder wake (Barkley & Henderson 1996) and in jets (Shaabani-Ardali et al. 2019).While this approach is easily combined with existing (matrix-free) time-stepping frameworks, it has the major drawback that the cost is directly proportional to the pulsation period  ∝ Re/Wo 2 , where Wo is the Womersley number or frequency parameter encapsulating the pulsation time-scale, which can be prohibitively large for low frequencies.The alternative and the method adopted in this work (see also Pier & Schmid (2017)) is to solve the generalised eigenvalue problem directly which is nominally independent of , although, as we will see, there is an intrinsic dependence via the bandwidth of the temporal Fourier expansions of the baseflow and perturbation.
In this work, we study the temporal Floquet stability of the fully-developed pulsatile flow in a torus at two values of the curvature over a wide range of pulsation frequencies and amplitudes.The curvatures of  = 0.1 and  = 0.3, which exhibit structurally different linear instability mechanisms in the steady case, are considered at slightly supercritical (steady) Reynolds numbers in order to assess the effect of added pulsations on the temporal stability of the flow.The baseflows are computed via a novel Newton-Raphson iteration on a Poincaré section of the periodic orbit that is described and validated in detail in Kern (2023).The generalised eigenvalue problems are solved directly using an iterative shift-and-invert strategy to extract only the dominant Floquet exponents and associated eigenfunctions.

Governing equations
We consider the pulsatile, incompressible flow of a Newtonian fluid in a toroidal pipe driven by a time-periodic external forcing governed by the incompressible Navier-Stokes equations

𝜕𝑼 𝜕𝑡
which are expressed in orthogonal toroidal coordinates (Germano 1982;Hüttl et al. 1999).Therefore,  is the velocity vector with streamwise, radial and azimuthal components (  ,   ,   ),  is the reduced pressure and  is a time-periodic external volume force that mimics the effect of a streamwise pressure gradient.The flow is contained in a rigid toroidal pipe which are modelled by imposing a no-slip boundary condition on the pipe walls.The detailed structure of the Navier-Stokes equations can be found in Kern (2023).
The pulsatile flow oscillating at frequency Ω imposed by the volume force is a solution of 2.1.As the geometry of the torus is invariant along the pipe, the resulting flow field is independent of the streamwise coordinate .The velocity components as well as the pressure that are only functions of (, , ) can therefore formally be Fourier transformed in time to yield which are driven by a volume force of the form Note that since the original velocity, pressure and forcing fields are real-valued, we have , where the superscript * denotes complex conjugation.Introducing the Fourier expansions into (2.1), the governing equations for the  th Fourier component can be written as iΩ  () + ∑︁ where we have omitted the explicit dependence on (, ) for notational conciseness.
The equations 2.4 are exact for arbitrary choices of time-periodic forcing functions  .In order to reduce the parameter space, we restrict ourselves to purely sinusoidal forcing, i.e. we choose  () = 0, || ⩾ 2. Due to the non-linearity of the equations, the imposed simple harmonic temporal dependence of the forcing does not carry over to the velocity and pressure fields.Instead, the convective term in 2.4 lead to a broadband response of the flow which requires the retention of more Fourier modes in the computation of the baseflow profiles.Note however that in most cases the amplitudes of the higher harmonic components quickly decay since they are driven exclusively by the intrinsic non-linear interactions of the flow (Kern 2023).

Non-dimensional parameters
The governing equations are non-dimensionalised with the bulk velocity   , the radius of the cross-section of the pipe   and the kinematic viscosity , encapsulated in the Reynolds number In literature it is common practice to use the Reynolds number Re  based on the pipe diameter  = 2  instead of the radius with Re  = 2Re.The pipe geometry is parameterised using the curvature  defined as the ratio of the pipe cross-sectional radius and the radius of curvature of the pipe at the pipe centreline   , i.e.  =   /  .The external volume force imposes a time-scale on the flow via the frequency Ω or, equivalently, the period of pulsation  = 2/Ω.The oscillation frequency is nondimensionalised by the kinematic viscosity  and the pipe cross-sectional radius to produce the Womersley number Wo =   √︁ Ω/ .
(2.6) The last parameter that defines the pulsatile baseflow is the amplitude of the pulsation compared to the steady component.Following the approach in Pier & Schmid (2017); Kern et al. (2021) used for the (linear) problem of pulsating channel flow, we choose to consider the time-periodic streamwise flow rate   given by Since the flow rate is real-valued we have  (0) ∈ R and we can equally choose  (1) ∈ R without loss of generality.Gathering all higher harmonic components in the term  ℎ.ℎ.

𝑠
(), (2.7) can then be rewritten as where  = 2 (1) / (0) is used as a governing parameter to parameterise the solution.A drawback of this choice is that since the amplitude of the higher harmonics is a function of the input parameters, the total flow rate varies over the parameter range although their contribution is negligible in most cases considered in this work.

Linearised stability equations
In order to study the stability of the pulsatile baseflow  = (, ) obtained from 2.4 to small disturbances, we consider an infinitesimal perturbation on top of the baseflow whose evolution is governed by the linearised Navier-Stokes equations given by where  = (  ,   ,   ) is the perturbation velocity field,  is the associated perturbation pressure.The equations use the same coordinate system and non-dimensionalisation as their non-linear counterpart (2.1).Combining all perturbation components into a single vector  as (, , ) = (  ,   ,   , )  , (2.10) where the superscript  denotes matrix transposition, (2.9) can be written concisely in operator form as where R is the singular restriction operator due to the fact that there is no evolution equation for pressure in incompressible flows and L is the full linear operator acting on .For notational convenience, L has been split into a diffusive/pressure part (L 1 ) and a convective part (L 2 ()) according to whether the operator explicitly depends on the baseflow velocity  or not.The linearised equations are complemented with the same boundary conditions as their non-linear counterparts.

Floquet analysis
The linearised equations (2.11) are the appropriate framework for the study of any type of perturbation, irrespective of the temporal variation of , and can be integrated in time starting from a given initial condition.In the special case of a periodic basic state, an alternative approach to study its linear stability utilises the results from Floquet theory.The Floquet formalism guarantees that the perturbations have the same temporal periodicity with frequency Ω as the baseflow, which allows us to expand the perturbation  in a temporal Fourier series (, , ) = ∑︁   ( )  iΩ (2.12) which carries over to each of the vector components.Combining this expansion with a normal mode ansatz assuming exponential behaviour in the homogeneous streamwise direction and in time (the intracyclic temporal variation is fully described by the temporal Fourier expansion), we obtain where the (•) symbol denotes the spatial Fourier transform yielding the streamwise wavenumber .In the following we focus on a single wavenumber.Introducing this ansatz into (2.11)leads to an infinite system of coupled ordinary differential equations where the evolution of the th Fourier component û() of the perturbation can be expressed as (2.14) Note that the coupling introduced by the baseflow modulation between adjacent temporal Fourier modes happens exclusively via the convection term L 2 since it depends directly on the baseflow due to the non-linearity in the convective term of the Navier-Stokes equations and does not alter the coupling between the velocity components within the same mode, i.e. the operators L 1 and L 2 are unchanged.Note that equations (2.14) are an exact description of the asymptotic behaviour of a normal mode on top of the time dependent baseflow  in the limit of long times (when the initial transients have passed).
Stacking the Fourier components of the perturbation into a single (doubly infinite) vector   we can write the Floquet problem in more compact form where R  is the global restriction operator and F is the Floquet operator.The asymptotic behaviour of the linear solution is governed by the complex growth rate  and is often expressed using the Floquet multiplier that encapsulates the evolution of the normal mode over one period of pulsation.If the real temporal growth rate is positive, i.e.Re() > 0, then the Floquet system is linearly unstable.This corresponds to a Floquet multiplier that moves outside of the unit circle (|| > 1).

Perturbation energy growth
In order to measure energy growth, we use the standard energy norm ⟨•, •⟩  which naturally arises from the incompressible Navier-Stokes equations and is computed as the integral over the cross-section of the torus (2.17) where the superscript  denotes the complex transpose.
The Floquet eigenfunctions are exact solutions of (2.11) and the intracyclic variation of the perturbation can therefore be reconstructed for any time  using equation (2.12).The instantaneous energy growth rate  of the perturbation  can then be computed as where we have taken advantage of the linearity of the inner product, used the notation ⟨, ⟩  =  and dropped the explicit dependence on time for notational convenience.
Realising that  = R via the definition of the restriction operator, we can write the perturbation energy in terms of the full vector  and use (2.11) to write which is the Rayleigh quotient of the linear operator with ⟨R, R⟩  = ⟨, ⟩  .
The variation of the perturbation energy  (), relative to an initial amplitude  0 =  (0), is then computed as  () It is clear that the energy growth over one pulsation period  is exp where  0 = (0) and we have used the definition of the Floquet multiplier  relating the perturbation across the period.

Discretisation and solution methodology
3.1.Spatial discretisation and symmetrisation In order to solve equations (2.4) and (2.15) numerically, the differential equations must be discretised in the cross-section of the torus.The same discretisation is employed for both the baseflow computations and the solution of the stability problem.A combined Fourier-Chebyshev collocation technique is used in which the differential equations are solved on a spectral grid with   Chebyshev points in radial direction and   equispaced points in azimuthal direction at which the differential equations are enforced exactly.The spatial resolution used for the unsteady computations at each of the curvatures considered in this work was determined based on a spectral analysis of the steady solution shown in A. The resolution chosen for this study is   ×   = 30 × 110 ( 1 = 6600) for  = 0.1 and   ×   = 30 × 150 for  = 0.3 ( 1 = 9000).The corresponding computational grids used for both linear and non-linear calculations are shown in figure 1 for each curvature.
The two-dimensional arrays containing the values of the discretised velocity components and the pressure at the collocation points are combined into a stacked column vector v (placeholder for either the baseflow (Q) or perturbation (q) fields) such that where  1 =   ×   is the size of each component u  , u  , u  , p and 2 = 4 1 is the size of a full field.  is the number of radial Chebyshev gridpoints and   is the number of equispaced azimuthal gridpoints.Due to the assumed symmetry, the computations are performed on half the grid.
In the process, we replace the continuous differential operators L 1 , L 2 and the restriction operator R with their discrete counterparts denoted L 1 , L 2 and R, respectively, with the following macro-structure where 1  and 0  are the identity and zero matrix of order , respectively, and U indicates the dependence on the (discretised) basic state.In each operator, the first three rows correspond to the components of the momentum equation and the last row corresponds to the discretised continuity equation.For the discretised convective operator L 2 the dependence on the baseflow U is noted explicitly.The details of the individual operators withing L 1 and L 2 are summarised in B. The zero sub-blocks are due to the fact that the pressure is not coupled to itself and is not involved in the computation of the convective terms.
Due to the assumed symmetry of the problem along the equatorial plane of the torus, we can restrict the computational domain to half the cross-section.This reduces the number of points in the azimuthal direction by half to    =   /2 at the cost of solving for symmetric and antisymmetric perturbations separately.The symmetrisation is introduced at the level of the operators L 1 and L 2 before the Floquet operator is assembled.For each operator L  , the symmetrised version L , is computed as where T  ∈ R  ×2 and T ∈ R 2 ×  are symmetrisation operators whose structure is given in C. The same symmetrisation is also applied on the solution fields such that In order to complete the problem, the equations need to be supplemented with the appropriate boundary conditions on the perturbations.In this case, we impose no-slip boundary conditions for the velocity at the walls ( =   ) and symmetry conditions at the centreline for all components.The pressure has no boundary conditions at the wall.The no-slip Dirichlet condition is directly imposed on the wall points on the full operators equations (3.2)-(3.3)after symmetrisation by adjusting the first   rows of each velocity component.

Temporal discretisation of the Floquet equations
The practical solution of the generalised eigenvalue problem (2.15) requires the truncation of the temporal Fourier expansion of the eigenfunctions at some   ∈ N leading to a Fourier expansion with   = 2  + 1 modes.
Combining the variables into a column vector such that all   Fourier components of  are gathered and combined into a long vector q of length  0 =  ×   we obtain The discretisation of equation (2.15) leads to a generalised eigenvalue problem for the eigenvalue  of the form  R  q = F   q , (3.8) The discretised Floquet operator F   ∈ C  0 ×  0 has the following macro-structure where P = Ω • diag(−  , . . .,   ) is the diagonal matrix of complex temporal Fourier frequency shifts,

𝑘
is the shift matrix of order  with ones on the  th sub/super-diagonal and ⊗ is the Kronecker product.
The method can deal with arbitrary baseflow Fourier expansions but the system size grows rapidly with   .In order to resolve all baseflow interactions the temporal Fourier expansion must be wide enough, i.e. the number of temporal Fourier modes must be chosen such that   >   .In strongly coupled systems, converged results typically require   ≫   , especially for low values of the Womersley number.
The boundary conditions for each Fourier component of the Floquet eigenvalue problem carry over from the isolated linearised Navier-Stokes problem and are implemented in the same fashion.
The size of the Floquet operator as well as the degree of coupling is determined solely by the width of the baseflow Fourier expansion that is retained for the Floquet analysis.The use of spectral differentiation leads to high accuracy but also reduces the overall sparsity of the derivative operators.All matrices used in this study are constructed in sparse format only saving the non-zero elements.The effect of truncating the baseflow Fourier expansion is considered in section 5.4.The bottleneck for the computation of the spectra is the explicit construction, storage and decomposition of the Floquet matrix F   .Since we are not interested in the full spectrum but instead in the effect of pulsation on the most unstable eigenvalue or at most a few of the least stable modes, we employ a shift and invert algorithm to compute a small number of eigenvalues of equation (3.8) close to the dominant eigenvalue of the corresponding steady problem using the Matlab ® built-in function eigs.The tolerance for the iterative eigenvalue computation was uniformly set to 10 −9 .
Due to the quickly increasing problem sizes with   , the limit for this study is   = 19 at  = 0.1 and   = 13 at  = 0.3 which corresponds to Floquet matrices of order  0 = 138 600 (240 × 10 6 nonzero elements) and  0 = 117 000 (250 × 10 6 nonzero elements), respectively.For the largest problems the limiting factor is the available RAM on the clusters as the decomposition of the shifted Floquet matrix needed for the efficient shift-and-invert algorithm required more up to 1 TB of memory.If the number of retained Fourier modes needs to be further increased, this restriction could be circumvented by using a matrix-free approach to compute the matrix-vector products for the Arnoldi-iteration within the shiftand-invert algorithm, which could then be combined with appropriate preconditioning, but these extensions are beyond the scope of this work.

Baseflow computation
The laminar pulsating baseflows used for the stability analysis in this work are computed using a Newton-Raphson fixed-point iteration in Fourier space.For the sake of completeness, we sketch the solution procedure here.
The truncation of the symmetrised baseflow Fourier expansion at  =   leads to an expansion involving   = 2  + 1 fields of type (3.6) that are stacked into a long vector q of length  =  ×   yielding q = q (−  )  , . . ., q (0)  , . . ., q (  ) At the th iteration starting from the current solution q  , the Newton-Raphson method consists of solving a system of linear equations for the update of the current solution Δq  = q  − q −1 using the linearisation of the governing equations (2.9) around q  which can be written very compactly as where N (q  ) is a shorthand for the discretised non-linear Navier-Stokes equations (2.1) applied to the solution at iteration  and F   is the same linear Floquet operator as in (3.8) constructed for the   Fourier components of the baseflow due to the structural similarity of the Fourier expansions of the baseflow and the linear perturbation.This process is repeated to convergence or stopped at divergence.More details on the solution method can be found in Kern (2023).
For a given combination of input parameters (, Re,Wo, ), the number of Newton steps required for convergence is dependent on the quality of the initial guess and the chosen tolerance.Simultaneously, the cost per Newton iteration quickly increases with the number   of Fourier components retained in the baseflow computation which determine the size of the problem (3.11).Since the spatial structure of the pulsating flows varies considerably more with the Womersley number than with the pulsation amplitude (Kern 2023), the computations are organised as follows.
First, starting from the analytical solution for the steady flow in a straight pipe ( = 0), the steady state solutions (  = 0) are computed for the target Reynolds number of Re = 1700 and curvatures  = 0.1 and 0.3.Next, for each Womersley number, the Newton-Raphson method is initialised with the corresponding steady baseflow solution at the same curvature and Reynolds number expanded to include one unsteady Fourier component.This reduced system is solved for each pulsation amplitude sequentially until either the next higher amplitude cannot be converged or the maximum pulsation amplitude of  = 0.5 is reached.The solutions are then expanded to include the next Fourier component (  = 2) and the updated system is converged to tolerance.This process is repeated until a maximum of   = 3 unsteady Fourier components are included in the computations.Based on the spatial resolution chosen for this work, the linear problems at each Newton step at this temporal resolution involve matrices of order  = 46 000 and 63 000 for  = 0.1 and  = 0.3, respectively.

Pulsating baseflow
For both considered parameter combinations (Re, ) = (1700, 0.1) and (1700, 0.3) the steady solution of the flow in a toroidal pipe is supercritical and the spectrum of the linear operator at the critical streamwise wavenumber exhibits a single unstable eigenvalue that is either symmetric ( = 0.1,  = 1.7) or antisymmetric ( = 0.3,  = 2.1). Figure 2 shows for each case the corresponding steady baseflow profiles as well as the spatial structure of the unstable mode and the spectrum of the linear operator at the most streamwise wavenumber  at which the steady flow first becomes linearly unstable.The spectra also include the variation of the leading eigenvalue for variations of  around the respective critical value showing that at the considered mildly supercritical parameter combination the most dangerous streamwise wavenumber remains unchanged.For a direct comparison with Canton et al. (2016) note that the reference uses the diameter-based Reynolds number and the present results need to be multiplied by Re  /Re = 2.
Starting from the steady solutions, the pulsatile solutions are computed in the parameter range (Wo, ) ∈ [5, 100] × [0, 0.5].The highest considered pulsation amplitudes of  = 0.5 could not be reached for all Womersley numbers because the Newton-Raphson iteration did not converge.This could be due to the stiffness of the problem related to the strong nonlinearity of the Navier-Stokes equations for these parameters for which the present approach of using reduced order solutions as stepping stones to reach the final temporal resolution may be inappropriate.Figure 3 shows a map of the coarse sweep of the parameter space indicating which configurations were computed as well as the parameter region (shaded in red) where the baseflow solutions could not be converged.
Due to the truncation of the Fourier series for the baseflow at   = 3, not all baseflow solutions that have been converged to tolerance can be considered accurate.In fact, the width of the Fourier spectrum as well as the decay rate of the higher harmonics is highly dependent on both pulsation frequency and amplitude.In Kern (2023), the authors computed baseflows at lower Reynolds number (and thus less restrictive spatial resolution) up to   = 6 and showed that for most of the considered Womersley number range the flow rates of the higher harmonics decayed nearly exponentially.Only some cases at intermediate Womersley numbers with strong non-linear coupling exhibited more complicated, nonmonotonic variation of the flow rate with   , especially at higher pulsation amplitudes.In this work, we assume that the amplitude of the flow rate of the higher harmonics is exponential even at the higher Reynolds number considered here and therefore consider a baseflow solution accurate if an exponential extrapolation based on the resolved harmonics predicts that the flow rate   of the next higher harmonic is below the chosen tolerance indicating that increasing the temporal resolution would not alter the solution.Mathematically, we assume that the amplitude of the flow rates of the higher harmonics decays exponentially, i.e.
and consider the case converged if Although this criterion is likely less accurate for intermediate Womersley numbers, we use it for the entire parameter range given that the pulsation amplitudes reached in the highly coupled parameter regions are very low. Figure 4 shows the relative flow rates of the higher harmonics  (2) and  (3) (full and dashed lines, respectively) for specific values of  (colours) over the considered frequency range, where available.The irregular variations of  The baseflow for three Womersley numbers at  = 0.3 and  = 0.01 is shown in figure 5 computed with   = 3.The chosen Womersley numbers are representative of the different flow regimes that exist for the pulsating flow in toroidal pipes.At Wo = 5 (top row) we see that the spatial structure of the steady and fundamental is very similar.Furthermore, the steady and the first unsteady component are essentially in phase with each other indicating that the superposition of the two fields will be comparable to a sinusoidal modulation of the steady component in time, indicative of the flow approaching the quasi-steady limit.At the other extreme, Wo = 60 (bottom row), the first unsteady component has a plug-flow type spatial uniformity with a thin shear layer close to the wall whereas the higher harmonics have a low amplitude and are strongly localised close to the of wall on the inner side of the pipe.At this parameter combination, the temporal variation of the baseflow is comparable to a simple sinusoidal shift of the steady profile.At intermediate Womersley numbers (Wo = 25, centre row), the flow structure of the first unsteady Fourier component is considerably more complicated with a large degree of spatial non-uniformity and the emergence of strong secondary vortices shifted toward the sides of the pipe.Compared to the other plots, the amplitudes of the higher harmonics are orders of magnitude larger although all cases are computed at the same value of  which showcases the increased non-linear interaction between the components.The overall structure and variation of the baseflow for  = 0.3 is similar, albeit exhibiting stronger local velocity gradients.

Physical Floquet spectrum and temporal convergence
Reconsidering the definition of the Floquet eigenfunctions from equation (2.12), we observe that for a purely imaginary shift of the eigenvalue  by an integer multiple of the base frequency Ω ( ∈ N), we have (5.1) The expansions in the square brackets are equivalent since Fourier series are not affected by finite shifts in the coefficients.This means that the eigenfunctions of F   are invariant with respect to a purely complex shift by an integer multiple of the base frequency Ω and its spectrum is periodic along the imaginary axis with the period Ω.This means that the physically relevant spectrum is contained in a band of width Ω along the imaginary axis that contains the full information about the system.
When solving the Floquet eigensystem numerically, the equivalence of the shifted Fourier series is lost in general due to the numerical truncation at finite .This also destroys the exact periodicity of the spectrum since the fringe Fourier components of the problem lack the proper coupling with the adjacent modes.Figure 6a shows the most unstable part of the Floquet spectrum (circles) for a representative case computed with   = 11 together with the spectrum of the corresponding steady case (red crosses).We can clearly see the spurious eigenvalues appearing at the fringes that are erroneous in both frequency and growth rate.The existence of these spurious modes precludes the use of the "largest real part"option in the iterative eigenvalue solver instead of the shift-and-invert strategy to find the dominant eigenvalue.In order to assess the convergence of a given eigenpair, we can compute the kinetic energy  ( ) , −  ⩽  ⩽   of the Fourier modes of the eigenvector using equation (2.17) where the numerical evaluation of the integral over the Chebyshev points in the radial direction involves the Clenshaw-Curtis quadrature.From Fourier theory we know that the energy of the modes decay exponentially to zero as the wavenumber increases since they are smooth functions of time.In this work, the criterion for convergence of the Floquet eigenpairs  is max ⩽ 10 −6 . (5. 3) The eigenstates corresponding to the central sub-problem, i.e. whose Fourier spectrum is centred in the resolved frequency range, are typically better resolved in Fourier space due to the shift-equivalence of equation (5.1).The eigenstates that are converged according to (5.3) are marked with full symbols in figure 6a and the corresponding physical spectrum is shown in figure 6b.For each converged eigenpair, the kinetic energy spectrum is shown in figure 6c showing that the spectra are not symmetric since the eigenfunctions are complex Fourier series.The kinetic energies of the Fourier components also show the expected exponential decay towards to higher harmonics.In the computation of the Floquet spectra at consecutive pulsation amplitudes, the previous eigenvalue is used as a complex shift in order to track the variation of the most unstable eigenvalue.

Variation of the leading Floquet exponent with the pulsation amplitude
For small values of the pulsation amplitude , the baseflows can be computed accurately for all considered Womersley numbers.Close to the steady limit, the cycle-to-cycle growth rate has a quadratic dependence on , which is reminiscent of the studies of pulsating Poiseuille flow where this relation has been proved analytically (see e.g.Hall (1975); Von Kerczek (1982)).In this work, we extract the initial variation from the Floquet eigenvalues for 0.001 ⩽  ⩽ 0.003 using a least squares fit.To cross-check the quadratic fits, the growth rates are extrapolated to the steady case ( = 0) and compared to the exact value Re( 0 ).The error in the growth rate was found to be below 10 −5 for all considered cases (< 1% Re( 0 )).The results for the initial curvature of the temporal growth rate as a function of  are shown in figure 7 for each  over the considered Womersley number range.For the frequency range with stronger variation, the range of Womersley numbers was traversed with a higher resolution of ΔWo = 1.We observe that the initial variation of the real part of the Floquet eigenvalue is not monotonous across the considered frequency range.Although the initial effect of pulsation (a) (b) Figure 7: Variation of the initial curvature ( → 0) with respect to the pulsating amplitude of the real temporal growth rate computed using least squares from data for  ⩽ 0.003.Figure 7a: Overview and full range of data.Figure 7b: Close-up of the values close to zero (box in figure 7a).
is overwhelmingly stabilising (i.e. 2 Re()/ 2 →0 < 0), there exist frequency ranges for both curvatures in which the pulsations lead to a slightly increased instability.For  = 0.3 only very low Womersley numbers (Wo ⩽ 12) are destabilising.For lower curvatures, these frequencies are stabilising and we see instead two instability-promoting frequency bands, 13 ⩽ Wo ⩽ 18 and 32 ⩽ Wo ⩽ 45.In spite of the existence of some destabilising frequency bands for both curvatures, the dominant effect of pulsations at intermediate frequencies is to stabilise the flow.The highest stabilisation rates are 1-3 orders of magnitude stronger, for  = 0.1 and  = 0.3 respectively, than the corresponding maximum destabilisation.For  = 0.1, the most stabilising frequencies are 24 ⩽ Wo ⩽ 30 with a maximum at Wo = 25 and the width of this frequency band increases to 13 ⩽ Wo ⩽ 75 for  = 0.3 where the most stabilising frequency is Wo = 38.Assuming the stability characteristics follow the initial quadratic dependence for higher values of , these stabilisation rates suggest that pulsation amplitudes of  = 0.022 and  = 0.005 are sufficient to render the pulsating flow linearly stable at  = 0.1 and  = 0.3, respectively.The departure of the growth rate variation with  from the initial quadratic trends is considered in detail below.
Starting from the steady spectrum, the variation of the leading Floquet eigenvalue is tracked as  increases for all computed baseflows.As the pulsation amplitude increases, the Fourier spectrum of the Floquet eigenfunctions broadens.This general feature of Floquet problems has also been noted in the case of pulsating Poiseuille flow (Pier & Schmid 2017), who also noticed that this effect is exacerbated at lower Womersley numbers.Also in the present study, the number of Fourier modes that need to be retained in order to converge the leading Floquet eigenpair increases particularly sharply for the lowest Womersley number, as shown in figure 8.The data shows that, while the eigenvalue computations at high Womersley number required at most   = 7 Fourier modes even for  = 0.5, the spectra atWo = 5 could not be converged past  = 0.01 even considering the largest problem size for each curvature.An interesting detail to note is that the spectral width of the eigenfunctions, on the one hand, grows faster for lower curvatures, and, on the other hand, does not increase monotonically with decreasing Womersley number.Indeed, the pulsation frequencies at which the baseflow problems are the stiffest and the temporal growth rates are most damped (in particular for (,Wo) = (0.1, 25) and (0.3, 40), respectively), also the temporal complexity of the linear  solutions increases the faster.This being said, the cases at the lowest considered Womersley numberWo = 5 required the largest   to converge for a given pulsation amplitude.
For the cases where the dominant eigenpair could be converged, the variation of the temporal growth rate as a function of Wo and  are shown in figure 9 for  = 0.1 and  = 0.3.In each figure, the left plot shows the computed real temporal growth rates in the full parameter space whereas the right plots zoom in on the region of low pulsation amplitudes and compare the computed values to the initial quadratic dependence extrapolated to higher values of  (lines).The data confirms the results from the analysis for  → 0 showing that the pulsations have little effect on the intercyclic growth rate for both very low and very high Womersley numbers.Although for  = 0.3 the cases at Wo = 100 are more stabilised than the corresponding case at  = 0.1, high pulsation amplitudes are necessary to make the effect noticeable.For all but the two lowest considered Womersley numbers the current analysis suggests that the pulsating flow at (Re, ) = (1700, 0.3) is linearly stable for  between ≈ 0.006 (Wo = 40) and  = 0.25 (Wo = 100).For the lower curvature case, the data for higher pulsation amplitudes shows that the frequency range Wo ⩾ 45 is also generally stabilising but it seems that increasing the pulsating amplitude will lead to linearly stability only in a very thin frequency band around 20 ⩽ Wo ⩽ 25.Indeed, already for Wo = 30, the initial downward trend of the growth rate stops at  ≈ 0.08.It is likely that the cases at higher curvature will deviate in a similar way from their initial trends, especially in the frequency range Wo = 20 − 45, but the present data cannot corroborate this because the Newton-Raphson iteration could not converge the baseflow past  = 0.01 in this parameter range.
This conjecture is supported by measuring how quickly the available data diverges from the initial quadratic variation for different Womersley numbers.To compare the trend to the computed data, we define   as the quadratic extrapolation of the growth rate using the data from figure 7ato higher pulsation amplitudes.The magnitude of the difference in the growth rates |Re( −   )| is shown for both curvatures in figure 10, relative to the initial growth  rate Re( 0 ) at  = 0.The data shows that the deviation scales uniformly with  4 in the entire parameter range.Only the data for very high Womersley numbers and low pulsation amplitudes does not seem to follow this trend, but this could be due to more noisy data because of the tolerance on the eigenvalue computation and the small initial stabilisation rates compared to the steady case in this parameter range.Unsurprisingly, the cases that are more sensitive to the pulsations also deviate from the initial trend earlier.

intracyclic variation of the linear perturbation
The Floquet eigenvalue and the corresponding multiplier describe the intercyclic evolution of the linear perturbation, i.e. relates two states that are one period apart, integrating the instantaneous variations.In order to understand the intracyclic variations, we need to consider the associated Floquet eigenfunction which encapsulates all the temporal information within the period.

Instantaneous growth rate and perturbation energy
Reconstructing the perturbations u the Floquet using (2.12), the intracyclic variation of the instantaneous growth rate () of the perturbation is computed using (2.19) where W encapsulates both the spatial integration weights stemming from the Chebyshev-Fourier discretisation and the action of the restriction operator, L = L 1 + L 2 () is the full linear operator constructed using the instantaneous baseflow  and it is understood that all quantities, with the exception of W , are time-dependent and evaluated at the same instant .
The instantaneous growth rate itself gives an indication of the intracyclic behaviour of the perturbation but in order to have a reasonable comparison between cases, we need to consider the aggregate energy growth over the pulsation period since the energy growth is determined by both the instantaneous growth rate and the time span over which it acts.The practical computation of the perturbation energy needs to be approximated since the integrand in (2.20) is time-dependent.In practice we approximate the growth curve () as a piece-wise linear constant over the pulsation period in order to compute the perturbation energy as an ordered product of exponentials (5.5)where   = /Δ,   = Δ and   = (( +1 ) + (  ))/2.Since the temporal complexity () is given by the Fourier of the baseflow and the Floquet eigenfunctions, the instantaneous growth rate be obtained by explicitly computing at a number   of discrete points the period where  ⩾   , where the limit   is given by Nyquist-Shannon theorem based on the highest retained Fourier mode.The evaluation (5.5) via the ordered product of exponentials is very sensitive to the chosen time interval Δ.Therefore, a high temporal resolution is required for () which can be obtained via interpolation of the discrete data given that all relevant frequencies are resolved.In this work, we choose Δ ≈ 10 −3 so that   also depends on the period .
The instantaneous growth rate of the dominant Floquet eigenfunctions over the period was computed for each curvature for a range of Womersley numbers where the effect of pulsation on the temporal growth rate was largest.Based on these growth rates, the evolution of the perturbation energy can be traced over time and compared to the growth of the corresponding steady case.The perturbation growth rates and energies are shown in figure 11 for the highest pulsation amplitude for which the Floquet problems could be converged in the entire considered Womersley number range.Comparing the instantaneous growth curves with the sketch of the pulsating pressure gradient (thick black line, not to scale) we see the typical feature of pulsatile flows that are stabilised ( < 0) when the baseflow is accelerated and destabilised ( > 0) when it is decelerated.The shifts with respect to the baseflow pulsations are related to the phase lag that exists between the pressure and the oscillating boundary layer on which the instabilities develop.Furthermore, we observe that even at low pulsation amplitudes of  = 0.05 and  = 0.01 for  = 0.1 and  = 0.3, respectively, the maximal instantaneous growth rates are roughly an order of magnitude larger than the growth rate of the unstable mode of the corresponding steady problem (dashed black line).Interestingly, even though some of the instantaneous growth rates over the pulsation period exhibit similar excursions, the resulting energy growth curves vary considerably.This is due partly to the details of the intracyclic variations but mainly to the pulsation period  ∝ Re/Wo 2 that quickly increases with decreasing Womersley number leading to longer intervals of growth and decay.The instantaneous growth curves, in particular for  = 0.1,Wo = 20 − 25, also clearly show the appearance of higher harmonics in the perturbation that lead to a deviation of the curves from the quasi-sinusoidal variation (e.g. at  = 0.3,Wo = 20).
In order to put the data shown for a single pulsation amplitude  for each curvature into context, we compute the amplitude of Δ =  − 2 0 , the variation of the intracyclic energy growth rate compared to the steady case.The last row of figure 11 shows how the extremal intracyclic growth/decay rates vary as a function of pulsation amplitude and frequency.Indeed, the extremal growth rates vary much more strongly with  in case of higher curvatures.The data for higher pulsation amplitudes available only for Wo = 45 indicates the reason for the immense stabilisation: the maximal decay rate increases linearly with  while the maximum growth rates a plateau already at  = 0.03.
In the parameter ranges considered in this work, the intracyclic energy growth is likely not large enough to trigger transition.The extreme intracyclic perturbation growth documented for low pulsation frequencies in Pier & Schmid (2017) (see figure 3 in the reference) are not found in the present case.Nevertheless, given the intracyclic growth rate variations found in the present study for small values of  and their increasing trend with , it is possible that considerably larger intracyclic growth rates are possible also in the pulsating torus that might lead, when non-linearities are considered, to a "ballistic"regime similar to the one described in pulsating channel flow (Pier & Schmid 2017).It should be noted in this context that a large intracyclic energy growth is a product of both large instantaneous growth rates and long periods they can act over.Considering that the period  ∝ Re/Wo 2 , the considered Reynolds number in this work (Re = 1700) leads to considerably shorter periods and thus lower integrated energy growth than the Re = 7500 considered in Pier & Schmid (2017) assuming the same Womersley number and intracyclic growth rate variation.The crosses correspond to the symmetric (s) and antisymmetric (a) part of spectrum Λ 0 of the corresponding steady problem for reference.The dominant eigenvalue moves long the orbit in anti-clockwise direction during the period.

eigenvalue orbits in the spectrum of the local operator
In a recent study, Kern et al. (2022) documented the formation of subharmonic eigenvalue orbits in the periodic spectrum of pulsating plane Poiseuille flow as the pulsation amplitude is increased which appear due to the existence of spectral degeneracies called exceptional points at which two or more eigenpairs coalesce.In the reference, it could be shown that for low pulsation frequencies the instantaneous growth rate of the asymptotic linear perturbation follows that of the leading local eigenvalue branch closely.If that particular eigenvalue orbit merges with another at an exceptional point to form a subharmonic orbit, a branch transition becomes necessary since the linear perturbation is aligned with the dominant Floquet eigendirection and thus has the same periodicity as the baseflow.The branch transition processes were shown to involve strong non-normal growth bursts via the Orr mechanism.
Since the formation of subharmonic eigenvalue orbits is related to the non-normality of the linearised Navier-Stokes operator, they were predicted to appear in any time-periodic system involving a non-normal linear operator and thus also in the present case.
To confirm this prediction, the leading eigenpair { 1 (, ),  1 (, )} of the local problem (i.e. the local in time linear operator constructed from a snapshot of the time-dependent pulsatile baseflow) was tracked over the pulsation cycle for increasing values of the pulsation Since the eigenvalue problems are too large for many eigenvalues to be computed, the eigenpair was tracked over the period in small steps of Δ  = /100 using standard Rayleigh-quotient iteration initialised with the eigenvector 1 from the previously computed eigenpair which converges quickly to the next eigenpair.The smoothness of the spectrum as  is varied guarantees that, given the step size Δ  is small enough, the eigenpair can be tracked (Kern et al. 2022).In order to always follow the leading eigenvalue for different values of , the value of  1 (, ) was first computed  =  and varying  in small steps starting from the steady spectrum using the same iterative technique.The orbit was then traversed backwards since we aim to follow the mode during the maximal local amplification phase located at the end of the period due to the phase shift.In case of a suhbarmonic orbit, going forward in time would instead follow the orbit into the damped part of the spectrum.
The variation of the leading local eigenvalue orbit over one period for Re = 1700,  = 0.1, Wo = 35 is shown in figure 12 where the pulsation amplitude is varied in the range  ∈ [0, 0.103].We see that the eigenvalue traces form closed orbits for  ⩽ 0.1 and then, starting from  = 0.103, the eigenvalue does not return to its initial locus at the end of the period indicating that the orbit has merged with at least one other orbit.Just before the subharmonic orbit forms, at  = 0.1, we can clearly see the increased sensitivity of the damped part of the eigenvalue orbit, which is a hallmark of the proximity to an exceptional point in the spectrum.No attempt was made to trace the full resulting subharmonic eigenvalue orbit due to the fact that it involves very damped eigenvalues which cannot be computed with high accuracy with the present spatial resolution (see A).In spite of the existence of subharmonic eigenvalue orbits, no clear evidence of non-normal growth bursts comparable to those documented for pulsating plane channel flow was found.This is likely due to the restriction of the number of retained Fourier modes   as well as the fact that the could not be converged at high values of .The tell-tale in the growth rate curves appear at a pulsation amplitude of around 10% and a Fourier analysis of the example case analysed in detail in figure 11 of the reference (Re = 7500, Wo = 18, Q = 0.18) indicates that, using the present method, a temporal resolution of   ≈ 33 would be required for convergence of the dominant eigenpair (note that in the the data was obtained using a time-stepper approach).Furthermore, also here the lower Reynolds number suggests that the occurrence of the described phenomena would be shifted to even higher pulsation amplitudes.Nevertheless, based on the fact that subharmonic eigenvalue orbits exist and the Floquet solutions exhibit a similar behaviour to the dominant linear perturbations computed for pulsating plane channel flow in the literature (Pier & Schmid 2017;Kern et al. 2022) regarding the intracyclic growth rate variation and the widening of the Fourier spectrum for low pulsation frequencies, it is to be expected that similar non-normal growth bursts will also appear in the present case, albeit for parameter combinations inaccessible to this study.

Spatial structure of the Floquet eigenfunctions
The data concerning the instantaneous growth rates as well as the perturbation energy growth should be considered in combination with the spatial structure of the eigenfunctions.We choose a representative case for each curvature at the pulsation frequency that is most stabilising and at the highest pulsation amplitude for which converged results are available, i.e.Wo = 25,  = 0.05 and Wo = 25,  = 0.05, for  = 0.1 and  = 0.3, respectively.For these parameter combinations, the evolution of the dominant Floquet eigenfunction over the period for each curvature is shown in figure 13 together with the corresponding instantaneous growth rates.Given the low pulsation amplitudes, the structure of the eigenfunctions is very similar to the corresponding steady case (figure 2).
For  = 0.1, the dominant mode is symmetric and spread over the pipe cross-section.The peak decay rates are associated with a more localised mode on the Dean vortices close to the inner wall with strong in-plane velocity perturbations.As the baseflow accelerates and in particular the Dean vortices strengthen the perturbation spreads out.With the deceleration of the flow, the mode localises on the Dean vortices again.For  = 0.3, the eigenmode is antisymmetric and there is barely any spatial variation over the pulsation cycle.Instead, as the growth rate of the perturbation increases, the mode focuses even more on the part of the Dean vortices close to the inner wall and the in-plane variation is reduced.During the baseflow acceleration phase, this process is reversed.Also in this case, decay is associated with increased in-plane velocities in the perturbation.
Considering the comparatively small spatial variations of the eigenmodes over the pulsation cycle, the instantaneous growth rates are surprisingly high, in particular for  = 0.3 where the considered perturbation corresponds to a pulsation amplitude of only  = 0.01, five times less than for the perturbation shown for  = 0.1.This sensitivity is likely due to the fact that the instabilities feed off of the high local shear of the Dean vortices close the inner pipe wall that undergo the strongest variations even for small pulsation amplitudes.Especially at higher curvatures, where the dominant perturbation is extremely localised, its growth/decay rates are particularly sensitive to minute changes in the baseflow.The fact that the dominant instability mode for  = 0.3 is much more localised in the cross-section is probably the main factor explaining the extreme stabilisation rates we observe for this configuration, even for very small pulsation amplitudes.

Effect of baseflow temporal fidelity on the Floquet exponent
The size of eigenvalue problems encountered in this study are heavily dependent on the spatial and temporal resolution not only of the perturbation eigenfunctions but also of the baseflow itself.Given the low amplitudes of the higher harmonics of the baseflow for many of the parameter combinations, it is reasonable to consider whether a reduced order model of the baseflow can be used for the stability calculations which would lead to a strong reduction of the problem size.In order for the problem size effectively, we need to reduce the width of the baseflow Fourier spectrum.It is not clear a priori how this will affect the Floquet computations since this step not only means a reduction on the fidelity of the baseflow data but also a structural change within the Floquet problem via the reduction of the coupling width between adjacent modes.
In this work we consider two crude methods to reduce the problem size.First, we consider the baseflow solution computed with   = 3 and then simply discard the higher harmonics ( ⩾ 2) for the Floquet computations.The number of block sub-diagonals in the Floquet problem is directly proportional to the retained baseflow components and the resulting reduced-order Floquet problem has roughly 25% fewer non-zero elements.The second option considered in this work is to avoid the comparatively expensive accurate baseflow computation with many Fourier modes entirely and to use the baseflow computed with   = 1 directly in the Floquet problem.We stress here that the first option introduces no error in the baseflow computation but reduces the coupling withing the Floquet problem via the truncation of the baseflow expansion whereas the second option introduces errors both via the reduced coupling but also due to neglect of higher harmonics in the computation of the baseflow itself.Especially the second option, although computationally efficient, cannot be expected to yield very accurate results but may be a tool to get a quick overview of the parameter space to choose the regions for which to refine the data.
In order to assess the error introduced by the reduction of the baseflow expansion, we track the dominant Floquet eigenvalue in the parameter range and compare the resulting temporal growth rates with the data from figure 9.The results of the comparison are shown in figure 14 where the baseflow truncation is denoted with the subscript  = 1 (full lines) whereas the computations using the baseflow computed at   = 1 are denoted with subscript  = 2 (dashed lines).The data is shown only for a subset of the considered Womersley numbers for clarity but the trends are similar throughout the parameter range.The data show that the error in the eigenvalues scales as  4 in both cases, the same as the deviation from the initial quadratic variation (albeit at a lower level).This indicates that the scaling is defined by the errors due to the reduced coupling, while the inaccuracies due to the baseflow computed with low temporal resolution affect only the proportionality factor.We see that for high frequencies (Wo > 60) and low pulsating amplitudes ( ⩽ 0.1 and  ⩽ 0.01 for  = 0.1 and  = 0.3, respectively) the error is low and is dominated by the tolerance of the eigenvalue solver (horizontal line).In this parameter range, the baseflow is essentially sinusoidal.For intermediate Womersley numbers, the error quickly increases by several orders of magnitude.While the error of the data computed with the truncated baseflow is usually 5-10 times smaller than for option 2, for high curvatures close to the peak stabilisation at 40 ⩽ Wo ⩽ 50, the difference between the two reduced-order models becomes negligible.

Conclusions
In this paper we analyse the linear stability of fully-developed pulsatile flow in toroidal pipes using Floquet theory.The two-dimensional streamwise invariant baseflows driven by a sinusoidal pressure gradient are computed for two curvatures,  = 0.1 and  = 0.3, at Re = 1700 for Womersley numbers and pulsating amplitudes in the range (Wo, ) ∈ [5, 100] × [0, 0.5] using a Newton-Raphson iteration in Fourier space applied to the full non-linear Navier-Stokes equations expressed in toroidal coordinates.Due to the strong nonlinearities in the equations for certain values of the Womersley number, the baseflow could not be converged for all pulsation amplitudes.At the chosen Reynolds number the steady flow at both curvatures is supercritical and the effect of pulsation on the linearly stability is analysed by tracking the leading Floquet eigenpair as the pulsation amplitude is increased.
The study shows that pulsations have an overwhelmingly stabilising effect at both curvatures.Only thin frequency bands are found for which the baseflow modulation is slightly destabilising.It was found that for the parameter combinations, for which the governing equations exhibited increased non-linear coupling leading to complicated baseflow dynamics, the effect of the temporal baseflow modulation on the Floquet exponent was strongest and is linked with an higher spatio-temporal complexity of the corresponding eigendirection.The stabilisation due to the pulsations is particularly salient for  = 0.3, which is attributed to the strong localisation of the unstable eigenmode of the steady problem on a region of the baseflow with strong gradients which therefore is more sensitive to intracyclic modulations of the baseflow.At  = 0.1, the unstable mode of the steady torus is much more spread over the cross-section and is thus less affected by the pulsations.For the same reason the instantaneous intracyclic growth rate excursions, which are found to be orders of magnitude larger than the steady growth rate, are particularly strong at high curvature.For very high and very low Womersley numbers the baseflow structure is much simpler and the effect of the pulsations is strongly reduced.In these parameter regions, reasonably accurate results can be obtained at considerably reduced cost for small  by neglecting the higher harmonics in the computation of the Floquet exponents.
Comparing the present results with studies in literature of the stability of pulsating Poiseuille flow, the general trends are similar but the effects are stronger in the toroidal pipe.
The increased sensitivity is likely due to the spatial localisation of the modes as well as the stronger coupling of adjacent perturbation Fourier modes via the baseflow which is broadband in frequency via the non-linearity of the Navier-Stokes equation in toroidal coordinates in contrast to pulsating Poiseuille flow where the governing equations become linear.While the existence of subharmonic eigenvalue orbits predicted in Kern et al. (2022) was confirmed, the associated transient growth bursts during the branch transition processes were not clearly seen.Given that these phenomena appear unambiguously in pulsating channel flow only at pulsation amplitudes far beyond the formation of subharmonic eigenvalue orbits, it can be expected that the Floquet eigendirections will exhibit similar dynamics also in toroidal pipes at higher values of  outside the scope of this study.Further studies will need to consider higher pulsation amplitudes, in particular for intermediate Womersley numbers, in order to corroborate this conjecture.At the same time, non-linear simulations could shed light on the question whether the intracyclic perturbation energy growth, which is seen to increase sharply with , will eventually lead to a ballistic regime similar to that documented for pulsating channel flow (Pier & Schmid 2017).

Acknowledgement
The authors thank Ardeshir Hanifi for helpful discussions on the implementation of the Floquet solver.

Funding
Financial support for this work was provided by the European Research Council under grant agreement 694452-TRANSEP-ERC-2015-AdG and the Swedish Research Council (VR) through the grant no.2017-04421.The computations were performed on resources by the Swedish National Infrastructure for Computing (SNIC) at the PDC Center for High Performance Computing at the Royal Institute of Technology (KTH).

Appendix A. Code validation and spatial convergence study
In order to validate the linear operators for the flow in the torus, we compute the critical Reynolds number for the considered radii of curvature  and compare the results with the data from Canton et al. (2016) computed using the finite element code PaStA based on an unstructured grid.For each case, the radial and azimuthal resolution is the same as in the other computations in this work.Due to differing length scales used for normalisation, the Reynolds number in the reference (Re  ) differs from the one in this work (Re) with the relation Re  = 2Re and the streamwise wavenumber  relates to the variable  in the reference as  = .The Floquet operator is constructed using the same linear operators in its sub-blocks.The results of the comparison are shown in Tab. 1 showing a very good agreement between the codes (relative error below 0.2%), the present method systematically yielding a slightly lower critical Reynolds number.
With the discretisation and linear operators validated against literature, we need to assess the required spatial discretisation.Preliminary studies have shown that the spectrum of the linear operator was particularly sensitive to the azimuthal resolution (  ) once an adequate radial resolution of   = 30 was chosen which is confirmed via the comparison of the growth rate of the leading eigenvalue shown in Tab. 1. Figure 15 shows the variation of the spectra for the steady case for different choices of   for the two considered curvatures.For  = 0.1, the spectra for   = 110 (squares) and   = 130 (crosses) are indistinguishable such that a resolution of   = 110 was chosen for the cases at low curvature.For  = 0.3 on the other hand, we observe that even an azimuthal resolution of   = 150 is not sufficient to accurately capture the portion of the spectrum shown here.Since the leading eigenvalue, which is the main focus of this study, barely changes between   = 130 and   = 150, the latter was considered an adequate resolution for the present study.

Appendix B. Discretised operators for the linearised Navier-Stokes equations in toroidal coordinates
This section presents the operators appearing in the Navier-Stokes equations linearised around a baseflow U and expressed in toroidal coordinates and discretised using the combined Chebyshev-Fourier collocation method.The discretised velocity components used in the definition of L 2 are denoted U  , U  and U  , where the pre-multiplication with 1  1 to obtain diagonal matrices is omitted for notational clarity.The baseflow pressure does not enter in the computation of the linear operators.We denote D ()  , D ()  , and D ()  , the  th order two-dimensional azimuthal and radial differentiation matrices.Note that the radial differentiation matrices for the streamwise D ()   , and radial/azimuthal D ()  , components are not identical.The two-dimensional differentiation matrices are computed using the standard one-dimensional radial (Chebyshev) and azimuthal (Fourier) differentiation matrices from DMSUITE by Weideman & Reddy (2000), combined based on the chosen data structure.The streamwise differentiation matrices are given in terms of the chosen streamwise wavenumber  and the curvature  as The operator L 1 representing the diffusive/pressure terms has the following block structure With the definition of ℎ  =   sin() + 1 , (B 3) the sub-blocks are given by: (i) -component of the momentum equation: and with and where  =    2 and  =  − 1,  ∈ {−1, 1} denotes the chosen symmetry, antisymmetric or symmetric, respectively, and J  is the exchange matrix of order  with ones on the antidiagonal.Empty parts of the matrices are zeros.In the case even    , similar but more convoluted expressions for T and T  can be derived.

Figure 1 :
Figure 1: Computational meshes for the pipe cross-section for each considered curvature. is the number of radial Chebyshev gridpoints and   is the number of equispaced azimuthal gridpoints.Due to the assumed symmetry, the computations are performed on half the grid.

Figure 2 :
Figure 2: Steady baseflows, dominant eigenvector and spectra at Re = 1700 at  = 0.1 (figures 2a-2c) and  = 0.3 (figures 2d-2f), with the real part of the eigenvector corresponding to the most unstable eigenvalue for  = 1.7 and  = 2.1, respectively.The left half of each flow field shows   and the right half shows the in-plane velocity   .Left: Steady baseflow.The left half of the figure is overlaid with contours of the streamfunction.Center: Leading eigenvector.Right: Spectrum of the linear operator with symmetric (sym) and antisymmetric (asym) modes, red and blue, respectively.The line indicates the smooth variation of the dominant eigenvalues for variations of  in the range of roughly ±20%.The circles are realisable values (/ ∈ Z).

Figure 3 :Figure 4 :
Figure 3: Parameter map for the pulsating baseflow for  = 0.1 (left) and  = 0.3 (right), both computed with three pulsating flow components (  = 3).The circles correspond to computed cases with filled (open) symbols indicating whether the baseflow solution is considered accurate (approximate) following equation (4.2).The red crosses indicate where the Newton solver diverged (no attempts at higher  were made).

Figure 5 :Figure 6 :
Figure 5: Examples of the baseflow profiles at Re = 1700. = 0.1 and  = 0.01 including   = 3 unsteady Fourier components with, from left to right in order, the steady component ( = 0) and the real (top) and imaginary (bottom) parts for the unsteady components  ∈ [1, 2, 3].The left half of each figure shows the   and the right half shows the in-plane velocity    .For the steady component the left half of the figure is overlaid with contours of the streamfunction.Figures 5a-5d:Wo = 5.Figures 5e-5h:Wo = 25.Figures 5i-5l:Wo = 60.Note the different colourbars for each component.
Figure 6: The least stable part of the Floquet spectrum for symmetric modes at Re = 1700,  = 0.1,Wo = 25,  = 0.02 and  = 1.7 with   = 11 Fourier components (circles) compared to the corresponding symmetric steady spectrum (red crosses).Figure 6a: Floquet spectrum (circles) with the converged eigenvalues marked with filled symbols.Figure 6b: Physically relevant spectrum showing only converged eigenvalues, coloured by their growth rate.Figure 6c: Fourier spectrum of the kinetic energy of each converged eigenpair normalised by their respective sum.The colour code is the same as in figure 6b and the spectrum of the dominant eigenvector is emphasised.The dashed line indicates the chosen convergence threshold.

Figure 8 :
Figure 8: Required number of retained Fourier modes   for convergence of the leading Floquet eigenvalue as a function of  andWo.The area shaded in red indicates where the baseflow was not computed and the red crosses indicate cases that could not be converged with the respective maximum   .

Figure 9 :
Figure 9: Variation of the growth rate of the leading Floquet eigenvalue as a function of Wo and  at Re = 1700 for  = 0.1 (top row) and  = 0.3 (bottom row).In the frequency range where strong variations occur, the amplitude range is traversed in smaller steps.A variation in the line weight is used for improved readability.Left: Real part of the Floquet eigenvalues for  ⩾ 0.005 compared to steady growth rate ( 0 , red cross).The full (dashed) lines and filled (open) symbols indicate that the baseflow computed at   = 3 is accurate (approximate).Right: A closeup of the area in the black box.The line indicates the extrapolation   .

Figure 10 :
Figure10: Difference between the exact Floquet eigenvalue  and the quadratic extrapolation   relative to steady eigenvalue (only cases with converged baseflow shown).Note that only  ⩾ 0.005 is shown since the data for smaller  was used to compute the quadratic dependence.

Figure 11 :
Figure 11: Intracyclic variation of the dominant Floquet eigenfunction compared to the unstable mode for  = 0 for  = 0.1 (left) and  = 0.3 (right) for intermediate Womersley numbers.The colourmap is the same as in previous figures.Top row: Instantaneous growth rate Re() over the pulsation period .Symbols = computed, lines = interpolated.Black line: Baseflow pulsation (not to scale).Centre row: Integrated perturbation energy over time.Bottom row: Variation of the absolute value of the extremal instantaneous intracyclic growth (decay) rates of Re( − 2 0 ) shown with full (dashed) lines as a function of .The vertical dashed red line indicates the parameter combinations for the data in figures 11a-11d.

Figure 12 :
Figure 12: Evolution of the eigenvalue orbit of dominant  1 (, ) as a function of  ∈ [0, ]  for Re = 1700,  = 0.1,Wo = 35.The black line indicates the variation of the leading eigenvalue with  at  =  from which the tracking was initialised.The crosses correspond to the symmetric (s) and antisymmetric (a) part of spectrum Λ 0 of the corresponding steady problem for reference.The dominant eigenvalue moves long the orbit in anti-clockwise direction during the period.

Figure 14 :
Figure14: Difference in the growth rates computed using a reduced order model   ( = 1: full lines,  = 2: dashed lines) compared to the accurate data .For  1 the Floquet problem was solved retaining only one unsteady baseflow Fourier component from the solution computed at   = 3.For  2 the approximate baseflow solution computed at   = 1 was used.The thin dashed horizontal line indicates the threshold for the eigenvalue computation.

Figure 15 :
Figure 15: Spectra for the steady flow in the torus at different azimuthal resolutions.The red (blue) symbols correspond to symmetric (antisymmetric) modes.

Table 1 :
Critical Reynolds numbers for the onset of linear instability for the considered curvatures from literature compared to the present work.