Nils Tilton

Denis Martinand

Taylor-Couette-Poiseuille flow with a weakly permeable inner cylinder: absolute instabilities and selection of global modes

Keywords: absolute/convective instability, nonlinear instability, Taylor-Couette flow

Variations in the local stability of the flow in a Taylor-Couette cell can be imposed by adding an axial Poiseuille flow and a radial flow associated with one or both of the cylinders being permeable. At a given rotation rate of the inner cylinder, this results in adjacent regions of the flow that can be simultaneously stable, convectively unstable, and absolutely unstable, making this system fit for studying global modes of instability. To this end, building on the existing stability analysis in absolute modes developing over axially invariant base flows, we consider the case of axially varying base flows in systems for which the outer cylinder is impermeable, and the inner cylinder is a weakly permeable membrane through which the radial flow is governed by Darcy's law. The frameworks of linear and nonlinear global modes are used to describe the instabilities and assess the results of direct numerical simulations using a dedicated pseudospectral method. Three different axially evolving set-ups are considered. In the first, fluid injection occurs along the full inner cylinder. In the second, fluid extraction occurs along the full inner cylinder. Besides its fundamental interest, this set-up is relevant to filtration devices. In the third, fluid flux through the inner cylinder evolves from extraction to injection as cross-flow reversal occurs. In agreement with the global mode analyses, the numerical simulations develop centrifugal instabilities above the predicted critical rotation rates and downstream of the predicted axial locations. The global mode analyses do not fully explain, however, that the instabilities observed in the numerical simulations take the form of axial stacks of wavepackets characterized by jumps of the temporal frequency.

Introduction

Linear stability analysis has matured substantially since its first successful predictions of thermoconvective instabilities in Rayleigh-Bénard convection by [START_REF] Rayleigh | On convection currents in a horizontal layer of fluid when the higher temperature is on the under side[END_REF] and centrifugal instabilities in Taylor-Couette flow by [START_REF] Taylor | Stability of a viscous liquid contained between two rotating cylinders[END_REF]. Nevertheless, these two classical systems remain actively studied to advance our understanding of hydrodynamic instabilities and test new approaches of stability † Email address for correspondence: denis.martinand@univ-amu.fr analysis [START_REF] Koshmieder | Bénard Cells and Taylor Vortices[END_REF], provides a thorough presentation). Though academic, Rayleigh-Bénard and Taylor-Couette flows bear sufficient resemblance to practical systems, and they can be readily studied experimentally to verify and motivate analytical approaches. The current work demonstrates that Taylor-Couette-Poiseuille (denoted TCP hereinafter) flow with a permeable inner cylinder provides a new opportunity to investigate how instabilities develop and tend to synchronize themselves in spatially varying open flows, a class of flows by which numerous industrial and geophysical systems are described. The current work, for example, is motivated by applications of supercritical TCP flow to membrane separation [START_REF] Hallström | Decription of a rotating ultrafiltration module[END_REF][START_REF] Beaudoin | Plasma filtration in Couette flow membrane devices[END_REF]. More broadly, however, providing rigorous theoretical frameworks for synchronized, global, instabilities is pivotal to understand and predict how spatially varying open flows lose their stability and transition to turbulence. By comparing the results of theoretical analyses to the results of dedicated direct numerical simulations (DNS) of TCP flows, the current work shows that the dynamics of these global instabilities are still not fully understood beyond their critical conditions.

Early linear stability analyses of Taylor-Couette flow with a stationary outer cylinder considered the ideal case of infinitely extended concentric cylinders. Temporal stability analyses found standing wave instabilities that vary sinusoidally in the axial direction. Excellent agreement between experimental, numerical, and theoretical results are observed for the transition to supercritical flows, owing to the pitchfork nature of the bifurcation. The addition of a mean axial flow, producing TCP flow, can similarly be studied theoretically by considering infinitely extended cylinders between which instabilities vary sinusoidally in the axial direction. This shows that a Hopf bifurcation produces instability in the form of propagating waves. Mean axial flow also stabilizes the appearance of centrifugal instability in comparison to Taylor-Couette systems (see [START_REF] Diprima | The stability of a viscous fluid between rotating cylinders with an axial flow[END_REF][START_REF] Donnelly | Experiments on the stability of spiral flow between rotating cylinders[END_REF][START_REF] Chung | Hydrodynamic instability of viscous flow between rotating coaxial cylinders with fully developed axial flow[END_REF][START_REF] Johnson | Hydrodynamic stability of flow between rotating porous cylinders with radial and axial flow[END_REF], among many publications). For small mean axial velocity, instabilities appear as toroidal vortices travelling axially. With increasing mean axial velocity, these are superseded by travelling helical vortices. For sufficiently large axial flows, shear instabilities can also develop [START_REF] Meseguer | On the competition between centrifugal and shear instability in spiral Poiseuille flow[END_REF][START_REF] Cotrell | The connection between centrifugal instability and Tollmien-Schlichting-like instability for spiral Poiseuille flow[END_REF]. The current work, however, will consider only centrifugal instabilities.

The addition of mean axial flow motivates two new approaches to modelling instabilities, due to the fact the system is now 'open'. First, the mean flow-induced non-normality of the unstable modes can produce transient instabilities (see [START_REF] Chomaz | Global instabilities in spatially developing flows: non-normality and nonlinearity[END_REF], for a review). Second, considering the impulse response of the flow (the Green's function) produces instabilities of finite axial extent. This motivates a distinction between convective and absolute instabilities [START_REF] Huerre | Absolute and convective instabilities in free shear layers[END_REF], for instabilities in fluid flows). Convective instabilities grow in space and amplitude, but are advected downstream such that they eventually wash out of finite-length domains. In contrast, absolute instabilities grow sufficiently to overcome mean advection and invade an entire system, both up-and downstream. Convective/absolute stability analysis has been applied to Taylor vortices [START_REF] Büchel | Pattern selection in the absolutely unstable regime as a nonlinear eigenvalue problem: Taylor vortices in axial flow[END_REF] and extended to helical vortices [START_REF] Martinand | Absolute and convective instability of cylindrical Couette flow with axial and radial flows[END_REF][START_REF] Altmeyer | Islands of instability for growth of spiral vortices in the Taylor-Couette system with and without axial through flow[END_REF].

Most instabilities of practical interest develop on 'imperfect' base flows that vary spatially in more than one direction. In addition, end effects in open flows are more critical experimentally and numerically than in closed flows. In open flows, inlet conditions are advected by the mean flow and impact the whole domain, as observed by [START_REF] Babcock | Noise-sustained structure in Taylor-Couette flow with through-flow[END_REF][START_REF] Babcock | Noise amplification in open Taylor-Couette flow[END_REF] for TCP flows. Outlet conditions can also impact the instabilities developing within the system, as observed in rotating-disk flows [START_REF] Pier | Transition near the edge of a rotating disk[END_REF][START_REF] Appelquist | On the global nonlinear instability of the rotating-disk flow over a finite domain[END_REF]. Therefore, experiments and numerical simulations often use specific inlet and outlet configurations that modify physical or geometrical characteristics along the mean flow direction to isolate instabilities developing in the core of the system. For example, experimental set-ups of Rayleigh-Bénard-Poiseuille convection often use a tailored 'thermal entry' where a temperature difference between the plates is slowly imposed to a developed Poiseuille flow [START_REF] Nicolas | Bibliographical review on the Poiseuille-Rayleigh-Bénard flows: the mixed convection flows in horizontal rectangular ducts heated from below[END_REF]. Therefore, introducing a degree of inhomogeneity in the base flow of the stability analysis, in the form of variations of its physical or geometrical characteristics, is relevant, if not required.

When considering the stability of an axially varying base flow with autonomous boundary conditions, instabilities take the form of synchronized modes called 'global' modes. A classical example is the von Kármán vortex street [START_REF] Provansal | Bénard-von Kármán instability: transient and forced regimes[END_REF]. The direction of the mean flow is then an extra eigendirection and such global modes can be computed numerically (see [START_REF] Theofilis | Advances in global linear instability analysis of non-parallel and three-dimensional flows[END_REF], for a review). Computationally more affordable analyses of global modes are possible when the base flow varies slowly in the mean flow direction. Such base flows are often called 'weakly non-parallel'. In this case, analyses can consider a convective/absolute stability analysis that approximates the local base flow as axially invariant, albeit with the physical and geometrical parameters set at the local values of the axially varying base flow. For flows that are everywhere stable or convectively unstable, perturbations can grow while in convectively unstable regions, but they either decay when reaching a stable region or they advect out of the system. Global instability arises when self-sustained localized perturbations occur in a region of local absolute instability. Analytical approaches to global modes consequently build on how absolute instabilities travel and saturate over axially varying base flows and look for a specific axial location that can act as a wavemaker and govern the global mode by imposing its critical conditions and frequency.

Two classes of wavemakers have emerged (see [START_REF] Huerre | Local and global instabilities in spatially developing flows[END_REF][START_REF] Huerre | Open shear flow instabilities[END_REF][START_REF] Chomaz | Global instabilities in spatially developing flows: non-normality and nonlinearity[END_REF], for reviews). The first, based on the linear dynamics of the absolute instability, locates the wavemaker where this absolute instability reaches its maximum local amplification. If the local amplification decreases monotonically in the axial direction, global modes can be obtained in a semi-infinite domain and the wavemaker will occur at the inlet [START_REF] Monkewitz | Global linear stability analysis of weakly non-parallel shear flows[END_REF]. Conversely, the wavemaker will occur at the outlet if the local amplification increases monotonically [START_REF] Healey | Model for unstable global modes in rotating-disk boundary layer[END_REF]. If the local amplification increases then decreases axially, the wavemaker will settle within the domain [START_REF] Monkewitz | Global linear stability analysis of weakly non-parallel shear flows[END_REF]. Critical conditions and frequencies of global modes are computed from the linear dispersion relation for absolute instabilities. As the flow departs from criticality, the related nonlinear dispersion relation leads to the saturated state of the global mode, including a modification of its frequency [START_REF] Pier | Fully nonlinear global modes in spatially developing media[END_REF][START_REF] Couairon | Absolute and convective instabilities, front velocities and global modes in nonlinear systems[END_REF][START_REF] Chomaz | Against the wind[END_REF][START_REF] Pier | Bifurcation to fully nonlinear synchronized structures in slowly varying media[END_REF].

The second class of wavemakers is based on the nonlinear counterpart of an absolute instability. It considers the behaviour of the front formed in the amplitude of the instability at the boundary between the upstream region of unperturbed base flow and the downstream region of saturated instability. Following [START_REF] Dee | Propagating pattern selection[END_REF], this front can be 'pulled' by the unperturbed flow and propagate upstream until it reaches the boundary between the convectively and absolutely unstable regions where it gets pinned and acts as a wavemaker. As these dynamics are governed by the unperturbed flow, the critical conditions and frequencies of global modes are computed from the linear dispersion relation for absolute instabilities [START_REF] Pier | Steep nonlinear global modes in spatially developing media[END_REF][START_REF] Pier | Bifurcation to fully nonlinear synchronized structures in slowly varying media[END_REF]. According to [START_REF] Van Saarloos | Fronts, pulses, sources and sinks in generalized complex Ginzburg-Landau equations[END_REF], the front could also be 'pushed' by the saturated instability. The present study implements three selection mechanisms introduced by [START_REF] Monkewitz | Global linear stability analysis of weakly non-parallel shear flows[END_REF], [START_REF] Pier | Steep nonlinear global modes in spatially developing media[END_REF] and [START_REF] Healey | Model for unstable global modes in rotating-disk boundary layer[END_REF]. These rely exclusively on the linear local dispersion relation to determine the location of the wavemaker, and the critical condition and frequency of the global mode.

The theoretical approaches above were initially developed in the literature for model scalar equations, including the one-dimensional complex Ginzburg-Landau equation or the two-dimensional vorticity equation for shear flows. Our work fits within a series of works extending these approaches to the three-dimensional Navier-Stokes equations in combination with DNS or experiments in configurations such as rotating disks [START_REF] Pier | Finite-amplitude crossflow vortices, secondary instability and transition in the rotatingdisk boundary layer[END_REF][START_REF] Appelquist | Global linear instability of the rotating-disk flow investigated through simulations[END_REF][START_REF] Appelquist | On the global nonlinear instability of the rotating-disk flow over a finite domain[END_REF]. Global modes in Rayleigh-Bénard-Poiseuille flow have received theoretical [START_REF] Carrière | Transverse-roll global modes in a Rayleigh-Bénard-Poiseuille convection[END_REF], numerical [START_REF] Martinand | Three-dimensional global instability modes associated with a localized hot spot in Rayleigh-Bénard-Poiseuille convection[END_REF]) and experimental [START_REF] Grandjean | Experimental investigation into localized instabilities of mixed Rayleigh-Bénard-Poiseuille convection[END_REF] attention, with excellent agreement. For that case, the thermoconvective instability can be spatially altered by modifying the local temperature between the plates. Global modes in TCP flows are less studied because it is more difficult to spatially vary the centrifugal instabilities. Consider the case of a fixed outer cylinder of radius r 2 and a rotating inner cylinder with radius r 1 and angular velocity Ω. Experimentally, the angular velocity cannot be smoothly varied along the axial direction. Consequently, the local stability must be altered indirectly.

