Laminar–Turbulent Intermittency in Annular Couette–Poiseuille Flow: Whether a Puff Splits or Not

Direct numerical simulations were carried out with an emphasis on the intermittency and localized turbulence structure occurring within the subcritical transitional regime of a concentric annular Couette–Poiseuille flow. In the annular system, the ratio of the inner to outer cylinder radius is an important geometrical parameter affecting the large-scale nature of the intermittency. We chose a low radius ratio of 0.1 and imposed a constant pressure gradient providing practically zero shear on the inner cylinder such that the base flow was approximated to that of a circular pipe flow. Localized turbulent puffs, that is, axial uni-directional intermittencies similar to those observed in the transitional circular pipe flow, were observed in the annular Couette–Poiseuille flow. Puff splitting events were clearly observed rather far from the global critical Reynolds number, near which given puffs survived without a splitting event throughout the observation period, which was as long as 104 outer time units. The characterization as a directed-percolation universal class was also discussed.


Introduction
The discontinuous reverse transition of wall-bounded turbulence into a laminar flow is a fundamental problem that has been studied for many years, while the laminar-to-turbulent transition is rather smooth, or its critical point is often well predicted by linear stability theory. Subcritical flows in the reverse transition are known to feature two regimes in competition, namely, laminar and turbulent, in which there occurs large-scale intermittency that coexists spatially with a laminar flow. The large-scale nature of localized turbulence often forms a regular pattern once established. The intermittent structure or formation pattern of localized turbulence varies depending on the flow system, and a number of studies have been conducted on canonical flows, such as a circular pipe flow (CPF) and planar flows. In the CPF, a so-called equilibrium turbulent puff, or simply a "puff," is localized in the streamwise direction, resulting in uni-directional intermittency. The puff turbulence is sustained within a Reynolds-number range based on the bulk velocity U and the pipe diameter D of Re D = 2000-2700 [1]: Although there are some differences depending on the experimental conditions, such as the disturbance introduction method and the pipe length [2][3][4], studies have indicated that the puff's nature is deeply related to the determination of the lower-limit Reynolds number (the global critical Reynolds number, Re g ), above which turbulent motions can survive globally. Streamwise-localized solutions underlying the puff have been found, and Hopf bifurcations to new branches including unstable periodic orbits are expected to cover the turbulent attractor [5]. It is also known that puffs can split (or proliferate) more frequently than their decay and have a finite lifetime even at Re > Re g [6][7][8][9]. Avila et al. [10] identified Re g = 2040 ± 10 for the CPF by monitoring both the puff-splitting time and the decay time. Recent attempts have been made to elucidate the puff-driving flow in a plane Couette-Poiseuille flow (pCPf), allowing localized turbulence to be tracked for long periods of time while stationary in the observation window. A recent quench experiment on the decay of Couette-Poiseuille turbulence is likely to approach the crossover of the decay rate, that is, the quantitative identification of Re g [50]. However, to the best of our knowledge, except for in limited studies [51,52], there is no DNS available for the subcritical transition process of the aCPf. The present DNS is the first to explore laminar-turbulent intermittency in a low-η aCPf.
The remainder of this paper is organized as follows. Section 2 presents the flow configuration, dimensionless parameters, and equations used in our simulations. In Section 3, which is dedicated to the preliminary results, we validated the current code and illustrated the parameter dependence of the base flow in terms of the mean friction on the inner cylinder. Section 4 begins with a puff characterization of the observed turbulent patches. Space-time diagrams of a turbulent quantity revealing the puff splitting and decay are then presented. All results are summarized and discussed in Section 5.

Problem Setup and Methods
The problem under consideration is the turbulent annular flow of an incompressible Newtonian fluid, for which the governing equations are the equation of continuity and the Navier-Stokes equation, as described in the classical cylindrical coordinate system of (x, r, θ): Here, the velocity vector is represented by u, or (u x , u r , u θ ), which are the respective components in (x, r, θ); p is the pressure, and t is the time. These quantities are non-dimensionalized and are marked by a * superscript: u * = u/u w , p * = p/ρu 2 w (ρ, density), t * = tu w /h, x * = x/h, and r * = r/h, where u w and h are the inner-cylinder axial velocity and the gap between the two cylinder radii, respectively, as illustrated in Figure 1. The Reynolds number Re w is therefore based on u w , h, ρ, and the fluid viscosity µ, whereas another definition using one-quarter of Re w is more conventional for studies on the pCf [15,16,22]. In only the axial-direction component of Equation (2), a constant pressure gradient in x is added as an external force term, −dP/dx, with the axial unit vector e x . In addition to the imposed pressure gradient, the flow is driven by an axial translation of the inner rod with a constant velocity of u w > 0. The x-axis corresponds to the central axis common to both cylinders, and the radius ratio of η is an important geometrical parameter and is set to 0.1 for the main analysis in this study. Periodic boundary conditions are imposed in both the x and θ directions, and no-slip boundary conditions are enforced at the wall surfaces of the cylinders. In the following sections, the imposed pressure gradient is re-defined as the pressure gradient function F(p), which is normalized as and can be interpreted as the ratio of the imposed pressure gradient (i.e., the Poiseuille-like driving force) against the wall-bounded viscous shear stress (the Couette-like driving force). As introduced above, there are two control parameters for the flow under consideration, i.e., Re w and F(p). Poiseuille-like flows are realized for a large F(p), whereas Couette-like flows are obtained for a small F(p), and a specifically pure Couette flow corresponds to F(p) = 0. As indicated in [44], the ratio of the shear stress at the two walls, which can be defined by γ = τ in /τ out in an aCPf, is another candidate of the control parameter relevant to a Couette-Poiseuille flow. Flows with γ ≈ 0, or a shear-less inner cylinder wall, are of special interest because they exhibit nearly zero mean shear at the moving rod, and can thus be a model for an understanding of the puff dynamics in a pure pipe. Under such conditions, the inner cylinder practically affects the core flow only as an impermeable thin rod, and the coherent turbulent structures and turbulent production that dominantly occur near the static outer cylinder wall mimic those found in a canonical system of the CPF. Although the system chosen here is closer to a CPF than to an aPf or an aCf, it should be noted that the different boundary conditions regarding the inner rod preclude a mathematical homotopy continuation with the CPF. Except for a fully laminar flow state, F(p) providing γ = 0 is not explicitly obvious, and thus, a parametric survey must be conducted for each given Re w . In this study, we conducted a preliminary survey of the F(p) dependence of τ in for several Re w values using a DNS with a medium-scale computational domain, as reported in Section 3. Based on these results, we selected F(p), which will provide τ in ≈ 0 (γ ≈ 0) at each tested Re w , and accordingly applied the main DNS using a large-scale domain to reduce the spatial limitation on the laminar-turbulent coexistence. Figure 1. Couette-Poiseuille flow in an annular channel between two concentric cylinders with a radius ratio of η ≡ r in /r out = 0.1, driven by a constant pressure gradient and an inner-cylinder axial movement. In this study, the pressure gradient is adjusted such that the mean velocity gradient on the inner-cylinder surface is approximately zero; that is, The numerical conditions of the preliminary and main simulations for η = 0.1 are summarized in Tables 1 and 2, respectively. Long domain sizes of 51.2h and 409.6h were employed in the axial direction to capture a single turbulent puff and expected multiple puffs, whereas the radial and azimuthal domain lengths were of geometrically determined values of h and 2π, respectively. The grid resolutions have been confirmed to be fine such that fine-scale eddies in turbulent patches are well resolved, at least for the particularly interesting transitional regime of Re w ≤ 1600.
Equations (1) and (2) were discretized using a staggered central finite-difference method, where the fourth-order central difference scheme was used in both x and θ, along with the second-order scheme in r on a non-uniform radial grid. A time advancement was performed using a fractional-step second-order Adams-Bashforth scheme in combination with a Crank-Nicolson scheme for the radial viscous term. The Courant-Friedrichs-Lewy (CFL) condition was continuously monitored in all directions, and accordingly, the time-step ∆t constraint for the nonlinear terms was enforced to ensure stability. The details of the numerical method were reported in the literature [34,53]. The code validation carried out is discussed in the next section.