The stability of TCP flow is governed by three control parameters. The first is the Taylor number Ta = Ωr 1 d/ν (or rotating Reynolds number), where ν is the kinematic viscosity of the fluid and d = r 2 -r 1 the gap width. The second is the radius ratio η = r 1 /r 2 . The third is the axial Reynolds number β = Wd/ν, where W is the mean axial velocity. One method of locally increasing the instability of the TCP flow is to locally increase r 2 . This increases the gap width and the Taylor number Ta, while decreasing the radius ratio η and axial Reynolds number β owing to the conservation of the axial mass flow rate. In § 5, we show that these changes can indeed promote centrifugal absolute instabilities. This approach is nevertheless impractical, as it requires a set of different outer cylinders that cannot be tuned easily. Though Taylor-Couette systems with varying radii have been studied experimentally and numerically in the absence of mean axial flow [START_REF] Soward | The linear stability of the flow in the narrow gap between two concentric rotating spheres[END_REF], no results with mean axial flow are available in the literature, to our knowledge.

The current study shows that the local stability of the TCP flow can be altered in a more controlled fashion by replacing one or both of the impermeable cylinders with permeable cylindrical membranes. As fluid exits or enters the annular region through these membranes, the mean axial flow increases or decreases and thus alters the local stability. This configuration for TCP flow is used in an application called rotating filtration, in which a cylindrical semi-permeable membrane rotates within a concentric fixed impermeable cylinder. As suspension is introduced at one end of the annulus, fluid is extracted through the inner cylinder, and concentrate exits the opposite end of the annulus (see [START_REF] Lee | Rotating membrane filtration and rotating reverse osmosis[END_REF], and references therein for applications). This approach reduces the accumulation of contaminants on the membrane surface because the centrifugal instabilities wash them away from the membrane [START_REF] Schwille | Design parameters for rotating filtration[END_REF]. Due to the small permeability of the membrane, the base flow in rotating filtration systems varies slowly in the axial direction, thus opening the door to a global mode analysis.

Previous analytical and numerical studies have considered the related case of TCP flow with uniform radial through-flows applied at the inner and outer cylinders. The radial flows are applied such that the mass flux entering one cylinder equals that exiting the other cylinder. The resulting base flow is axially invariant and amenable to convective/absolute stability analysis (see Martinand et al. 2009, and references therein). Until recently, stability analysis of global modes in rotating filtration has been hindered by analytical and numerical difficulties related to the fact that the filtrate flow through the membrane is driven by the transmembrane pressure difference. Considering the base flow, this leads to a simultaneous coupling between the axial pressure drop and the transmembrane pressure. A generalized analytical solution, sketched in figure 1, has been found, however [START_REF] Tilton | Pressure-driven radial flow in a Taylor-Couette cell[END_REF], as an asymptotic expansion that couples the transmembrane flow and pressure using Darcy's law. Centrifugal instabilities in rotating filtration systems are also challenging to simulate numerically because the membrane introduces an additional coupling between the velocity and the pressure on the permeable cylinder. As a result, traditional pressure solvers for integrating the Navier-Stokes and continuity equations in time can fail. [START_REF] Tilton | A 3D pseudospectral algorithm for fluid flows with permeable walls. Application to filtration[END_REF] have addressed this issue by developing an efficient, non-iterative, projection method tailored to simulating membranes in unsteady three-dimensional fluid flows.

The current study is organized as follows. In § 2, we present the geometry, governing equations, base state and global mode analyses. In § 3, we compare the predictions of our global mode analyses with dedicated DNS. In § 4, we further address the frequency jumps, a specific feature observed in the DNS, and test the robustness of these instabilities by applying more intrusive inlet conditions or higher Taylor numbers to the DNS. After summarizing our results, their implications for rotating filtration and possible future work are then discussed in § 5.

Global modes in axially varying Taylor-Couette-Poiseuille flow

Hereinafter, we use cylindrical coordinates (r, θ , z), in which the velocity vector is denoted V = (U, V, W) t . As sketched in figure 1, we consider an open Taylor-Couette cell with a stationary outer cylinder of radius r 2 and an inner rotating cylinder of thin permeable membrane that is shaded grey in figure 1. The rotating membrane surrounds a cylindrical cavity that is maintained at a constant pressure P ref . The radial fluid velocity U(r 1 , θ, z) through the permeable membrane depends on the transmembrane pressure difference, as described by Darcy's law:

U(r 1 , θ, z) = - κ µh [P(r 1 , θ , z) -P ref ], (2.1)
where κ is the membrane permeability (units m 2 ), h is the membrane thickness and µ is the fluid dynamic viscosity. Note that the transmembrane flow in rotating filtration is typically quite small, such that it flows along a thin film on the inner rotating membrane.

The assumption of constant P ref is consequently well justified, and without loss of generality, we set P ref = 0. Fluid flow in the annular region is governed by the incompressible Navier-Stokes and continuity equations. We non-dimensionalize these equations with respect to the gap width d = r 2 -r 1 , characteristic time d 2 /µ, characteristic pressure ρν 2 /d 2 , and characteristic velocity ν/d, with ρ the density of the fluid. The flow is characterized by a Taylor number Ta, radius ratio η and non-dimensional permeability σ :

Ta = Ωr 1 d ν , η = r 1 r 2 , σ = κ hd . (2.2a-c)
Hereinafter, all expressions are non-dimensional, unless otherwise stated. The governing Navier-Stokes and continuity equations then take the form:

∂V ∂t + V • ∇V = -∇P + ∇ 2 V and ∇ • V = 0. (2.3a,b)
On the outer cylinder, no-slip conditions are applied. On the permeable inner cylinder, the azimuthal and axial components of the velocity also satisfy the no-slip condition, while U is related to the transmembrane pressure difference through Darcy's law:

U(r 1 , θ, z) = -σ P(r 1 , θ, z) and U(r 2 , θ , z) = 0 V(r 1 , θ, z) = Ta and V(r 2 , θ , z) = 0 W(r 1 , θ, z) = W(r 2 , θ , z) = 0.    (2.4)
The no-slip condition on the inner cylinder is guided by practical applications to membrane filtration, which use strongly anisotropic membranes with negligible slip. Focusing on the filtration of solutions, the dimensional ratio κ/h of membranes used in filtration units typically ranges from 10 -10 m for reverse osmosis to 10 -4 m for microfiltration. Though the gap d might vary substantially with the size of the unit, the non-dimensional permeability σ is thus a very small quantity in real systems.

Base flow

The transmembrane pressure difference between the annulus and cavity drives a radial extraction or injection of fluid through the permeable cylinder. This in turn drives radial and axial flows that vary axially in the annular region. An analytical solution for the axially varying steady axisymmetric base flow V b (r, z) has been found by [START_REF] Tilton | Pressure-driven radial flow in a Taylor-Couette cell[END_REF] in the form of an asymptotic expansion about the small parameter radius r 1 . The outer cylinder is impermeable, while the inner cylinder is made of a z-coordinate, or equivalently, a function of a slow variable εz. Based on this solution, a local axial Reynolds number is introduced:

β(z) = W(z) = 1 π(r 2 2 -r 2 1 ) r 2 r=r 1 2π θ=0 W b (r, z)r dθ dr = A 1 exp(z ζ σ ) + A 2 exp(-z ζ σ ), (2.5)
with W b (r, z) as in appendix A and the geometrical factor

ζ = 16 log η (r 4 2 -r 4 1 ) log η + (r 2 2 -r 2 1 ) 2 .
(2.6)

Constants A 1 and A 2 in (2.5) are first order in ε. They can be set by prescribing the axial Reynolds number β(z) at two arbitrary locations z = 0 and z = L.

To characterize the transmembrane flow, we introduce a local radial Reynolds number, related to the axial one β(z) by the conservation of mass:

α 1 (z) = r 1 U b (r 1 , z) = r 2 2 -r 2 1 2 dβ(z) dz = √ ζ σ (r 2 2 -r 2 1 ) 2 [A 1 exp(z ζ σ ) -A 2 exp(-z ζ σ )]. (2.7)
This radial Reynolds number is second order in ε and, consequently, small. Owing to the exponential terms in (2.5) and (2.7), β(z) and α 1 (z) only make sense over a finite axial length, assuming end effects can be ignored. The constants A 1 and A 2 in (2.5) and (2.7) are selected so as to consider three different types of base flows over a finite axial length, as shown in figure 2. Panels (ac) all show an increase along the axial direction of the radial Reynolds number on the inner cylinder α 1 (z) (dashed line). This occurs due to the axial pressure drop driving the mean axial flow. Figure 2(a) demonstrates an 'injection' configuration in which the axial Reynolds number β(z) (solid line) increases monotonically with z because fluid is entering the annulus from the membrane along the entire axial length L of the domain. Figure 2(b) demonstrates an 'extraction' configuration in which β(z) decreases monotonically with z due to the fact that fluid is leaving the annulus through the membrane along the entire axial length. Finally, figure 2(c) shows 'cross-flow reversal' in which the transmembrane flow transitions from extraction to injection.

Global modes and their selection mechanisms

We now proceed to analytically address the centrifugal instabilities as global modes. The velocity and pressure fields are decomposed as

[V, P] = [V b , P b ] + [V p , P p ],
where [V b (r, z), P b (r, z)] is the base flow presented in § 2.1 and appendix A, and [V p (r, θ , z, t), P p (r, θ, z, t)] is a small time-dependent perturbation. Linearizing the problem (2.3)-(2.4) about the base flow [V b , P b ] produces an autonomous system of partial differential equations admitting synchronized solutions, i.e. solutions presenting a common, complex, temporal frequency ω glob such that [V p , P p ] = [V glob p (r, θ , z), P glob p (r, θ, z)] exp(-iω glob t). Though the spatial shape functions [V glob p , P glob p ] can be obtained by solving a multi-dimensional eigenproblem numerically, we leverage the slow axial variation of V b to reduce the numerical cost to a one-dimensional eigenproblem. Based on the weak z-dependence of V b (r, z), the instabilities are expressed in the form of the Wentzel-Kramers-Brillouin (WKB) ansatz of a wavepacket: The leading order (ε 0 ) of the WKB expansion satisfies a local generalized eigenvalue problem