Preliminary Simulations
The reliability of the current simulation code may be demonstrated through a comparison with the existing pCPf DNS database at a comparable Reynolds number. Kasagi and coworkers [43,54] applied a DNS of several pCPfs and released their database obtained, from which a condition of (Re w , (p)) = (6000, 15.96) was chosen for the code validation during this test. At this Reynolds number and the mean pressure gradient, the pCPf is under a fully turbulent state throughout the channel, and no large-scale intermittency occurs. Its friction Reynolds number Re τ , normalized by the friction velocity on the fixed wall and the half width of the gap, is 154. Kuroda et al. [43] adopted a spectral method with a 128 × 128 Fourier series in the horizontal directions and Chebyshev polynomials up to order 96 in the wall-normal direction. Their domain size was 2.5πh × h × πh, whereas our counterpart simulation on an aCPf of η = 0.9 employed a nearly equal domain size of 8h × h × π/8 (≈ 3h at the gap center) in (x, r, θ). Figure 2 shows comparisons of the present mean and second-order statistics. An overbar, such as u x (r), denotes an ensemble-averaged quantity with respect to t, x, and θ, and subscript 'rms' indicates a root-mean-square value. The present control parameters of Re w and F(p) for our aCPf are 6000 and 16.0, respectively, and the resulting friction velocity and the friction Reynolds number on the fixed outer cylinder wall are 0.050u w and 151. The present results shown in Figure 2 are in reasonable agreement with the reference study, despite the wall curvature of the aCPf. A noticeable difference is detected only near the fixed wall (y/h ≈ 0.9), where the profile of u x exhibits a steep gradient, and thus, those of the streamwise turbulent intensity u xrms and Reynolds shear stress u x u r have peaks. The rather coarse grid resolution and the low-order spatial discretization (our finite difference code versus the previous spectral code) might affect the accuracy of the present simulation. In addition to the peak values of u xrms and u x u r , the second-order statistics from the present DNS and those of Kuroda et al. [43] agree well, particularly considering the differences in the flow geometry. The current Fortran code has been employed in different studies for several different boundary conditions [34,35,37,38,55], and thus, no further validation will be shown here.  [43] for plane Couette-Poiseuille flow (pCPf), respectively. Here, the wall-normal coordinate y represents the distance from the inner (bottom) wall, that is, y = r − r in .
In this study, we simulated a low-η annular flow that mimics a CPF by approximating the base flow, or the mean velocity profile, to that in the CPF. In the CPF, the velocity profile reaches its maximum at the pipe center; the velocity gradient becomes zero at the pipe center, and is at maximum on the surface of the pipe (i.e., the outer-cylinder surface). To match the base-flow characteristics of a CPF in an annular system, it is necessary to conduct a parametric investigation on the appropriate magnitude of the pressure gradient applied in the annular channel. As a preliminary analysis, we employed a medium-scale computational domain to reduce the computational cost of the parametric study. The computational domain size is smaller in the x direction than the present main analysis shown in Section 4. The streamwise length of the domain, L x , was sufficient to capture one turbulent puff. The purpose of the preliminary analysis is to identify the value of the pressure gradient function F(p) at each Reynolds number such that the friction coefficient on the inner cylindrical wall, C f ,in , is practically zero, where C f ,in is defined by the following: where U is the bulk mean velocity obtained through a simulation. The positive/negative sign of C f ,in corresponds to the positive/negative velocity gradient on the wall surface of the inner cylinder. Given du/dr < 0, C f ,in < 0, and vice versa. Table 1 shows the calculation conditions and the ranges of Re w and F(p) in the preliminary analysis. In the preliminary analysis, the calculation area in the mainstream direction was set to a smaller calculation area than that for the main analysis, but can capture one turbulent puff. Figure 3 shows the F(p) dependence of C f ,in at several values of Re w near the global critical value. In each analysis plotted in the figure, a turbulent field with a high Reynolds number at equilibrium for each given F(p) was set as the initial flow field, and the ensemble-averaged C f ,in value was acquired after reaching a statistically steady state. Note that laminarization did not occur in any of the cases shown here. In general, as F(p) increases, C f ,in increases monotonically while changing from a negative to a positive value. This is consistent with the transition of the mean velocity profile from Couette-like to Poiseuille-like, and it can be confirmed that "the turning point" of F(p) indicating C f ,in = 0 increases with Re w . According to this Reynolds-number dependence, an extrapolation predicts a value of F(p) that brings C f ,in = 0 at a lower Re w , by which the main DNS analysis in the next section was executed.

Results and Discussion
The main DNS at Re w ≤ 1600 for which a laminar-turbulent intermittency was clearly confirmed through the preliminary analysis is presented in this section, and the characteristics of the localized turbulence are discussed. Table 2 summarizes the numerical conditions, including the friction Reynolds number Re τ that was obtained. As in the preliminary analysis, the radius ratio was η = 0.1, and the computational domain was extended only in x; however, the grid resolution was not changed. Because Re τ is lower than that of the preliminary analysis, the grid spacing in terms of the wall units was finer. The table also shows the grid resolutions based on the friction velocity u τ under each condition: where u τ,in and u τ,out are defined by the corresponding wall shear stress, τ in , and τ out , as well as by the relation τ = ρu 2 τ , from which inner units can be defined. For a low Reynolds-number regime of Re w < 1600, which is of interest in this study, the grid spacings of ∆x + < 8, ∆r + min < 0.2, ∆r + max < 3, and ∆z + < 4 are comparable to or higher in resolution than those in previous studies [35,37,53]. The initial conditions during each analysis adopted a turbulent field with a one-step-higher Reynolds number, but reduced the Reynolds number adiabatically. In other words, the study was carried out carefully such that the sudden drop in the Reynolds number will not be a proximate cause of laminarization. Table 2. Numerical conditions for the main DNS with a long domain. The grid resolutions of (∆x, ∆r, ∆θ) are described in their dimensionless form based on u τ and µ/ρ. The minimum and maximum ∆r of the radial direction, in which we used non-uniform grids, are shown. † Laminar values from a laminarized case.  Figure 4 presents a three-dimensional visualization of localized turbulence in the form of puffs, which is observed as an equilibrium state reached after a lengthy simulation under the condition of Re w = 1600 and F(p) = 6.5. The turbulent region can be clearly detected by showing the radial velocity fluctuations or the wall-normal velocity component. The threshold value of ±0.03u w for the iso-surface was arbitrarily chosen to extract its typical arrowhead shape similar to that of a puff. A slight change in this threshold value does not significantly affect the interpretation of the present results. In the snapshot, multiple turbulent patches, called 'puffs' hereafter, can be confirmed to be distributed intermittently with respect to the streamwise direction. The blank regions between neighboring puffs can be regarded as being in a laminar flow because of an insignificant fluctuating velocity, implying the well-established coexistence of laminar and turbulent regions in the aCPf. As is clear from the enlarged figure, the puff has an arrowhead shape, and the puff extends downstream in the center of the outer pipe. Although the average velocity gradient on the inner cylindrical wall surface is almost zero, this situation is considered to be due to the similar driving mechanism of the puff of the CPF. For the CPF, Shimizu et al. [11] reported that turbulence in the puff originates from low-speed streaks, as well as from streamwise vortices along the (outer-)pipe wall and across the trailing edge of the puff through the Kelvin-Helmholtz instability, which induces velocity fluctuations that propagate downstream faster than the puff itself in the core region. Such a driving mechanism of the puff is also common to the present aCPf with nearly zero C f ,in . The streamwise size of each puff is approximately 30 times the gap width h, which corresponds to 15 times the hydraulic diameter, and is consistent with that of the puff observed in the CPF [1,7,8,12]. The array of puffs seems variable in intervals, but is likely not less than 30h. The wavelength and periodicity of the puffs are examined using two-point correlation functions of the turbulence quantities. In the x and θ directions, the auto-correlation coefficients are defined as follows:

Puffs in Annular-Pipe Flow
and where i ∈ (x, r, θ). Figure 5 shows the two-point correlation coefficients of each velocity component for the case visualized in Figure 4. The statistical dataset was accumulated over the time of 5000h/u w after achieving a pseudo-equilibrium state of multiple puffs. From Figure 5a, the axial periodicity and interval of the puff can be estimated. First, we note that the three curves at different y * exhibit consistency, implying that flow state and patterning are only weakly dependent on y or r. As also plotted in (b) and (c) for the other directional components, fine-scale turbulent structures inside a puff should have a rather short streamwise extent, and indeed, the profiles of R rr and R θθ fall to almost zero at ∆x < 5h. The profile of R xx also decreases drastically for a small ∆x, although its significant oscillation for a long axial extent suggests a spatial coexistence of laminar and turbulent regions rather than turbulent structures, since these two flow states have different mean velocity profiles, particularly near the walls. The oscillations observed in Figure 5a are somewhat strong at both the inner and outer walls, relative to the gap center. The profile of R xx takes the first negative local minimum at ∆ ≈ 30h and shows regular spikes at intervals of approximately 60h. The correlation is not zero even at half the computational domain length (L x /2 = 204.8h). Peaks at 60h, 120h, and 180h manifest the presence of seven distinct puffs in L x on average. This suggests that the puffs at this Reynolds number tend to be arranged regularly throughout the axial extent. If the puff spacing is irregular, the correlation coefficient distribution should not show periodic fluctuations and should asymptotically approach zero. This regularity of the puff arrangement may differ from the characteristic of the DP universal class, which should exhibit a wide-scale invariant pattern close to the critical point [39,41]. Mukund and Hof [3] reported a similar aspect on multiple puffs in a CPF, where they referred to the wave-like fashion as 'puff clustering'; that is, the resultant pattern of clustering puffs was observed to propagate like waves. They also pointed out that interactions between puffs were responsible for the approach to the statistical steady state and strongly affected the percolation threshold. This may predict a difference in the global stability between a single puff (i.e., isolated puffs) and multiple puffs (puff clustering), as discussed in Section 4.2.  Figure 5d-f indicate the azimuthal intervals between fine-scale turbulent structures, such as low-speed streaks inside the puff. There exists no large-scale pattern in the azimuthal direction, unlike those of the helically shaped turbulent patches in high-η aPf [35] and aCf [37]. The blue curve in Figure 5d, measured near the outer cylinder wall, only has a peak at ∆θ = π/2. The cross-sectional flow pattern observed here consists of four low-speed streaks close to the outer wall spaced at π/2. This azimuthal configuration regarding turbulence inside the puff is in agreement with those found in the CPF [6].

The azimuthal two-point correlation functions shown in
The presence of turbulent equilibrium puffs was observed even at Re w < 1600, and the flow field finally reached the fully laminar state at Re w = 1500. Although the space-time diagram (STD) and the turbulent fraction, Ft(t) (the plot of which is shown later), reveal a tendency toward laminarization at Re w = 1525, one turbulent puff was maintained in the present computational domain at least during the present observation time of >1.3 × 10 4 , and the laminarization was not completed. If normalized by the hydraulic equivalent diameter 2h and bulk velocity U, the Reynolds numbers of Re w = 1600 and 1500 correspond to Re D = 2190 and 2045, respectively. This range of Re D = 2045-2190 is close to or slightly narrower than that for the counterpart of the CPF (Re D = 2000-2700 [1], 2040-2400 [10], 2300-3000 [2], and 2000-2200 [4]). In particular, a discrepancy in the lower bound value of the subcritical transition regime, that is, the global critical Reynolds number, is of interest, although the similarity with the results by Avila et al. [10] is rather surprising. A cause of this discrepancy remains unclear: One of the main causes may be the presence of the inner cylinder, which suppresses turbulent motions across the central axis in the case of an aCPf. Another cause may be the non-slip inner-cylinder surface, which prevents a puff from splitting into two puffs in the case of an aCPf. Shimizu et al. [8] proposed a model process of puff splitting in the CPF, which starts with an azimuthally isolated streak propagating downstream through the laminar-turbulent interface of the puff. An emitted streaky disturbance can be a seed of a "daughter puff," which spreads again in the azimuthal direction and grows into a turbulent puff after leaving the parent puff sufficiently far away. As a system even closer to the CPF, an ideal aPf with a stress-free boundary condition at the inner wall can be analyzed, although such an unpractical situation will be considered as a future task. In terms of the conjecture that puff splitting is unlikely in the aCPf relative to the CPF, we traced puffs with lengthy simulations, and their STDs are shown in Section 4.2.