[V p , P p ] = [v glob p (
A(k(z), n, V b,0 (r, z))[v loc p , p loc p ] = ω loc B[v loc p , p loc p ], (2.9)
where the superscripts on eigenvector [v loc p , p loc p ] and eigenvalue ω loc stress that (2.9) depends on the local base flow as discussed below. The operator A is a r-differential operator, parametrized by the axial coordinate z through the wavenumber k(z) and the base flow V b,0 (r, z), the leading order ε 0 of the base flow V b (r, z). To include axial variations at order ε 0 , the base flow reported in appendix A must be slightly modified, so that the axial Reynolds number β(z), given by (2.5), is considered as order ε 0 instead of ε 1 . Referring to the expressions in appendix A, this produces a leading-order base flow that is identical to axially invariant TCP flow, composed of the axial velocity W b = β(z)w b (r), and the leading order of the azimuthal velocity, Tav b,0 (r). Consequently, problem (2.9) amounts to the stability of axially invariant TCP flow with the parameter β set to β(z).

As the convective/absolute stability of this flow is well studied [START_REF] Recktenwald | Taylor vortex formation in axial throughflow: linear and weakly nonlinear analysis[END_REF][START_REF] Martinand | Absolute and convective instability of cylindrical Couette flow with axial and radial flows[END_REF][START_REF] Altmeyer | Islands of instability for growth of spiral vortices in the Taylor-Couette system with and without axial through flow[END_REF], it is only summarized in appendix B and results relevant for the global mode analysis in § 3 briefly presented. Figure 3 shows the results, in the form of the critical Taylor numbers and critical frequencies, of the local convective/absolute stability problem (B 2), as a function of β, for η set to 0.85. Beyond their respective critical conditions, figure 3 also sheds light on some differences and similarities between convective and absolute instabilities. Whereas absolute critical modes are always travelling toroidal vortices (n abs crit = 0), convective critical modes become helical vortices for axial Reynolds numbers above β = 13.8. The critical frequency is a growing function of the axial Reynolds number β for both convective and absolute instabilities. For axisymmetric vortices, absolute and convective instabilities exhibit very similar frequencies. It is therefore difficult to discriminate between convective and absolute toroidal instabilities on the basis of their temporal frequencies.

As the local stability problem (2.9) is parametrized by the axial coordinate z, its solution provides no selection criterion for a global mode frequency ω glob per se.

A selection criterion is usually related to specific axial location(s) from which the wavepacket can propagate both down-and upstream, i.e. to a 'turning point' z t.p. where the group velocity vanishes: ∂ω ∂k (z t.p. ) = 0.

(2.10)

Hereinafter, the (z t.p. )-dependence means that the group velocity is computed from the resolution of the local stability problem (B 2) where the β is set to β(z t.p. ). Condition (2.10) can be satisfied by an absolute mode that is unstable if the system is locally above the critical conditions for absolute instabilities. A turning point then acts as a wavemaker if it imposes the temporal frequency of the local absolute mode, i.e. ω abs (z t.p. ) to the rest of the domain. Critical conditions for the global mode, Ta glob crit , can be then deduced from the growth rate of the absolute mode at the turning point, i.e. the imaginary part of ω abs (z t.p. ).

Several analytical selection criteria for identifying wavemakers are available in the literature. They begin by identifying the axial location of a turning point. We implement two linear selection criteria introduced in [START_REF] Monkewitz | Global linear stability analysis of weakly non-parallel shear flows[END_REF] and extended in [START_REF] Healey | Model for unstable global modes in rotating-disk boundary layer[END_REF]. Both criteria consider a turning point located where the flow is locally the most unstable, i.e. where Im(ω abs (z)) is maximum. Elaborating on figure 3(a), the growth rate of the critical absolute instabilities is a decreasing function of the axial Reynolds number β:

∂ β Im(ω abs crit ) = -∂ Ta Im(ω abs crit )/∂ Ta β abs crit < 0. A local maximum of the growth rate Im(ω abs (z)) therefore corresponds to a local minimum of β(z).
The first selection criterion [START_REF] Monkewitz | Global linear stability analysis of weakly non-parallel shear flows[END_REF] addresses cases where an island of absolute instability is surrounded by a sea of stability, up-and downstream. Though this criterion assumes an infinite axial domain, it also applies to finite axial configurations in which the inlet is locally stable and the outlet is either locally stable or convectively unstable. The axial flow shown in figure 2(c) presents a minimum of β(z) located within the fluid domain, constituting a potentially unstable 'double turning point'. The second selection criterion addresses cases of semi-infinite axial extent that present a maximum absolute instability at the inlet or outlet, as in [START_REF] Monkewitz | Global linear stability analysis of weakly non-parallel shear flows[END_REF] and [START_REF] Healey | Model for unstable global modes in rotating-disk boundary layer[END_REF], respectively. Though this criterion assumes a semi-infinite axial extent, it also applies for systems of finite axial extent in which the outlet or the inlet are locally stable or convectively unstable. Figure 2 demonstrates base flows that present minimums of β(z) at the inlet (panel a) or the outlet (panel b). These constitute potentially unstable 'single turning points'. A linear selection criterion therefore leads to critical conditions in term of a critical Taylor number Ta glob crit , together with characteristics of the global mode at these critical conditions, such as the critical frequency ω glob crit .

Linear selection by a double turning point

Based on the local stability problem at a double turning point z d.t.p. , where the flow is locally the most unstable, the local growth rate also satisfies: ).

∂Im(ω abs ) ∂z (z d.t.p. ) =
(2.13)

The global mode is characterized by the real frequency:

ω loc crit = ω abs crit (z d.t.p.
).

(2.14) Following (4.21a) in [START_REF] Monkewitz | Global linear stability analysis of weakly non-parallel shear flows[END_REF], matching the perturbation around this double turning point to the perturbation in the rest of the domain mathematically imposes a small correction to the frequency in the form: and the correction to the azimuthal velocity dβ/dzv b,1 (r). These components being proportional to dβ/dz, they vanish at the double turning point. The frequency correction in (2.15) is complex, and the sign in the expression is chosen so that Im(ω glob ) < Im(ω loc ), i.e. so that this correction is stabilizing. This stabilizing effect increases the critical Taylor number above which such a global mode can be observed, compared to the leading-order critical value (2.13): (2.17)

ω glob = ω loc + i 1 2 ∂ 2 ω ∂k∂z abs ± 1 2   ∂ 2 ω ∂k 2 abs ∂ 2 ω ∂z 2 abs - ∂ 2 ω ∂k∂z abs 2   1/2 = ω loc ± 1 2 d 2 β dz 2 1/2 ∂ 2 ω ∂k 2 abs ∂ω ∂β abs 1/2 . ( 2
Ta glob crit = Ta loc crit - 1 2 d 2 β dz 2 1/2 Im   ± ∂ 2 ω ∂k 2
Beyond the axial variations of its local wavevector imposed by the dispersion relation with ω set to ω glob crit , the axial variation of the amplitude of the global mode is obtained in the vicinity of z d.t.p. as a hypergeometric function, matched up-and downstream to WKB expansions (2.8). As the perturbation travels downstream, it gains energy as long as it lies in a locally unstable domain. The linear global mode is thus expected to reach its maximum amplitude on the downstream boundary between the convectively unstable and the stable domains.

A nonlinear extension of global modes developing from a double turning point has been proposed (see [START_REF] Pier | Fully nonlinear global modes in spatially developing media[END_REF][START_REF] Pier | Bifurcation to fully nonlinear synchronized structures in slowly varying media[END_REF], yielding a correction to the frequency and the amplitude at saturation of the instability. This extension, nevertheless, is based on the nonlinear dispersion relation and has not been implemented here.

Linear selection by a simple turning point

For base flows without cross-flow reversal, the maximum of the local growth rate, and the corresponding turning point, is expected to occur at the axial boundary of the domain where β(z) is minimum. Figure 2(a) demonstrates that in cases of fluid injection the turning point (marked with a circle) is located at the inlet. ). An instability governed by such a turning point should thus be observed as the Taylor number overcomes the critical absolute Taylor number at the simple turning point:

Ta loc crit = Ta abs crit (z s.t.p.
).

(2.18)

The global mode is characterized by the real frequency:

ω loc crit = ω abs crit (z s.t.p.
).

(2.19) Following (5.17b) in [START_REF] Monkewitz | Global linear stability analysis of weakly non-parallel shear flows[END_REF], and developing on (2.8) in Healey (2010), matching the perturbation around this turning point to vanishing boundary conditions for the velocity of the global mode, at the inlet (injection) or at the outlet (extraction), imposes a correction to this frequency:

ω glob = ω loc - a 0 3 √ 2   ∂ω ∂β abs 2 ∂ 2 ω ∂k 2 abs   1/3 dβ dz 2/3 , (2.20)
where a 0 ≈ -2.3381 is the first zero of the Airy function. Again, though the (z s.t.p. )dependences are dropped for convenience in (2.20), all quantities inferred from the local stability problem and β(z) are to be computed at z s.t.p. . The frequency correction is complex and, together with choosing (dβ/dz) 2/3 real, the cube root in (2.20) is chosen so that Im(ω glob ) < Im(ω loc ), i.e. so that this correction is stabilizing. This stabilizing effect increases the critical Taylor number above which such a global mode can be observed, compared to the leading-order critical value (2.18): (2.22)

Ta glob crit = Ta loc crit + Im   a 0 3 √ 2   ∂ω ∂β
Beyond the axial variations of its local wavevector imposed by the dispersion relation with ω set to ω glob crit , the axial variation of the amplitude of the global mode is obtained in the vicinity of z s.t.p. as an Airy function, matched down-(injection) or upstream (extraction) to a WKB expansion (2.8). With z s.t.p. located at the inlet, as the perturbation travels downstream, it gains energy as long as it lies in a locally unstable domain. The linear global mode is thus expected to reach its maximum amplitude on the downstream boundary between the convectively unstable and the stable domains. With z s.t.p. at the outlet, the maximum amplitude of the linear global mode is at the outlet.

Even though the exact nature of the outlet boundary condition is only a secondorder effect in the present selection criterion, experimental or numerical outlet conditions usually differ from imposing vanishing Dirichlet conditions for the velocity and pressure fields of the perturbation at the outlet boundary, conditions upon which the analytical approach of [START_REF] Healey | Model for unstable global modes in rotating-disk boundary layer[END_REF] is based. This could hamper comparisons between this theoretical construction and experimental set-ups or numerical simulations.

A nonlinear extension of global modes developing from a single turning point at the inlet has been proposed (see [START_REF] Couairon | Absolute and convective instabilities, front velocities and global modes in nonlinear systems[END_REF][START_REF] Chomaz | Against the wind[END_REF], yielding a correction to the frequency and the length over which the amplitude of the instability reaches saturation. As for nonlinear double turning points, this extension is based on the nonlinear dispersion relation and has not been implemented here.

Nonlinear selection by the upstream convective/absolute boundary

The nonlinear selection criterion developed by [START_REF] Pier | Steep nonlinear global modes in spatially developing media[END_REF][START_REF] Pier | Bifurcation to fully nonlinear synchronized structures in slowly varying media[END_REF] is based on the front-forming dynamics of the nonlinear wavepacket. At the convective/absolute boundary, where the local stability evolves from convectively unstable to absolutely unstable as one moves downstream, the amplitude of the global mode is expected to present a steep front that acts as a wavemaker. This front matches an upstream region where the small perturbation is governed by linear dynamics to a downstream region where the perturbation is substantial and governed by nonlinear dynamics. The location of this front, and the selected complex frequency, are actually outcomes of the local linear stability problem and of the existence of a turning point satisfying (2.10) at this boundary. This selection mechanism allows to consider Taylor numbers above the critical conditions Ta At leading order, the frequency of the global mode based on this front, located at convective/absolute boundary z c/a , is then:

ω loc front = ω abs crit (z c/a ), (
∂ 2 ω ∂k 2 abs crit      1/3 2 Im           - 2 ∂ω ∂β abs crit ∂ 2 ω ∂k 2 abs crit      1/3      Im ∂ω ∂β abs crit , (2.25)
and moves the position of the front slightly downstream, to be located at: 3. Global modes in direct numerical simulations

z glob front = z c/a + a 0 dβ dz 2/3 Im       ∂ω ∂β abs crit      - 2 ∂ω ∂β abs crit ∂ 2 ω ∂k 2 abs crit      1/3       Im ∂ω ∂β abs crit      - 2 ∂ω ∂β abs crit ∂ 2 ω ∂k 2 abs crit      1/3 2 . ( 2

Numerical method

We perform DNS of centrifugal instabilities in a rotating filtration system using the multidomain pseudospectral method of [START_REF] Tilton | A 3D pseudospectral algorithm for fluid flows with permeable walls. Application to filtration[END_REF]. The governing Navier-Stokes equations are integrated in time using a second-order projection method that applies the Darcy boundary condition (2.1) on the membrane exactly at each time step. Together with boundary conditions (2.4) on the cylinders, laminar annular Taylor-Couette-Poiseuille flows are applied at inlet and outlet, prescribed by the global Taylor number Ta and the local axial Reynolds numbers β inlet and β outlet . To accommodate these inlet and outlet boundary conditions, the membrane permeability σ is forced smoothly to zero in buffer regions adjacent to the inlet and outlet, as illustrated in figure 4. The governing equations are discretized spatially using a multidomain spectral method using Chebyshev polynomials in the radial and axial directions, and in Fourier modes in the azimuthal direction.

Perturbations to the base flow occur from a combination of the numerical noise occurring throughout the system at all time steps and a controlled perturbation applied at the inlet on the first four azimuthal modes of the axial component of the velocity (see ( 52) and (53) in [START_REF] Tilton | A 3D pseudospectral algorithm for fluid flows with permeable walls. Application to filtration[END_REF]. In all cases but figure 15(b), the controlled perturbation is smoothly ramped up and down in time, starting at the initial time step and over a short duration 5/Ta (see ( 54 Typical results for simulations, performed at η = 0.85 and Ta = 120 with a short initial perturbation, are shown in figures 5 (fluid injection), 8 (fluid extraction) and 11 (cross-flow reversal). Panels (a) show snapshots, taken at the final simulation times, of the isosurfaces where the azimuthal velocity is 50 % that of the rotating cylinder, i.e., 0.5r 1 Ω. The three configurations show the development of instabilities in the form of axisymmetric vortices. Panels (b) show spatio-temporal diagrams of the radial velocity at the mid-gap U(r mid , z, t). These show that in all three cases, the systems achieved saturated steady states in which vortices are present in only portions of the domains. In figures 5 and 11, the vortices tend to organize themselves in several packets, with different temporal frequencies. Following saturation, we compute these frequencies from the spatio-temporal diagrams by performing fast Fourier transforms (FFTs) in

Comparisons between numerical and analytical results

In this section, we use the global mode selection criteria presented in § 2 to analyse DNS results of systems characterized by fluid injection (figures 5-7), fluid extraction (figures 8-10) and cross-flow reversal (figures 11-13). For each case, we begin by extracting the axial variation of β(z) from the DNS. To this we then fit the analytical expression (2.5) by evaluating the constants A 1 and A 2 using the mean axial flow in the DNS at the boundaries between the inlet and outlet buffers and the central subdomain. Panels (c) of figures 5 (fluid injection), 8 (fluid extraction) and 11 (crossflow reversal) show there is excellent agreement between the DNS (dashed light grey, green online, curves) and expression (2.5) (solid dark grey, blue online, curves). The only deviation occurs in the inlet/outlet buffers (marked by solid vertical black lines in figures 5-13) where the permeability σ of the DNS is smoothly set to zero.

Linear selection criteria are based on turning points located at minima of the axial Reynolds number β(z). For fluid injection (extraction), the minimum of β(z) is located at the inlet (outlet), where, despite the buffer, the value obtained from expression (2.5) only minutely differs from the mean axial flow in the DNS. For cross-flow reversal, the location and value of the minimum of β(z) are found from expression (2.5). As the permeability σ does not appear in the local stability problem (B 2), the presence of the inlet/outlet buffers does not impact the local linear analysis of the perturbations, beyond setting β(z). The whole analysis is therefore based on the knowledge of the radius ratio η, the permeability σ , the mean axial flow at two locations and the Taylor number, together with resolutions of the local convective/absolute stability problem (B 2). Boundaries where the flow transitions from stable to convectively unstable and from convectively unstable to absolutely unstable are found by solving the axially invariant stability problem (B 2) with Ta as in the DNS, and β so as to be at critical conditions for convective or absolute instabilities. Boundaries are then found as z such that β(z) in (2.5) is equal to the critical β. In figures 5-13, stable regions are labelled a convectively and an absolutely unstable domain by a dash-dotted dark grey (blue online) vertical line in figures 5-13.

Fluid injection: maximum instability at the inlet

We begin by considering systems with fluid injection along the full length of the membrane. As demonstrated in figures 2(a) and 5(c), fluid injection causes the axial Reynolds number, β(z), to increase monotonically with z, such that the local stability increases as fluid moves downstream. In the DNS, we impose fluid injection by applying β inlet = 2 and β outlet = 40, for a permeability σ = 5 × 10 -5 and L central = 120. Simulations are repeated for Ta = 112,114,115,120,130,140 and 150. Figure 5 shows DNS results for Ta = 120. Panels (a,b) demonstrate that an axisymmetric instability develops and propagates from the inlet until it reaches a saturated state in which vortices are present over most of the domain, except near the outlet. From the spatio-temporal diagram in figure 5(b), we observe that the saturated state is characterized by vortices organized in several packets with different frequencies.

For fluid injection cases, our global mode analysis predicts a linear selection mechanism based on a simple turning point located at the inlet z inlet , as addressed in § 2.2.2. At leading order, based on the values of η = 0.85 and β inlet = 1.88, condition (2.18) and the resolution of the local stability problem (B 2) for absolute instability, with β set to β inlet leads to a critical value Ta loc crit = 111.2, above which a global mode becomes unstable. Equation (2.19) predicts a corresponding critical frequency ω loc crit = 6.87. Taking the corrections at the next order into account, expressions (2.21) and ( 2 Figure 6 shows the axial variations of the wavepacket in the DNS slightly above the critical conditions, at Taylor number Ta = 114. Figure 7 further departs from these critical conditions at Ta = 130. Global instabilities originating at the inlet occur in both cases. Near the inlet, these modes have the frequencies ω = 7.50 for Ta = 114 and ω = 7.55 for Ta = 130. This is in agreement with the linear selection mechanism governed by a simple turning point at the inlet, leading to ω glob crit = 7.40 at critical conditions, as mentioned above. Nonlinear effects are expected in the dynamics of the instability, because the Taylor number departs from its critical value for a global mode, Ta glob crit , in both The asymmetry between the inward (U(r mid , z, t final ) < 0) and outward (U(r mid , z, t final ) > 0) phases of the radial velocity observed in figures 6(b) and 7(b) tends to confirm the existence of noticeable nonlinearities that could also impact the frequency of the instability in the DNS. Such a modified frequency could be obtained using the nonlinear extension of semi-infinite global modes proposed in [START_REF] Couairon | Absolute and convective instabilities, front velocities and global modes in nonlinear systems[END_REF]. The frequency, as measured in the leading subpacket of the DNS, nevertheless remains very close to the result of the linear selection mechanism.

Moving downstream from the inlet, the local dominant frequencies of the wavepackets in the DNS, ω max (z) shown in figures 6(a) and 7(a) exhibit abrupt changes in a sequence of step-ups that are not accounted for by the global mode analysis. As indicated by the values in the captions of figures 6 and 7, these frequency jumps correspond to frequency doublings and occur until the wavepacket enters the locally stable domain and diminishes. Comparing the frequency jumps in figures 6(a) and 7(a) with the Hilbert transform in figures 6(b) and 7(b), we observe that the perturbation amplitude decreases as it approaches a frequency jump, and then increases downstream from the jump. These results suggest that the frequency jumps allow perturbations to resettle to more natural frequencies. This interpretation is supported by noting that the jumps occur both in the locally absolutely and convectively unstable regions, which suggests that the jumps are related to the local convective behaviour of the stability problem. The dotted dark grey (blue online) curves in figures 6(a) and 7(a) show the frequency ω loc (z) of the most unstable axisymmetric convective mode obtained for the local axial Reynolds number β(z) and the respective Taylor numbers. We see that frequency doublings occur at locations such that the local coarsely adapts to ω loc (z). This feature will be further addressed in § 4.1. 113, 115, 120, 130 and 140. Figure 8 shows DNS results for Ta = 120. We observe in figure 8(b) that an axisymmetric instability originates at the outlet. A wavepacket then propagates upstream until it reaches a steady state, within the region where the flow is locally absolutely unstable (downstream of the dash-dotted dark grey, blue online, vertical line). The spatio-temporal diagram in figure 8 In this configuration, global modes can also be governed by the nonlinear front-forming mechanism described in § 2.2.3. This possibility arises because the flow transitions from convectively unstable to absolutely unstable, at the axial locations z c/a indicated by the dash-dotted dark grey (blue online) vertical lines in figures 8, 9 and 10. At leading order, this mechanism selects the frequency ω loc front as the critical absolute frequency at z c/a , as given by (2.24). Implementing the next-order correction leads to the corrected frequency ω glob front computed from (2.25). This also shifts the front slightly downstream from z c/a to z glob front , according to (2.26). Both the modified front location and selected frequency are then functions of the Taylor number and are shown in dash-dotted light grey (green online) vertical and horizontal lines, respectively, in figures 9 and 10.

3.2.2.

Slightly above the critical conditions, at Taylor number Ta = 113 (figure 9), the very good agreement between the linear frequency ω glob crit = 7.55 and the numerical result ω = 7.69, together with the propagation of the perturbation upstream of the modified front , suggest that linear dynamics prevail for these weakly supercritical conditions.

Increasing the Taylor number further from critical conditions, to Ta = 130 (figure 10), together with unshown results at Ta = 115, 120 and 140 exhibits the role of the nonlinear front-forming mechanism in the selection and dynamics of the global mode. The upstream edge of the perturbation in figure 10(b) almost exactly collapses on z glob front (the dash-dotted light grey, green online, vertical line in figure 10a) as given by (2.26). The frequency ω max (z) = 18.49 of the leading subpacket agrees very well with ω glob front = 17.91 (the dash-dotted light grey, green online, horizontal line in figure 10a) as given by (2.25).

Whereas the DNS at Taylor number Ta = 113 in figure 9 shows a synchronized global mode, the DNS at Ta = 130 in figure 10 shows that the local dominant frequency of the perturbation, ω max (z) in figure 10(a), changes along the axial direction. This change takes the form of two successive jumps to smaller frequencies as fluid travels downstream across the wavepacket. As in § 3.2.1, this change is again associated with a 'bump' in the Hilbert transform shown in figure 10(b). This suggests again that the frequency jumps allow perturbations to resettle to more natural frequencies, depicted as the dotted (blue online) curves which show the frequency ω loc (z) of the most unstable axisymmetric absolute mode obtained for the local axial Reynolds number β(z) and the Taylor number of the DNS. Frequencies observed in figure 10 and for similar cases at Ta = 120 and 140 suggests that the frequency downstream of the step-down is about two-thirds of the frequency upstream of the step-down. Section 4.1 will discuss these step-downs further.

Cross-flow reversal: maximum instability within the domain

Our final configuration considers systems with cross-flow reversal, as in figures 2(c) and 11(c). The axial Reynolds number first decreases with z, then increases after the transmembrane flow reverses. This set-up exhibits a maximum of local instability where β = β min is a minimum. We impose cross-flow reversal in the DNS by applying β inlet = β outlet = 31.2, with a permeability σ = 2.3 × 10 -4 and L central = 130. Simulations are repeated for Ta = 110,112,115,120,130 and 140. Figure 11 shows DNS results for Ta = 120. An axisymmetric instability develops in the form of a wavepacket originating slightly downstream from the location of β min (the dashed vertical dark grey, blue online, line), as seen in figure 11(b). This wavepacket propagates both down-and upstream until its trailing edge reaches the downstream stable domain (the dotted vertical dark grey, blue online, line on the right), where the instability diminishes, and its leading edge gets pinned downstream of the convective/absolute boundary (the dash-dotted vertical dark grey, blue online, line on the left). Figures 11(b,c) show that the steady wavepacket is again composed of several subpackets with different temporal frequencies.

A region of local absolute instability can exist around the location of β min . This location, where β is minimum and the flow locally the most unstable, constitutes a double turning point. At leading order, from the minimum axial Reynolds number β min = 1.61, the linear selection mechanism driven by the double turning point produces a critical Taylor number (2.13) Ta loc crit = 110.5 and a critical frequency (2.14) ω loc crit = 5.95. Including next-order corrections (2.16) and (2.17) shifts these slightly to Ta glob crit = 110.9 and ω glob crit = 6.04. The threshold for the instabilities in the DNS is in agreement with these analytical values. No instability is observed at Ta = 110, whereas increasing the Taylor number to Ta = 112 exhibits an instability, as seen in figure 12. front location z glob Owing to the existence of the upstream convective/absolute boundary z c/a , the selection by the front-forming mechanism, described in § 2.2.3, may also occur, in a similar fashion to the cases of fluid extraction. As in § 3.2.2, both the modified front location, as given by (2.26), and selected frequency, as given by (2.25), are then functions of the Taylor number and are shown in dash-dotted light grey (green online) vertical and horizontal lines, respectively, in figures 12 and 13.

Slightly above the critical conditions, it is difficult to discern the selection mechanism of the global mode from the numerical simulation at Ta = 112 in figure 12. The associated frequencies ω glob front = 6.00 by the nonlinear front and ω glob crit = 6.04 by the linear double turning point are nearly identical when Ta = 112, and are both close to the observed value ω = 5.98 for the leading subpacket. The fact that this subpacket does not exhibit a clear front could suggest a linear global mode, selected by the double turning point. The amplitude of such a linear global mode should grow until it reaches the downstream boundary between convectively unstable and stable domains, but in our case, such a growth is hampered by the local detuning and the frequency jump.

Increasing the Taylor number to Ta = 130 in figure 13, together with unshown similar DNS at Ta = 115, 120 and 140, shows that the leading subpacket of the instability tends to be governed by the nonlinear front-forming mechanism, with this subpacket pinned on the shifted location z glob front (the dash-dotted light grey, green online, vertical line in (a)) and the frequency of the subpacket ω = 16.84 closely matching the analytical value ω glob front = 16.47 (the dash-dotted light grey, green online, horizontal line in (a)) and clearly departing from ω glob crit = 6.04. As in § § 3.2.1 and 3.2.2, the wavepackets at Ta = 112 and 130 exhibit jumps in the dominant temporal frequency, ω max (z) in (a), as fluid travels downstream. These jumps are similar to those observed in cases of fluid extraction, with step-downs upstream of the cross-flow reversal, and fluid injection, with step-ups downstream of the cross-flow reversal. As for fluid injection, figures 12 and 13, together with unshown results at Ta = 115, 120 and 140, show that step-ups consist of frequency doublings between the successive frequencies. In a fashion similar to fluid extraction, step-downs observed in figure 13 tend to follow a two-thirds ratio between the successive frequencies, though this is not as rigorous as the frequency doublings of the step-ups, and the lowest frequency ω = 6.95 is close to the frequency selected by a double turning point at the minimum of β, ω glob crit = 6.04. Those jumps are further addressed in the next section.

Discussion

The frequency jumps and the subpackets

For almost all cases, our DNS exhibit instabilities composed of several subpackets, oscillating at different temporal frequencies that occur in sequences of step-ups or step-downs. Step-ups are observed in regions of fluid injection as in the simulations of § 3.2.1, and in regions downstream of cross-flow reversal in § 3.2.3. Step-downs are observed in regions of fluid extraction, as in the simulations of § 3.2.2 and in regions upstream of cross-flow reversal in § 3.2.3. To investigate the dynamics of the subpackets and frequency jumps observed in our DNS, we consider the case of crossflow reversal shown in figure 11 for Ta = 120. This set-up simultaneously presents step-ups and step-downs of the frequency. Further information can be obtained from figure 14, which shows the complete Fourier transform in time of the instability, upon which the computation of ω max (z) is based, as a function of the axial coordinate z. As seen in figure 14, as well as in figures 6-7 and 12-13, step-ups clearly consist of successive frequency doublings. These frequency doublings occur irrespective of whether the upstream subpacket with the lowest frequency was imposed by a linear turning point, as in § 3.2.1, or if it resulted from a sequence of step-downs upstream, as in the cross-flow reversal cases of § 3.2.3. This suggests that the frequency doubling is governed by some local dynamics. In figure 14, for example, the step-up between the two subpackets pulsating at frequencies ω 2 and 2ω 2 occurs due to the local amplification in the subpacket downstream of the step-up of the first harmonic (2ω 2 ) of the frequency (ω 2 ) of the subpacket upstream of the step-up. The first harmonic (2ω 2 ) is fed by nonlinear dynamics and can be already observed in the subpacket upstream of the step-up. The frequency step-up occurs where this first harmonic is locally amplified by the mechanism of the centrifugal instability. This is shown by examining the local frequency ω loc (z), the dotted, blue online, curves in figures 14. Similar curves are also depicted in figures 6, 7, 12 and 13. This local frequency consists, in the domains where fluid injection occur and step-ups occur, of the frequency of the most unstable axisymmetric convective mode obtained at the local value of the axial Reynolds number β(z) and the Taylor number of the DNS. The FFT depicted in figure 14 shows that, in the region of fluid injection, the maximum amplitudes of the successive harmonics are located along this curve of local frequency. As seen in figures 6, 7 and 12, in cases of fluid injection and cross-flow reversal close to the critical conditions, the instabilities consist in sequences of step-ups. It seems then legitimate to consider the complete wavepacket as a global mode: the leading subpacket subdues the successive downstream subpackets and, moreover, subpackets are synchronized, as harmonics are synchronized to their fundamental.

While step-ups can be explained in terms of local dynamics, a physical interpretation for the step-downs is more elusive. Step-downs occur when the leading subpacket is governed by the nonlinear front-forming mechanism at the convective/absolute boundary. In these cases, a linear turning point is also present downstream of the front, where the axial Reynolds number β is minimum. In the simulations shown in figures 8, 10, 11, 13 and 14, the leading subpacket pulsates at the frequency ω glob front selected by the front (labelled ω 1 in figure 14). The leading subpacket is superseded by a subpacket pulsating at a frequency close to ω glob crit selected by the linear turning point (labelled ω 2 in figure 14). Furthermore, step-downs tend to follow a two-third ratio between the successive frequencies, though this ratio is not as rigorous as the frequency doublings observed in the step-ups. The step-down in figure 14 shows that the spatial overlap between the two subpackets involved, i.e. the z-range where both frequencies are noticeably present, is limited. It is therefore not clear if, across a step-down, the upstream subpacket drives the downstream subpacket. This motivates us to compare the local maximum frequency ω max (z) in regions where step-downs occur to a local frequency ω loc (z) defined as the frequency of the most unstable axisymmetric absolute mode, obtained at the local axial Reynolds number β(z) and the Taylor number of the DNS. This local frequency is shown as a dotted curve (blue online) in figures 10, 13 and 14. We observe that a new subpacket develops after a step-down where its frequency becomes the locally most unstable. In figure 14, the local frequency ω(z) of the DNS is compared to the analytical result ω loc (z). This latter is the frequency of the most unstable axisymmetric absolute mode where extraction and step-downs occur, while it is the frequency of the most unstable axisymmetric convective mode where injection and step-ups occur. This choice reflects the hypotheses that the perturbation downstream of a step-up is excited by the perturbation upstream, while the perturbation downstream of a step-down seems to arise irrespective of the perturbation upstream. Though we found no convincing mechanisms that clearly explain the selection of the frequency drop in step-downs, we report some possibilities. First, we note that step-downs cannot be unambiguously explained by a downstream subpacket governed by a linear turning point superseding an upstream subpacket governed by the nonlinear front-forming mechanism. Figure 13, for instance, exhibits two successive step-downs involving three subpackets, one downstream of the front, one near the double turning point, and an intermediate one, the dynamics of which is unclear. A second hypothesis is motivated by the observation that successive frequencies in a step-down (labelled ω 1 and ω 2 in figure 14) are in a ratio close to two-thirds. This ratio suggest that the first harmonic (2ω 1 ) of the upstream subpacket excites the downstream subpacket via its second harmonic (3ω 2 ). Unfortunately, this ratio is only loosely matched and figure 14 reveals that this common harmonic 2ω 1 ≈ 3ω 2 is weak. Third, a step-down could suggest that the locally most unstable absolute mode occurs along another branch of the dispersion relation. But the frequency of such a new branch would not conform to ω loc (z). Moreover, to date, no physically acceptable extra solution branches of absolute modes have been found for the dispersion relation. We cannot rule out, however, that an acceptable solution could exist, because finding extra pinch points in the complex k-plane is difficult without an analytical expression for the dispersion relation. Finally, a step-down could be the symptom of a secondary instability. However, secondary instabilities in Taylor-Couette systems usually take the form of vortices presenting a waviness in the azimuthal direction, and the frequency of the subpacket downstream of a step-down would again not conform to the local frequency ω loc (z) but to a different dispersion relation. Moreover, wavy vortices resulting from a secondary instability mechanism are observed in our DNS at higher Taylor numbers (see figure 16). Though we incline towards considering that the dynamics underlying the step-downs are intrinsic, i.e. the feature does not depend on external factors, they remain an open question. Related to these unsettled dynamics, to identify the whole wavepacket as a global mode is debatable.

The selection of the leading wavepackets

Besides the frequency jumps that cannot be explained in terms of global modes, the analysis presented § 2.2 leads to reliable and interesting results, summarized here. The thresholds in terms of Ta crit , beyond which instabilities develop in the DNS in § 3, agree with the analytical expressions based on a linear turning point, whether a simple turning point, as developed in § 2.2.2, in the case of fluid injection shown in § 3. As the Taylor number departs from the critical conditions, as in figures 10 and 13, the nonlinear selection by the front-forming mechanism accurately predicts the frequency and the location of the front of the leading subpacket in cases of fluid extraction and cross-flow reversal. This agreement, however, requires the corrected expressions (2.25) and (2.26). This mechanism is fully nonlinear, and the agreement between the DNS and the analysis actually improves as the Taylor number departs from its critical condition As the Taylor number increases, the front approaches to the absolute/convective boundary z c/a and the frequency approaches the leading-order expression (2.24).

The location and extension of the vortices in the domain can also be explained by the local analysis. In the case of fluid injection, the vortices develop between the inlet and the downstream boundary between the locally convectively unstable and stable regions (the dark grey, blue online, dashed vertical lines between the AU regions) and the outlet. In the case of cross-flow reversal, assuming again that the leading subpacket is selected by the front-forming mechanism, the vortices develop between the upstream boundary between the convectively and absolutely unstable regions (the dark grey, blue online, dashed vertical lines between the

How robust are the global modes?

In the DNS, the manner in which the flow is perturbed impacts the development of the vortices. In § 3, part of the perturbation is controlled as initial and/or boundary conditions, in the form of an initial non-axisymmetric perturbation of the axial component of the inlet velocity. The perturbation also includes numerical noise. Though very limited within the numerical subdomains of the multidomain simulations, some noticeable noise is expected at the boundaries between the subdomains, such as those between the inlet and outlet buffer and the central subdomain. Based on these different perturbations, we extrapolate that the initial perturbation at the inlet triggers the vortices in the case of fluid injection in § 3.2.1. The noise at the outlet triggers the vortices in the case of fluid extraction in § 3.2.2. The numerical noise in the whole domain triggers the vortices in the case of cross-flow reversal in § 3.2.3. The low level of this noise explains the longer time required to observe the instability in figure 11, compared to figures 5 and 8.

Different initial and/or boundary conditions could potentially lead to the selection of flow structures different from those we observed. Whereas the dynamics of the vortices observed in § 3 are based on local absolute instability, we now consider the possibility to observe structures driven by convective instabilities, in a fashion similar to the 'convective route' proposed to explain the transition to turbulence in the boundary layer over a rotating disk (see Imayama, Alfredsson & Lingwood 2016, and references therein on this topic). In our case, these convective modes are conceived of as a response to a forcing in the form of a non-axisymmetric continuous random noise on the axial component of the velocity at the inlet, instead of an initial pulse. In the case of fluid extraction and cross-flow reversal, as shown in figures 8-10 and 11-13, respectively, the inlet is located in a locally stable domain, isolating the dynamics of the vortices from the perturbations at the inlet. We confirmed that changing the perturbation at the inlet from an initial pulse to a continuous random noise had no impact on the instabilities. In the case of fluid injection, as shown in figures 5-7, the minimum axial Reynolds number β inlet = 2 at the inlet implies that the location is convectively unstable as soon as Ta > Ta conv crit ≈ 108.5. Changing the perturbation at the inlet from an initial pulse to a continuous random noise results in identical steady wavepackets, similarly decomposed into sets of subpackets oscillating at the same frequencies ω max (z). The initially forced and continuously forced cases only differ by the fact that, in the latter case, the global mode rises from a substantial background noise resulting from the forcing of convective instabilities. The convective modes are therefore present in the flow, but the dynamics of the global mode prevail.

The possibility of triggering convective instabilities is further addressed in the case of fluid extraction at Ta = 130, by imposing a lower axial Reynolds number at the inlet β inlet = 20, so that the inlet is now locally convectively unstable and by changing the inlet condition from an initial pulse to a continuous random noise. Figure 15 shows the results of these numerical simulations, with an initial perturbation (figure 15a) or with a continuous random noise (figure 15b) at the inlet. Whereas an initial perturbation evolves towards a leading subpacket governed by the front-forming mechanism, similar to the case shown in figure 10, a continuous noise leads to a more complex final steady wavepacket, the dynamics of which are hardly intelligible at this point, that could consist of convective instabilities upstream of the wavepacket. Despite these different upstream dynamics, the two downstream-most subpackets in both cases are of comparable axial extensions and frequencies, with ω = 11.90 and 8.28 for an initial perturbation (figure 15a) and ω = 12.18 and 8.30 for a continuous noise (figure 15b), probably denoting similar downstream dynamics related to the frequency step-downs addressed in § 4.1.

To investigate the sensitivity of the global modes to the numerical treatment of the DNS outlet conditions, we repeated all simulations shown in the current using a 'sponge layer' that attenuates perturbations in a near-outlet region (see Tilton et al. 2014, for details). The results were qualitatively similar, except in cases wavemaker. In this DNS, the flow was found to directly transition at Ta = 117 to a global mode governed by a front forming at the boundary between the convectively unstable and absolutely unstable domains. Modifications of the outlet are known to impact the dynamics of instabilities in rotating-disk boundary layers [START_REF] Pier | Transition near the edge of a rotating disk[END_REF][START_REF] Appelquist | On the global nonlinear instability of the rotating-disk flow over a finite domain[END_REF], a configuration where the local instability increases with the radius, in a fashion similar to the increase of the instability along the axis in TCP with fluid extraction. The impact of the nature of the outlet on the instability in TCP flow would require a substantial and thorough study.

Another aspect of the robustness of these axisymmetric global modes is their stability with respect to azimuthal waviness. Axisymmetric Taylor vortices, stationary or travelling, are expected to acquire some azimuthal waviness as the Taylor number is increased (see Lueptow, Docter & Min 1992, for experimental results). In a Taylor-Couette flow without axial and radial through-flows, wavy vortices are expected above Ta wavy crit ≈ 120, for η = 0.85. In an axially varying configuration of fluid injection, these wavy vortices are observed in the numerical simulation at Taylor numbers Ta = 160, shown in figure 16. In the configurations of fluid extraction and cross-flow reversal, these wavy vortices have not been observed at Ta = 150, and no DNS at higher Taylor numbers are available so far. This discrepancy in the critical conditions for the occurrence of the azimuthal waviness could be accounted for by the finite axial extension of the vortex stack, more limited in the cases of fluid extraction and cross-flow reversal, that could hamper the azimuthal waviness, though this is only speculative.

Conclusions and outlook

We investigated hydrodynamic instabilities developing in an open Taylor-Couette cell with a permeable inner cylinder. Our methodology considered both global mode theory and DNS using high-order spectral methods. The study was motivated by applications to membrane separations and by elucidating the physics of instabilities developing over spatially varying flows in open systems. The global mode analysis proved challenging because the selection criteria available in the literature have seen little application to the full three-dimensional Navier-Stokes equations (see § 2.2). The numerical simulations were also challenging due to the membrane Darcy conditions on the permeable cylinder and the high aspect ratio domains (see § 3.1 and [START_REF] Tilton | A 3D pseudospectral algorithm for fluid flows with permeable walls. Application to filtration[END_REF]. Our results nevertheless show that a global mode analysis based on the local linear dispersion relation successfully explains the critical conditions and locations of the instabilities observed in the DNS (see § 4.2). However, our DNS also clearly and reproducibly demonstrates a new feature of global modes to which we found no previous record in the literature. Specifically, the DNS produced global modes composed of stacked subpackets that are each synchronized at specific frequencies related to the spatial evolution of the most unstable local mode (see § 3.2). Cases in which the frequency steps up as one moves downstream between two subpackets can be explained by the local amplification of nonlinearly driven harmonics. A physical explanation for cases in which the frequency steps down, however, remains elusive (see § 4.1). We expect this new result will motivate future numerical and theoretical work to confirm and explain similar behaviours in other open systems with global modes.

As a final point and outlook, we consider how the mechanics of these instabilities could improve the design of rotating filtration devices. The set-up of such devices is similar to the case of fluid extraction, shown in § 3.2.2. The efficiency of rotating filtration is based on the idea that centrifugal instabilities, in the form of vortices, developing in the gap between the cylinders, will promote the mixing of the substance (whether a solute or a suspension) and thereby reduce the accumulation of the filtered retentate near the semi-permeable membrane surface. This reduces membrane fouling and improves performances. Crucial to the efficiency of rotating filtration devices is the transport and mixing of the solute or suspension by these vortices. These transport and mixing properties depend on the radial shape functions, temporal frequencies, wavenumbers and amplitude of the vortices in a complex fashion which is a current research topic for axially invariant configurations.

Owing to the large number of physical and geometrical parameters, a parametric study of the flow and related transport and mixing in axially varying set-ups remains a tremendous task, if addressed experimentally or numerically. Analytical approaches can offer a way to reduce this task. Together with local convective/absolute stability analyses, the analytical framework of global modes helps foresee some aspects of the dynamics of the vortices, such as their critical conditions, frequencies and wavenumbers, and in the case of global modes, their axial location and extension. These characteristics are determined by controlled parameters such as the solvent viscosity, the permeability and rotation rate of the inner membrane, the inlet feed rate of the device and the working pressure inside the inner cylinder. By setting the base flow and performing a local linear stability analysis as presented in appendix B, the flow configuration is amenable to parametric studies.

Present DNS and analyses in § 3 have covered cases with η = 0.85. To illustrate the interest of a parametric study, figure 17 shows the effect of varying the radius ratio η on the convective (figure 17a) and absolute (figure 17b) critical Taylor numbers, Ta conv crit and Ta abs crit , for the centrifugal instabilities, parametrized by the axial Reynolds number β. Results for convective instabilities are reproduced from [START_REF] Martinand | Weakly nonlinear analysis of cylindrical Couette flow with axial and radial flows[END_REF]. By favouring the occurrence of vortices, moderate radius ratios (η ∼ 0.5) could increase the efficiency of filtration devices. Besides the critical Taylor numbers, the critical azimuthal wavenumbers n conv crit and n abs crit are also indicated in figure 17 by the colourscale; non-zero values denoting the selection of helical travelling vortices instead of toroidal ones. Depending on their respective efficiencies for mixing and transport, which remain to be assessed, fine tuning the physical parameters to select toroidal or helical vortices could be of practical interest for filtration. Figure 17 (b) shows that for most of the radius ratios, the critical absolute instabilities remain toroidal. Helical absolute vortices, though, can be selected in a narrow band of the ηβ space, where η is close to unity, i.e. in narrow gaps, where right-handed helical modes (n abs crit = -1) are selected. In contrast, convective instabilities of figure 17(a) are found to select left-handed helical modes with azimuthal wavenumbers up to n conv crit = 22 as β increases over the same range. From the aspect of the azimuthal wavenumber, one might conclude that favouring convective instabilities instead of absolute ones present some interest.

Besides the selection of convective, absolute or global modes, DNS presented in § 3.2.2 show that the stack of vortices decomposes into several subpackets, with different temporal frequencies, axial wavenumbers and phase speeds and, consequently, different transport properties. Understanding and predicting the frequency jumps would be of clear practical interest, either to account for them in models of rotating filtration or to design systems that avoid them. We are currently working along these lines. Moreover, predicting the amplitude at saturation of the vortices remains a complex and open question. This amplitude varies in the axial direction, is impacted by the frequency jumps, and results from nonlinear dynamics that have not been address to date.

  FIGURE 1. (a) Radial and axial velocities of the laminar base flow for the axially varying problem in a meridional plane. The outer cylinder is impermeable while fluid flows across the inner permeable cylinder, according to Darcy's law, leading to a base flow varying along the axial direction. Panel (b) shows the azimuthal velocity of the base flow in an equatorial plane.

√FIGURE 2 .

 2 FIGURE 2. Three possible base flow variations of the axial Reynolds number β(z) (solid lines) and radial Reynolds number α 1 (z) (dashed lines). (a) Results obtained when β inlet = 2, β outlet = 40 and σ = 5 × 10 -5 . The axial base flow β(z) increases monotonically as fluid is injected along the membrane. The circle at the inlet marks the most unstable simple turning point as defined in § 2.2. (b) Results obtained when β inlet = 40, β outlet = 2 and σ = 5 × 10 -5 . The axial base flow β(z) decreases monotonically as fluid is extracted along the membrane. The circle at the outlet marks the most unstable simple turning point as defined in § 2.2. (c) Results obtained when β inlet = 20, β outlet = 20 and σ = 2.3 × 10 -4 . The base flow β(z) undergoes cross-flow reversal as α 1 (z) transitions from extraction to injection at z = 0. The circle at the centre marks the most unstable double turning point as defined in § 2.2.

  FIGURE 3. (Colour online) Local linear stability results for the critical convective (solid curves) and absolute (dashed curves) modes of instability, as functions of the axial Reynolds number β, for radius ratio η = 0.85. (a) Critical Taylor numbers Ta conv crit and Ta abs crit . (b) Critical frequencies ω conv crit and ω abs crit . Digits in (a) indicate the azimuthal wavenumbers of the critical convective modes n conv crit .

  .15) Though the (z d.t.p. )-dependences are dropped for convenience in (2.15), all quantities inferred from the local stability problem and β(z) are computed at z d.t.p. . Equation (4.21a) in Monkewitz et al. (1993) also includes a term due to the corrections to the base flow at order ε 1 . This term involves here the radial velocity U b = dβ/dzu b (r)

  Ta abs crit ).

  Figure 2(b) demonstrates that in cases of fluid extraction the turning point is located at the outlet. In the case of fluid injection (extraction), the simple turning point is located at the inlet (outlet), z s.t.p = z inlet (z s.t.p = z outlet ), as illustrated in figure 2(a) (injection) and figure 2(b) (extraction), where z s.t.p. are indicated by circles.At leading order, the frequency of a global mode set by a simple turning point is ω loc = ω abs (z s.t.p. ), computed by solving the local stability problem (B 2) for the most unstable absolute instability with β set to β(z s.t.p.

  glob crit , as increasing Ta results in moving the boundary upstream. Using this selection criterion thus leads to the global frequency, as a function of the Taylor number ω glob front (Ta). At arbitrary Taylor number Ta, the location of the convective/absolute boundary z c/a is found by solving the local stability problem (B 2) with β set such that Ta = Ta abs crit . The location z c/a is then found so that β(z c/a ) = β, and is a function of the Taylor number Ta. Such a location requires dβ dz (z c/a ) < 0, (2.23) and can be obtained in the case of fluid extraction and cross-flow reversal.

  2.24)computed from the resolution of the local stability problem (B 2) for the most unstable absolute instability with β set to β(z c/a ). Following (70) in[START_REF] Pier | Bifurcation to fully nonlinear synchronized structures in slowly varying media[END_REF], matching the linear perturbation upstream of the front to the nonlinear perturbation downstream of the front introduces a correction to the leading-order frequency (2.24):

  .26) Again, though the (z c/a )-dependences are dropped for convenience in (2.25) and (2.26), all quantities inferred from the local stability problem and β(z) are computed at z c/a . The axial evolutions of the wavevector and amplitude of the global mode are now governed by the fully nonlinear stability problem, which has not been obtained so far.

FIGURE 4 .

 4 FIGURE 4. Sketch of the numerical domain.

  ) in[START_REF] Tilton | A 3D pseudospectral algorithm for fluid flows with permeable walls. Application to filtration[END_REF]. For figure15(b), the controlled perturbation is continuously applied as the sum of 50 random frequencies (see (55) in[START_REF] Tilton | A 3D pseudospectral algorithm for fluid flows with permeable walls. Application to filtration[END_REF]). Simulations used a numerical domain decomposed into two central subdomains of length L central = 2 × 60 (L central = 2 × 65 for cross-flow reversal in § 3.2.3), bounded by an inlet buffer L in = 10 and an outlet buffer L out = 15. All simulations are integrated in time until a saturated state is achieved.

  FIGURE 5. (Colour online) Numerical simulations for η = 0.85, Ta = 120 with fluid injection. (a) Isosurface where the total azimuthal velocity is equal to one half of the velocity at the rotating inner cylinder, at final time. (b) Spatio-temporal diagram of the radial velocity at mid-gap U(r mid , z, t). (c) Axial Reynolds number β(z): numerical (dashed light grey, green online, curve) and analytical (solid dark grey, blue online, curve), this latter reconstructed from expression (2.5) and numerical values of the mean axial flow at the boundaries between the inlet buffer and central domain and the central domain and outlet buffer. These boundaries are shown in solid black circles (in panel a) and vertical lines. The dash-dotted dark grey (blue online) circle (in panel a) and vertical line locate the absolute/convective boundary; and the dotted dark grey (blue online) circle (in panel a) and vertical line show the convectively unstable/stable boundary. Stable regions are labelled S, convectively unstable regions are labelled CU and absolutely unstable regions are labelled AU.

  FIGURE 6. (Colour online) Numerical and analytical results for η = 0.85, Ta = 114, with fluid injection. (a) Solid piecewise horizontal (red online) curve: dominant frequency ω max (z), exhibiting three steps at ω = 7.50, 15.01 and 30.02 (with spectral resolution ω fft = 0.17) to be compared with the dashed dark grey (blue online) horizontal line: ω glob crit = 7.40. The dotted dark grey (blue online) curve in panel (a) shows the local frequency ω loc (z) at Ta = 114. (b) Radial velocity at mid-gap U(r mid , z), at final time (solid black curve), enveloped by the amplitude U mid (z) computed from the temporal Hilbert transform (solid, red online, curve). The solid black, dash-dotted dark grey (blue online) and dotted dark grey (blue online) vertical lines are as in figure 5.

  FIGURE 7. (Colour online) Similar to figure 6, with Ta = 130. The dominant frequency ω max (z) (solid, red online, piecewise horizontal curve in panel a) exhibits four steps at ω = 7.55, 15.26, 30.51 and 61.03 (with ω fft = 0.21), to be compared with the dashed dark grey (blue online) line: ω glob crit = 7.40.

  FIGURE 8. (Colour online) Similar to figure 5 (η = 0.85 and Ta = 120), but for fluid extraction.

  FIGURE 9. (Colour online) Similar to figure 6, with η = 0.85, Ta = 113 and fluid extraction. The dominant frequency ω max (z) (solid piecewise horizontal, red online, curve in panel a) exhibits a single step at ω = 7.69 (with spectral resolution ω fft = 0.19), to be compared with the dashed dark grey (blue online) horizontal line: ω glob crit = 7.55 and the dash-dotted light grey (green online) horizontal line: ω glob front = 7.17. The solid black dotted and dash-dotted dark grey (blue online) vertical lines are as in figure 8. The dash-dotted light grey (green online) vertical line in panel (a) shows the modified front location z glob front .

  FIGURE 11. (Colour online) As in figure 5 (η = 0.85, Ta = 120), for cross-flow reversal. The dashed dark grey, blue online, circle (in panel a) and vertical line locate the double turning point where the axial Reynolds number β(z) is minimum.

  FIGURE 12. (Colour online) Similar to figure 6, with η = 0.85, Ta = 112 and cross-flow reversal. The dominant frequency ω max (z) (solid piecewise horizontal, red online, curve in panel a) exhibits two steps at ω = 5.98 and 11.96 (with spectral resolution ω fft = 0.24), to be compared with the dashed dark grey (blue online) horizontal line: ω glob crit = 6.04 and the dash-dotted horizontal light grey (green online) line: ω glob front = 6.00. The solid black, dashed, dash-dotted and dotted dark grey (blue online) vertical lines are as in figure 11. The dash-dotted light grey (green online) vertical line in panel (a) shows the modified front location z glob front .

  .22) slightly shift these values to Ta glob crit = 112.7, above which a global mode with frequency ω glob crit = 7.40 becomes unstable. The corrected critical Taylor number Ta glob crit is in agreement with the DNS, which shows no instability at Ta = 112 and an instability at Ta = 114.

  Fluid extraction: maximum instability at the outlet Next, we consider systems in which fluid is extracted along the full length of the membrane, as in figures 2(b) and 8(c). Fluid extraction causes β(z) and the local stability of the flow to decrease monotonically as fluid moves downstream. In the DNS, we impose fluid extraction by applying β inlet = 40 and β outlet = 2, for a permeability σ = 5 × 10 -5 and L central = 120. Simulations are repeated for Ta = 111,

  (b) shows that the steady wavepacket includes two subpackets, with different temporal frequencies.For fluid extraction cases, our global mode analysis predicts a linear selection criterion based on the simple turning point at z outlet , described in § 2.2.2. For the values η = 0.85 and β outlet = 1.90, condition (2.18) predicts a critical value Ta loc crit = 111.4, above which a global mode becomes unstable. Equation (2.19) predicts the critical frequency ω loc crit = 7.04. Considering the corrections at the next order, expressions (2.21) and (2.22) slightly shift these values to Ta glob crit = 112.3, above which a global mode with frequency ω glob crit = 7.55 becomes unstable. This linear critical condition Ta glob crit is in agreement with DNS, which shows no instability at Ta = 111 and an instability at Ta = 113.

FIGURE 14

 14 FIGURE 14. (Colour online) Frequency spectrum of the radial velocity at mid-gap obtained for cross-flow reversal at η = 0.85 and Ta = 120, as in figure 11, as a function of the axial coordinate z. The dash-dotted (green online) vertical line shows the location of the front z glob front . The dashed (blue online) vertical line shows the location of the double turning point. The dotted (blue online) curve the local most amplified frequency ω loc (z).

  FIGURE 15. (Colour online) Numerical simulations for η = 0.85, Ta = 130 with fluid extraction and the inlet locally convectively unstable. Radial velocity at mid-gap U(r mid , z), at final time (solid black curve) and amplitude computed from the temporal Hilbert transform (solid red curve). The perturbation at the inlet is at initial time (a) or a continuous random noise (b).

  Figure 3(b) suggests that convective and absolute frequencies weakly differ. Comparing the value of ω(z) to the local

FIGURE 16 .

 16 FIGURE 16. (Colour online) Numerical simulation for η = 0.85, Ta = 160 with fluid injection: isosurface where the total azimuthal velocity is equal to one half of the velocity at the rotating inner cylinder, at final time. The whole domain is unstable and the dash-dotted (blue online) circle shows the absolute/convective boundary.

  decayand vanish. In the case of fluid extraction, assuming that the front-forming mechanism selects the leading subpacket, the vortices develop between the upstream boundary between the locally convectively unstable and absolutely unstable regions (the dark grey, blue online, dashed vertical lines between the

  and the downstream boundary between the locally convectively unstable and stable regions (the dark grey, blue online, dashed vertical lines between the

  FIGURE 17. (Colour online) Local linear stability results, as functions of the radius ratio η and axial Reynolds number β. (a) Convective critical Taylor number Ta conv crit . (b) Absolute critical Taylor number Ta abs crit . The colourscale of the surfaces corresponds to the critical azimuthal wavenumbers n conv crit and n abs crit . Solid (blue online) curve in (a) and dashed (blue online) curve in (b) show the specific results for convective and absolute, respectively, critical conditions at η = 0.85, as in figure 3.

  This condition is achieved in the case of cross-flow reversal, as illustrated in figure2(c), where z d.t.p. is indicated by a circle. At leading order, the global mode frequency imposed by a double turning point is ω loc = ω abs (z d.t.p. ), computed by solving the local stability problem (B 2) for the most unstable absolute instability, with β set to β(z d.t.p. ). A global mode governed by a double turning point should thus be observed as the Taylor number overcomes the critical absolute Taylor number at the double turning point:

	Ta loc crit = Ta abs crit (z d.t.p.
	0 and	∂ 2 Im(ω abs ) ∂z 2	(z d.t.p. ) < 0,	(2.11a,b)

in addition to condition (2.10). As the z-dependence of the dispersion relation amounts to the axial variations of the axial Reynolds number β, condition (2.11) requires: dβ dz (z d.t.p. ) = 0 and d 2 β dz 2 (z d.t.p. ) > 0. (2.12a,b)

  2.1 and fluid extraction shown in § 3.2.2, or a double turning point, as developed in § 2.2.1, in the case of cross-flow reversal shown in § 3.2.3. The agreement is already satisfying with the leading order of the critical Taylor number Ta loc crit as given by (2.13) and (2.18) and further improved when corrected to Ta glob crit as given by (2.16) and (2.21). Concerning the frequency selection, the instabilities observed in the DNS are almost always composed of stacks of several subpackets of waves oscillating at different frequencies. Nevertheless, the leading subpackets always exhibit a frequency related to a global mode selection criterion. The linear selection by a simple turning point at the inlet (2.19) works accurately in the case of fluid injection, even when the Taylor number substantially departs from its critical value, and its dynamics downstream are clearly nonlinear, and without requiring to proceed to the next order (2.22). The linear selection by a simple turning point at the outlet (2.19) works accurately in the case of fluid extraction for Taylor numbers close to the critical conditions, as in figure 9, where, moreover, the instability is fully synchronized. The linear selection by a double turning point located where the axial Reynolds number β(z) is minimum (2.14) works accurately in the case of cross-flow reversal for Taylor numbers close to the critical conditions, as in figure 12. For fluid extraction and cross-flow reversal, proceeding to the next orders (2.22) and (2.17) slightly improves those results.

Appendix A. Axially varying base flow

As explained in [START_REF] Tilton | Pressure-driven radial flow in a Taylor-Couette cell[END_REF], a base flow for the axially varying configuration obtained when the rotating inner cylinder is permeable while the stationary outer one is impermeable is obtained by an asymptotic expansion. This expansion requires one to assume that both the non-dimensional permeability √ σ and the pressure drop along the axial direction P/ z, primarily related to the axial flow, are small. The axially varying base flow is then obtained as an asymptotic expansion up to the second order in ε = √ σ :

(A 1a-c) with the varying axial Reynolds number β(z) as given in (2.5). According to the asymptotic expansion, the axial Reynolds number β(z) is first order in ε, so that U b (r, z) is second order in ε, the first term in V b (r, z) is zeroth order while the second one is second order, and W b (r, z) is first order.

The radial shape functions in (A 1) read as

where

and

Since the pressure field of the base flow is not required for the stability analysis, it is not included here.

Appendix B. Local convective/absolute stability analysis

The local linear convective/absolute stability analysis is similar to that in [START_REF] Martinand | Absolute and convective instability of cylindrical Couette flow with axial and radial flows[END_REF], and is only summarized here. The velocity and pressure fields are decomposed as

, where V inv b (r) = (0, Tav b,0 (r), βw b (r)) t is the axially invariant TCP base flow, quantified by Taylor number Ta and constant axial Reynolds number β, and [V inv p (r, θ, z, t), P inv p (r, θ , z, t)] is a small time-dependent perturbation. Equations (2.3) and (2.4) are linearized about the base flow V inv b , and the perturbation is sought as

where ω is the temporal frequency, k the axial wavenumber, and n the azimuthal wavenumber. This yields an eigensystem of r-differential equations satisfied by the radial shape functions [v p (r), p p (r)], written in condensed form as:

for which [v p (r), p p (r)] is an eigenvector and ω the associated eigenvalue. Operators A and B can be obtained in appendix B of [START_REF] Martinand | Weakly nonlinear analysis of cylindrical Couette flow with axial and radial flows[END_REF]. System (B 2) is solved for a given set of parameters (Ta, β, η, k, n) using a Chebyshev spectral method. The solvability conditions of the set of systems of differential equations obtained by differentiating (B 2) with respect to k, n and the parameters Ta and β lead to the several partial derivatives of the complex frequency ω required by the convective/absolute stability analysis and the global mode selection criteria (besides the axial group velocity

n ω and ∂ k ∂ n ω are also used in the analysis), and the corresponding partial derivatives of the eigenvector. As detailed in [START_REF] Martinand | Absolute and convective instability of cylindrical Couette flow with axial and radial flows[END_REF], critical conditions (k conv crit , n conv crit , Ta conv crit ) for the most unstable convective mode ([v conv crit , p conv crit ], ω conv crit ) are found by imposing Im(ω) = 0 and Im(∂ k ω) = 0, with k real and n integer. Critical conditions (k abs crit , n abs crit , Ta abs crit ) for the most unstable absolute mode ([v abs crit , p abs crit ], ω abs crit ) are found by imposing Im(ω) = 0 and ∂ k ω = 0, with k complex and n integer. Those critical conditions and modes are found as functions of the other parameters β and η, as shown in figures 3 and 17.