Space-Time Diagrams
As for the puff turbulence in CPF, it is well known that a turbulent puff can split into two puffs over time, the turbulence between puffs should attenuate and become a laminar pocket, and one or both puff(s) should decay quasi-stochastically because of their finite lifetime [3,6,8,9]. Avila et al. [10] observed the puff-turbulence sustainment only due to puff-splitting events that have time scale shorter than the puff-decay time scale. These features may be identified from the temporal development of the puff spatial distribution. The STDs of the present aCPf are shown in Figures 6-8, where the horizontal axis is the streamwise coordinate in a frame of reference moving at a certain velocity, and the vertical axis represents the dimensionless time at each Reynolds number. The frame-moving velocity is nearly the mean gap-center velocity, which also corresponds to the propagation velocity of an observed single puff. The color contour shows the azimuthal average of the radial velocity at mid-gap, u r θ = 2π 0 u r (x, h/2, θ, t)dθ/2π, such that the laminar and turbulent regions can be clearly distinguished. Although the apparent length of each turbulent puff depends on the criterion used to discriminate it from the surrounding laminar flow, a different choice does not change the qualitative conclusions obtained.  Figure 6. Space-time diagram for (a) Case 1 at Re w = 1600 and F(p) = 6.5, (b) Case 2 at Re w = 1600 and F(p) = 6.5 with a different initial condition from that in Case 1, and (c) a typical discrete expansion process of turbulence in a subcritical transitional pipe flow at Re = 2300, cited from Avila et al. [10]. In (c), the contour color indicates the cross-sectional average of the streamwise vorticity squared, where red and blue correspond to turbulent puff and laminar regions, respectively, and the Reynolds number Re is based on the mean velocity U and the pipe diameter D. In (a,b), the contour shows the azimuthally averaged radial velocity u r θ at the gap center. The axial distribution is monitored from a moving frame of reference with a speed close to the puff propagation. The temporal development is monitored from t = 0, that is, the beginning of each DNS with a higher-Re w field with more puffs; therefore, some initial puffs decayed immediately after the start of the simulation. Not to scale (aspect ratio x:y = 10:1).
We first present the results for Re w = 1600 and F(p) = 6.5, as discussed in Section 4.1. The flow field visualized in Figure 4 was first achieved through an adiabatic decrease in Re w (with a change in F(p), accordingly) from a fully turbulent regime, and was then used as the initial condition for the following simulation to trace the behavior of the puff in the phase diagram of the L x -space and time for as long as possible. The STD obtained is shown in Figure 6b, which monitors the pattern starting from an initial state with several puffs-the isolated turbulent patch featured as a red and blue segment at given time t. The overall puff pattern remains intrinsically spatiotemporally intermittent and exhibits both puff decay and splitting very frequently. These individual puffs have statistically well-defined lengths, similar to those in a CPF [8]. The number of puffs captured in the present domain is roughly constant between 5 and 7, and it is again confirmed that the puff intervals tend to be constant even if splitting or attenuation occurs in each individual puff. For this reason, Figure 5a reveals regular oscillations in the correlation function R xx (∆x), while a snapshot visualized in Figure 4 happens to have no periodicity in the puffs when considering the complete pipe length. Figure 6b may invoke an STD obtained from experimental and numerical observations of a DP-like feature in other flow systems [39,42]. Another DNS labeled as Case 1 was repeated for the same parameter set of (Re w , F(p)), but with a different initial condition with a single puff, which was prepared from a lower-Re w DNS (see Figure 6a). The initial single puff is sustained for a long period of >5000h/u w , during which it splits irregularly, the first time at tu w /h ≈ 3000 and the second time at tu w /h ≈ 4000, but both newborn puffs decay after they are separated from their parents. A newly emitted daughter puff by the third splitting tu w /h ≈ 5500 grows and successively produces grandchild puffs. In addition, there are many signs of puff splitting. The puff turbulence eventually covers the entire domain, yet is intrinsically patchy, as in Case 2. It can be concluded that, in an aCPf similar to a CPF, the turbulent puff can split, regardless of the initial field, in qualitative agreement with a typical STD sample of a CPF [10], as displayed in Figure 6c. Re w = 1500 and F(p) = 5.3 (Re D ≈ 2045). The contour shows u r θ at the gap center. The axial distribution was monitored from a moving frame of reference. The same initial condition was applied for all cases presented here. Not to scale (aspect ratio x:y = 10:1). Figure 7a shows the STD of Re w = 1600 and F(p) = 6.5 (Case 2), but the speed of the moving frame of reference is modified such that puffs appear to be stationary with respect to space. With this adjustment, the propagation speed of the puff can be estimated as approximately 0.625u w . According to Figure 6a, when tracking a single puff in Case 1, the propagation speed is slightly faster and ≈0.65u w . The result is reasonable because the bulk velocity generally decreases with the expanding turbulent region. Figure 7b is an STD at Re w = 1500 with the same horizontal coordinate of (x − 0.625u w )/h, showing the eventual return to laminar flow. Once a puff starts to decay, its turbulent patch seems to accelerate slightly and takes approximately 300u w /h to attenuate completely. Before that, it took more than 4200h/u w before the system settled to the fully laminar state. While the flow at Re D = 2190 of Figure 7a Table 1 for each given F(p) value. The contour shows u r θ at the gap center. The axial distribution is monitored from a moving frame of reference. The same initial condition was applied for all cases presented here. Not to scale (aspect ratio x:y = 10:1).
We further investigated the intermediate range between the two above-discussed cases (2045 < Re D < 2190) to elucidate the trends in the frequency or time of the puff-splitting events. Figure 8 presents an STD at each control-parameter set. In all DNSs presented in the figure, the initial conditions are exactly the same. In the figure, six puffs can be seen initially, but two or three of them decay immediately, particularly in the lower-Reynolds-number cases. At the lowest Re w shown in Figure 8d, puffs disappear one after another on a time scale of O(1000h/u w ), and finally, one puff remains. There is no sign of decay in the surviving puff even after 13,000h/u w , but it is likely that the puffs will stochastically disappear and laminarize if a much longer simulation is available. This might also be true for the other cases presented here. At Re w < 1600, no puff splitting was observed, resulting in only puff damping. In only Figure 8b, a sign of puff splitting is detected at tu w /h ≈ 7500, although the "daughter puff" is not perfectly formed, and is finally attenuated before leaving the parent puff. Note here that a further DNS indicates no qualitative change in the flow pattern at least until tu w /h = 11,500 also for Re w = 1540, although not shown in the figure. According to a similar type of study [10], the puff splitting in the CPF was observed both numerically and experimentally for Re D > 2200, whereas clear splitting was measured in their experiments down to Re D = 2025 < Re g (=2040). If our observations were continued as long as 10 7 outer time units, as Avila et al. [10] experimentally did, the current system of the aCPf could exhibit a puff-splitting event even below the true Re g , which is not exactly determined as of now. At least, it can be said that the puff decay and splitting rates at this stage differ strongly from those observed at Re w = 1600 (Re D = 2190). As for this regime, a conclusion similar to an experimental study on a CPF [3] can be drawn, i.e., the cluster of puffs in a wave-like fashion results in fewer puff-splitting events in the STD, whose visual appearance differs from the STD for a DP universality class. Such well-organized distances between active sites (corresponding to the puffs) and the absence of splitting events are different features from those of the DP. Figure 9 shows the temporal change in the turbulent fraction, F t (t), which is the spatial ratio of the turbulent region to the entire calculated region, including both the turbulent and laminar regions. Here, F t (t) ≈ 1 indicates a fully turbulent state, and F t (t) = 0 is a fully laminar state. We set a threshold v th to distinguish between laminar and turbulent regions such that F t (t) ≈ 0.5 in the Reynolds number region where turbulent puffs densely appear in an axial extent, as in the case of Re w = 1600 and F(p) = 6.5 visualized in Figure 7a. Figure 9a shows the temporal change of F t (t) at Re w = 1500 and F(p) = 5.3, that is, the case diagnosed as a laminar regime by a visualization in Figure 7b, employing three different threshold values (v th , v th2 , and v th3 ). It can be confirmed that the time change of F t (t), particularly the gradient of the curve, does not depend on the threshold value.
When v th = 0.005, the temporal changes in F t (t) at several Reynolds numbers below Re w = 1600 are plotted in Figure 9b. In the vicinity of the critical point, a (1+1)-D DP universality class should obey a power law of F t (t) ∝ t −0.159 over time. From the figure, the current data at Re w = 1575-1550 seem to be consistent with (1+1)-DP, although more data and more exponents will be needed to properly confirm this trend. However, it should be noted that, for Re w ≤ 1550, none of the puffs split and turbulent puffs were only attenuated, as shown in Figure 8a. This result suggests that a value close to the critical exponent of DP can be obtained even under a non-DP phenomenon of a simple decaying process without splitting. We should regard this result as a 'spurious' DP feature because the puff splitting (or an active site that creates offspring) is a requisite for the critical point and, hence, DP behavior. In other words, this reminds us to take caution regarding the judgment of a DP within the laminar-turbulent intermittency.

Conclusions
We performed direct numerical simulations (DNSs) of the concentric annular Couette-Poiseuille flow (aCPf) and investigated the laminar-turbulent intermittent field of the so-called puff turbulence, particularly during its subcritical transition. From previous studies, the laminar-turbulent intermittency in annular flows (a pure Couette flow [37] or Poiseuille flow [34]) exhibits the helically shaped turbulent pattern with bi-directional spatial intermittency and puff turbulence with uni-directional intermittency, depending on the radius ratio. This fact leads to a unified understanding of the formation of localized turbulence patterns of different systems, including planar and circular pipe flows (CPFs); however, these analyses were conducted under conditions in which the basic velocity profiles do not qualitatively match those of the CPF. In this study, the radius ratio (of the inner/outer radii) was as low as 0.1, and the mean pressure gradient was imposed such that the inner-cylinder surface had a zero velocity gradient on average, so that the CPF was alternatively simulated by an annular system. Multiple puffs were demonstrated using a long computational domain in the axial direction, and the presence or absence of the puff-splitting event and its onset Reynolds number were investigated using a long-term DNS. The Reynolds number was reduced adiabatically from the fully turbulent field, and the following results were obtained.

•
At Re w = 1600, puff-splitting events occur along with stochastic puff decay, resulting in wave-like fashion of multiple puffs with constant intervals.
• At Re w < 1600, no puff-splitting event occurs, but initially given individual puffs survive over a present observation time of at least 10 4 h/u w , maintaining the intervals among the puffs.
• At Re w ≈ 1550, a 'spurious' feature of (1+1)D-DP was detected during the quenching process even without puff splitting, and a lower Re w deviates from the DP critical exponent.
• At Re w = 1500, the flow becomes fully laminar after the non-trivial finite lifetime of the puff.

•
The range of Re w = 1500-1600 (with the accordingly changed F(p)) corresponds to the bulk Reynolds-number range of Re D = 2045-2190 based on the hydraulic diameter and bulk velocity.
The question considered in this study was whether puff splitting can occur in an aCPf, which essentially has a non-slip inner cylinder. In fact, puff splitting was clearly observed at Re w = 1600, and a sign of splitting was detected at Re w = 1540, which may be close to the global critical point, Re g . This result guarantees that the planar system and the in-pipe system can be linked via the annular system. Near the criticality, oblique turbulent stripes grow or split in the longitudinal direction of the band, but the mainstream directional splitting, as seen in the CPF, is less pronounced in the planar flows. Our results suggest that the localized structures seen in both the planar and pipe flows can cause mainstream directional splitting. However, we should note that no completed puff splitting was detected near Re g . The puff splitting could be observed for Re w < 1600 and even below Re g by increasing both the observation time and domain by orders of magnitude. Such a task to explore the exact Re g value as well as the Reynolds-number dependence of the puff-splitting time scale near Re g is a challenging one that is almost impossible at present. Another possible approach is to study lifetimes of single puffs [56] and time scales of splitting [10] at conditions away from Re g , as in earlier studies on the CPF. This may allow us to discuss whether the current system behaves more like quasi-1D Couette flow or like pipe flow, as done by Shi et al. [57] for a pCf. We would like to report on this issue in another paper. Moreover, the characterization of a DP universal class remains skeptical. Similarly to the critical phenomena of the DP universal class, the region of the absorbing state (laminar-flow gap among puffs) should increase as the criticality approaches. From these facts, it is important to verify the DP feature after further expanding the axial computational domain. In addition, since the transition process of an aCPf has a dependence on the radius ratio and F(p), a parametric study will also be addressed in the future. Acknowledgments: Numerical simulations were performed on SX-ACE supercomputers at the Cybermedia Centre of Osaka University and the Cyberscience Centre of Tohoku University.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript: