2. Formulation of the problem

We consider an incompressible Newtonian fluid of kinematic viscosity $\nu$ , filling the gap between two infinitely long concentric rotating cylinders (figure 1 a). The angular velocities of the inner and outer cylinders, of radii $r^*_{i}$ and $r^*_{o}$ , are denoted as $\Omega _{i}$ and $\Omega _{o}$ , respectively. A complete set of independent dimensionless physical parameters characterising the problem are the radius ratio $\eta =r^*_{i}/r^*_{o}$ and the two Reynolds numbers $R_i={\textrm{d}}r^*_{i}\Omega _{i}/\nu$ and $R_o={\textrm{d}}r^*_{o}\Omega _{o}/\nu$ , where $d=r^*_{o}-r^*_{i}$ is the gap.

All variables are rendered dimensionless using $d$ and $d^2/\nu$ as units for space and time, respectively. In cylindrical coordinates $(r,\theta ,z)$ , the velocity $\boldsymbol {v}=(v_r,v_{\theta },v_z)=v_r\,\hat {\boldsymbol {r}} + v_{\theta }\,\hat {\boldsymbol { \theta }} + v_z\,\hat {\boldsymbol {z}}$ and pressure $p$ of the fluid are governed by the Navier–Stokes equations, the incompressibility condition and the zero axial net massflux condition, i.e.

(2.1) \begin{align} &\partial _{t}\boldsymbol {v} + (\boldsymbol {v}\cdot \nabla )\boldsymbol {v} = -\nabla p + \nabla ^{2}\boldsymbol {v} + f\,\hat {\boldsymbol {z}}, \end{align}

(2.2) \begin{align} &\qquad\qquad\qquad \nabla \cdot \boldsymbol {v} = 0, \end{align}

(2.3) \begin{align} &Q(\boldsymbol {v}) = \int _0^{2\pi }{\int _{r_{i}}^{r_{o}}{(\boldsymbol {v}\cdot \hat {\boldsymbol {z}})\,r\,{\textrm{d}}r\,{\textrm{d}}\theta }} = 0, \\[6pt] \nonumber \end{align}

where the axial forcing term $f=f(t)$ in (2.1) is instantaneously adjusted to fulfil the constraint imposed by (2.3). The no-slip boundary conditions at the cylinder walls are

(2.4) \begin{equation} \boldsymbol {v}(r_i,\theta ,z)=(0,R_i,0),\quad \boldsymbol {v}(r_o,\theta ,z)=(0,R_o,0), \end{equation}

with $r_i=r_i^*/d=\eta /(1-\eta )$ and $r_o=r_o^*/d=1/(1-\eta )$ . The base, laminar and steady circular Couette flow, henceforth referred to as CCF, has

(2.5) \begin{equation} \boldsymbol {v}_{{b}} = V_{{b}} \hat {\boldsymbol {\theta }} = \displaystyle \left (Ar+\frac {B}{r}\right )\,\hat {\boldsymbol {\theta }}, \;\; p_{{b}}(r)=\int \frac {V_{{b}}^2}{r} \,\textrm {d}r, \;\; f_{{b}} = 0, \end{equation}

with $ A=(R_o-\eta R_i)/(1+\eta )$ and $B=\eta (R_i-\eta R_o)/ [(1-\eta )(1-\eta ^2) ]$ . In what follows, we express the velocity and pressure fields as

(2.6) \begin{equation} \boldsymbol {v} = \boldsymbol {v}_{{b}}(r) + \mathbf {u}(r,\theta ,z;t), \quad p=p_{{b}}(r) + q(r,\theta ,z;t). \end{equation}

The fields $q$ and $\mathbf {u}= u_r\,\hat {\boldsymbol {r}} + u_{\theta }\,\hat {\boldsymbol {\theta }} + u_z\,\hat {\boldsymbol {z}}$ are the deviations from the CCF solution. The nonlinear boundary value problem for $\mathbf {u}$ and $q$ is discretised using a solenoidal Petrov–Galerkin spectral scheme described by Meseguer et al. (Reference Meseguer, Avila, Mellibovsky and Marques2007). The unknown perturbation fields are approximated by means of a Chebyshev $\times$ Fourier $\times$ Fourier spectral expansion in $(r,\theta ,z)$ , but in the annular-parallelogram domain $(r,\xi ,\zeta ) \in [r_{i},r_{o} ]\times [ 0,2\pi ] \times [ 0,2\pi ]$ , where the $2\pi$ -periodicity is imposed in the transformed coordinates

(2.7) \begin{equation} \xi = n_1\theta + k_1 z,\quad \zeta = n_2\theta + k_2 z, \end{equation}

following Deguchi & Altmeyer (Reference Deguchi and Altmeyer2013). Figure 2(a) shows the computational domain corresponding to $\eta =0.883$ and the wavenumber set $(n_1,k_1,n_2,k_2)=(10,2,0,4.5)$ . The domain was specifically designed for the earliest onset of non-trivial solutions while keeping compatibility with the large-scale tilt of spiral turbulence. For this reason, the bifurcation scenario we investigate is fully compatible with the overall structure of intermittent patterns experimentally observed in the counter-rotating regime of the Taylor–Couette system and, at the same time, can be reasonably expected to precede any other mechanism potentially leading to the onset of chaos. All computations have been run in this domain. The same resolution of $[0,50]\times [-8,8]\times [-8,8]$ modes as used in W22 has been employed throughout. The system of ordinary differential equations that results from spatial discretisation has dimension $O(10^5)$ .

Figure 2. (a) Annular-parallelogram computational domain defined by the coordinates of (2.7) with wavenumbers $(n_1,k_1,n_2,k_2)=(10,2,0,4.5)$ and $\eta =0.883$ , adopted from W22. The axial line probe (red dashed vertical line) used in the production of space–time diagrams is located at mid gap $r_m=(r_i+r_o)/2 \approx 8.047$ . (b) Three-dimensional flow structure of DRW solution for $(R_i,R_o) = (450,-1200)$ . Positive (yellow, $u_{\theta } = 250$ ) and negative (blue, $u_{\theta } = -100$ ) isosurfaces of perturbation azimuthal velocity.

To explore the dynamically relevant invariant sets, we combine direct numerical simulation (DNS) with a Poincaré–Newton–Krylov (PNK) iterative scheme. The spectrally discretised Navier–Stokes equations are integrated in time using a fourth-order linearly implicit IMEX scheme. The PNK scheme, built on top of the time integrator, looks for zeroes of a map defined with a purposely devised Poincaré section by means of a Krylov-space-based Newton solver. For a detailed account of the numerical methods used, refer to Ayats et al. (Reference Ayats, Deguchi, Mellibovsky and Meseguer2020) and W22.

We characterise flow states by the normalised kinetic energy $\kappa$ of the perturbation velocity field, and by the corresponding inner and outer cylinders normalised torque, $\tau _{i}$ and $\tau _{o}$ ,

(2.8) \begin{eqnarray} \kappa = \frac {E(\boldsymbol {u})}{E(\boldsymbol {v}_b)},\qquad \tau _{{i}, {o}}= \left . 1+\frac {\partial _r(r^{-1}\langle u_{\theta }\rangle _{\xi \zeta })}{\partial _r(r^{-1}V_{b})}\right |_{r\,=\,r_{i},r_{o}}, \end{eqnarray}

where

(2.9) \begin{equation} E(\boldsymbol {v}) = \frac {1}{2\mathcal {V}} \iiint _{\mathcal {V}}{\boldsymbol {v}\cdot \boldsymbol {v}\,{\textrm{d}}\mathcal {V}} = \frac {1}{2\mathcal {V}}\int _{0}^{2\pi }\!\int _{0}^{2\pi }\!\int _{r_{i}}^{r_{o}} \boldsymbol {v}\cdot \boldsymbol {v}\,r\,\textrm{d}r\textrm {d}\xi \textrm {d}\zeta = \frac {1-\eta }{1+\eta }\int _{r_{i}}^{r_{o}}{\langle \boldsymbol {v}\cdot \boldsymbol {v}\rangle _{\xi \zeta }\,r\,\textrm {d}r} \end{equation}

is the volume-averaged kinetic energy of some velocity field $\boldsymbol {v}$ . The volume of the transformed computational domain is $\mathcal {V}=2\pi ^2(r_{o}^2-r_{i}^2)=2\pi ^2(1-\eta )/(1+\eta )$ and $\langle \ \rangle _{\xi \zeta }$ implies averaging in both parallelogram directions. With these definitions, $\kappa =0$ and $\tau _{i}=\tau _{o}=1$ for CCF.

The system possesses translational invariance in $\xi$ and $\zeta$ . We define the entire set of possible spectral coefficients as the phase space. All coefficient vectors that belong to the group orbit induced by translation of a velocity field represent the same solution. We have systematically factored out the group orbit invariance employing the method of slices (Budanur et al. Reference Budanur, Cvitanović, Davidchack and Siminos2015), using the same template as Wang et al. (Reference Wang, Mellibovsky, Ayats, Deguchi and Meseguer2023). To analyse the period-doubling cascade within the framework of discrete-time dynamical systems, we have devised a Poincaré section

(2.10) \begin{equation} \Sigma =\left \{{\mathbf{a}}\in \mathbb {X}\left |\;\tau _i({\mathbf{a}})=\tau _o({\mathbf{a}}),\;\dfrac {{\textrm{d}}\tau _i}{{\textrm{d}}t}\gt \dfrac {{\textrm{d}}\tau _o}{{\textrm{d}}t}\right . \right \}, \end{equation}

in phase space $\mathbb {X}$ , which consists of all possible sliced spectral expansion coefficients vector $\mathbf {a}$ . The condition defining $\Sigma$ is based on the equality of inner and outer cylinder torque (first condition), and the sign of the rate of change of their difference (second condition) to discard reverse crossings. It is easy to check from the governing equations that for statistically steady states, the time average of $\tau _{i}$ must coincide with that of $\tau _{o}$ . Consequently, long-lived (including permanent) time-dependent solutions are bound to regularly fulfil the condition $\tau _{i}=\tau _{o}$ , a property that comes in handy in defining a robust Poincaré section. To simplify notation, we will hereafter call $\tau =\tau _i=\tau _o$ the value of the torque, whether inner or outer, on $\Sigma$ .

There are other ways of sampling a continuous system to produce a discrete system. For example. Kreilos & Eckhardt (Reference Kreilos and Eckhardt2012) and Gao et al. (Reference Gao, Podvin, Sergent and Xin2015) investigated the relationship between consecutive local maxima of the time series of some physical quantity. However, it is not easy to discern which minima/maxima are actually related to the underlying (pseudo-)periodicity of the solutions and which are simply incidental to the dynamics of the particular time signal used. Thus, the use of the above tailor-made Poincaré section provides a more reliable approach.

3. Period-doubling cascade

Emulating the influential experiments by Andereck et al. (Reference Andereck, Liu and Swinney1986), we employ a radius ratio $\eta =0.883$ ( $r \in [r_i,r_o]=[7.547,8.547]$ ). The outer cylinder Reynolds number is fixed at $R_o = -1200$ . For sufficiently large $R_i$ , the centrifugal instability of the base flow develops into fully fledged turbulence. A reduction of $R_i$ , while still in the supercritical region, has the flow re-organise into a pattern of alternating laminar and turbulent winding helical bands, commonly known as the spiral turbulence regime (Coles Reference Coles1965; Andereck et al. Reference Andereck, Liu and Swinney1986; Dong Reference Dong2009; Meseguer et al. Reference Meseguer, Mellibovsky, Avila and Marques2009b ). At this particular value of $R_o$ , the CCF is linearly stable for $R_i \lesssim 447.35$ , but a number of subcritical nonlinear equilibrium states persist below this threshold (Meseguer et al. Reference Meseguer, Mellibovsky, Avila and Marques2009a ; Deguchi et al. Reference Deguchi, Meseguer and Mellibovsky2014; and W22), which are believed to act as precursors of spiral turbulence. Drawing from the helix angle of the turbulent stripes observed experimentally, W22 employed the domain and coordinate system shown in figure 2(a) to seek minimal flow unit solutions that are compatible with the tilt of the banded pattern. Subharmonic instabilities of some such solutions were proposed as the mechanism that might be responsible for the spatial intermittency characterising spiral turbulence. Within this minimal parallelogram-annular periodic domain, which is capable of sustaining turbulence at sufficiently high $R_i$ , W22 identified a stable finite amplitude travelling wave (shown in figure 2 b) that coexists with the stable CCF. This drifting rotating wave, DRW for short, embodies all the essential elements of the self-sustaining mechanism (Wang, Gibson & Waleffe Reference Wang, Gibson and Waleffe2007; Hall & Sherwin Reference Hall and Sherwin2010) and plays a central role in the transition process. In W22, two periodic solutions, P $_1$ , arising from a Hopf bifurcation of DRW, and P $_2$ , originating at a period-doubling bifurcation of P $_1$ , were also identified as the first steps in the route to chaos. This sequence of bifurcations suggested a period-doubling cascade as the most plausible scenario for the onset of chaotic dynamics, but the issue was not pursued further.

3.1. Onset of the period-doubling cascade

We now look into the bifurcation sequence the system undergoes in the range $R_i\in [395.43,395.79]$ eventually leading to sustained chaotic dynamics. Since we will only be varying $R_i$ , we shall drop the subscript and simply write $R$ to ease notation.

Figure 3(a) shows the DRW (square), P $_1$ (dashed blue line) and P $_2$ (solid green line) at $R=395.67$ in a three-dimensional projection of the full phase space $\mathbb {X}$ on the subspace spanned by the triplet of key quantities $(\tau _o,\tau _i,\kappa )$ . The Poincaré section $\Sigma$ appears in this representation as the transparent grey plane, which contains DRW and is pierced, in the direction defined by (2.10), at a single point by P $_1$ and at two different points by P $_2$ . In the torque time series of figure 3(b), the Poincaré crossings are conveniently identified as the intersections between the $\tau _i$ (thick lines) and $\tau _o$ (thin) signals for which the former is overtaking the latter. The two-dimensional projection of figure 3(a) onto the plane $(\tau _o,\tau _i)$ , as depicted in figure 3(c), further clarifies how the sequence of intersection points of figure 3(b) collapses in phase space onto a single point for P $_1$ (empty blue disk) and as two distinct points for P $_2$ (filled green circles), all contained in $\Sigma$ (dashed grey straight line).

Figure 3. Poincaré section $\Sigma$ and the periodic orbits P $_1$ (dashed blue) and P $_2$ (solid green). The black squares are DRW solutions. All solutions are computed at $R=395.67$ . (a) Projection of the phase space on the $(\tau _o,\tau _i,\kappa )$ coordinates. (b) Inner ( $\tau _i$ , thick lines) and outer ( $\tau _o$ , thin) torque time series of P $_1$ and P $_2$ . (c) Two-dimensional phase map projection on the $(\tau _o,\tau _i)$ plane. The Poincaré section is shown in transparent grey in panel (a) and as a dashed grey line in panel (c). The circles on the P $_1$ (empty blue) and P $_2$ (filled green) curves correspond to their representation on $\Sigma$ .

The drift dynamics of non-axisymmetric solutions in Taylor–Couette as we have here is inevitably masked when monitoring aggregate quantities such as torque or kinetic energy due to their spatial averaging properties. Point measurements are instead subject to drift-induced time dependence. Figure 4(ai) presents a space–time diagram of the axial vorticity distribution measured along a line probe that is fixed in the lab reference frame (red dashed line in figure 2 a). DRW, a relative equilibrium, is characterised by a solid-body motion composing axial translation and azimuthal rotation with constant phase velocities. As a result, the space–time diagram exhibits the repetition of a periodic pattern. The same line probe produces a quasi-periodic space–time diagram for P $_1$ , as shown in figure 4(bi).Here, P $_1$ is a relative periodic orbit and, therefore, requires suitable shifts in both the $z$ and $\theta$ directions to align the flow fields after every one period to reveal the space–time invariance. The effects of the drift can be suppressed by attaching the line probe to a moving reference frame defined with the method of slices (Budanur et al. Reference Budanur, Cvitanović, Davidchack and Siminos2015; Wang et al. Reference Wang, Mellibovsky, Ayats, Deguchi and Meseguer2023). In this reference frame, the space–time diagrams of DRW and P $_1$ (figures 4 aii and 4 bii) are greatly simplified, the former appearing as time-independent and the latter as purely time-periodic.

Figure 4. Space–time diagrams for (a) DRW, (b) P $_1$ and (c) P $_2$ , all computed at $R=395.67$ . The roman number labels denote measurements of radial vorticity $\omega _r(z;t)$ along axial probe lines at $(r,\theta )=(r_m,\theta _0)$ fixed to (i) the lab (stationary) reference frame, (ii) a reference frame co-moving with the solution and (iii) the same co-moving frame but with the temporal mean $\langle \omega _r\rangle _{t}$ subtracted. The azimuthal location, $\theta _0$ , is chosen consistently across reference frames and solutions to enable comparison. Colour shading according to $\omega _r\in [-1400,1400]$ or $\omega _r- \langle \omega _r\rangle _{t}\in [-300,300]$ , as need be. Dashed vertical lines indicate the natural period of the corresponding solution.

Further subtraction of the time average helps expose the true nature of the time dependence of the solutions. For P $_1$ , the time-horizon considered in figure 4(biii) allows for just over four repetitions, as indicated by the dashed vertical lines. Solution P $_2$ , shown in figure 4(ciii), has instead a natural period about twice that of P $_1$ , the period-doubling consisting in a modulation that shortens one of the half-periods while lengthening the other.

Let us now shift our focus to the dependence of the solutions on the parameter $R$ . Figure 5(a) shows the bifurcation sequence that DRW undergoes when varying the inner Reynolds number in the range $R\in [391.20,395.90]$ . DRW emerges from a saddle-node bifurcation (SN $_1$ ) at $R\simeq 391.50$ and, leaving aside subharmonic instabilities, the nodal (upper) branch (the one represented in figure 2 b and the square in figure 3 a, 3 c) remains stable to perturbations fitting the domain until the advent of a supercritical Hopf bifurcation (H) at $R\simeq 392.85$ . As reported by W22, this is the only known case of a stable non-trivial solution in the subcritical parameter regime of Taylor–Couette flow. Note that, here, we are only interested in the stability to perturbations that fit within the periodic box. Stable/unstable solutions are denoted by solid/dashed lines. The PNK method must be used to compute unstable solution branches. A stable relative periodic orbit (P $_1$ , blue line) bifurcates from DRW at H. All time-dependent solutions will be represented, from now on, through the collection of intersection points on the Poincaré section. The P $_1$ branch becomes unstable in a period-doubling bifurcation (PD $_1$ ), whence a branch of stable period-doubled relative periodic orbits (P $_2$ , pair of green lines) emerges. The P $_2$ loses stability, in turn, in a second period-doubling bifurcation (PD $_2$ ) issuing a branch of period-4 solutions. The $R$ value for which figures 3 and 4 were computed corresponds to just short of PD $_2$ , hence the instability of DRW and P $_1$ , in contrast with P $_2$ , which is stable.

Figure 5. Bifurcation scenario as recorded on the Poincaré section $\Sigma$ . (a) Initial steps of the bifurcation scenario. Shown are DRW (black), P $_1$ (blue) and P $_2$ (green), reported in W22. P $_4$ (green) emerges at the second period-doubling bifurcation point PD $_2$ . Both stable (solid line) and unstable (dashed) solution branches are shown. (b) Detailed view (close-up of the region bounded by a solid grey box in panel a) of stable solution branches across the period-doubling cascade and beyond. The accumulation point for the period-doubling cascade ( $R_\infty$ ) is to be computed in § 4.

The remainder of the period-doubling cascade and the onset of chaotic dynamics at larger values of $R$ is depicted in figure 5(b). Only stable solutions are shown, so DNS is all that has been required to produce the solution branches in the diagram. Checking agreement with $\delta _F$ from the approximate location on the bifurcation diagram of the period-doubling points may seem a straightforward task. However, confirmation of Feigenbaum universality by brute force poses challenges, even for simple systems such as the logistic map. To begin with, the parameter spacing between consecutive period-doubling bifurcations shrinks very fast as one progresses along the cascade, demanding the computation of bifurcation points with an ever increasing number of significant digits. Moreover, the orbits in the immediate vicinity of a period-doubling point, which are required for the accurate estimation of the bifurcations, are close to neutrally stable and, therefore, take massive computational time to reach convergence with DNS. Unfortunately, employing the PNK method becomes impractical for solutions of increasingly long periods. The inability of the method to discriminate between stable and unstable orbits, combined with the fact that stable orbits in a period-doubling cascade coexist with all unstable orbits of lower period, which are favoured unless very close initial guesses are produced in advance, curbs any advantage one might have expected to achieve from using PNK. Both complications combined render the careful appraisal of agreement with Feigenbaum’s first constant extremely hard. A rigorous and systematic analysis, as detailed in the coming section, becomes thus necessary.

3.2. Determination of period-doubling points

To analyse the period-doubling cascade, it is most convenient to focus on the sequence of $\bf a$ on $\Sigma$ . In the framework of dynamical systems theory, this corresponds to investigating the properties of the Poincaré map, also known as first return map. As we shall explain later, it is not necessary to monitor the full set of coefficients, and the sequence of torque values provides (nearly) all necessary information in the long run, once the initial transients are over. Let $\tau (\ell )$ , $\ell =1,2,3,\dots$ be the corresponding sequence of torque values on $\Sigma$ , with $\ell$ an index recording the chronological order. As an example, figure 6(a) shows $\tau (\ell )$ for the P $_4$ solution (red circles) at $R=395.711$ (between PD $_2$ and PD $_3$ ) and for the P $_8$ solution (blue) at $R=395.719$ (between PD $_3$ and PD $_4$ ). Both solutions are stable and well converged with DNS, such that the respective series sequentially repeat, always in the same order, the same four (P $_4$ ) and eight (P $_8$ ) distinct values, as indicated by the numbered horizontal dashed lines. The effect of the period-doubling bifurcation PD $_3$ can thus be portrayed as doubling point $j$ of the P $_4$ discrete-time orbit into points $j$ and $j+4$ of the P $_8$ cycle, the latter two points having sprung from the former and drifted away in opposite directions. The unstable P $_1$ (cyan) and P $_2$ (green) orbits in the immediate vicinity of PD $_3$ are shown for reference, to help generalise the rule that relates points $j$ and $j+N/2$ of the period- $N$ orbit to point $j$ of the period- $N/2$ orbit from which they bifurcate. A phase map analogous to that of figure 3(c) is shown in figure 6(b), but with the region where $\Sigma$ is pierced by the orbits suitably magnified to further clarify where the P $_4$ and P $_8$ cycles stand in relation to P $_1$ and P $_2$ .

Figure 6. All orbits up to period 8 at (P $_1$ and P $_2$ , unstable) and around (P $_4$ and P $_8$ , at $R=395.711$ and $395.719$ , respectively) PD $_3$ . (a) Torque ( $\tau =\tau _i=\tau _o$ ) of P $_1$ (cyan), P $_2$ (green), P $_4$ (blue) and P $_8$ (red) on $\Sigma$ as a function of the discrete time $\ell$ (crossing index). The dashed lines indicate the distinct values of $\tau$ . (b) Two-dimensional phase map projection on the $(\tau _o,\tau _i)$ plane. Shown are the phase map trajectories of all four orbits (dashed line for unstable, solid for stable) along with their representation on the Poincaré section (circles, open for unstable, filled for stable). The numbers indicate the order of the crossings.

The $n$ th period-doubling bifurcation PD $_n$ can be readily analysed by splitting the $\tau (\ell )$ sequences of interest into $J=2^n$ separate subsequences

(3.1) \begin{equation} \tau _j^J(k)=\tau (kJ+j), \end{equation}

each starting at one of $J$ consecutive Poincaré crossings $j\in \{1,\ldots J\}$ , and sampling every $J$ th crossing thereafter. The limits

(3.2) \begin{equation} (\tau _j^J)_\infty = \lim _{k\to \infty }\tau _j^J(k) \qquad j\in \{1,\ldots J\} \end{equation}

exist for all subsequences of all orbits of period $N=J=2^n$ or lower ( $J/2,J/4,\dots$ ) along the cascade, but fail for orbits of longer period ( $2J,4J,\dots$ ). Moreover, only period- $J$ orbits will produce $J$ distinct $(\tau _j^J)_\infty$ values, while shorter period orbits will have the limits coincide in pairs according to $(\tau _j^J)_\infty =(\tau _{j+J/2}^J)_\infty$ . This property is instrumental in telling apart orbits of period $N=J$ from orbits of period $J/2$ or smaller, as the amplitudes defined by

(3.3) \begin{equation} A_j^J = \left \lvert (\tau _j^J)_\infty -(\tau _{j+J/2}^J)_\infty \right \rvert \qquad j\in \{1\ldots J/2\}, \end{equation}

will all vanish for any orbit other than P $_J$ . For instance, suppose we slightly perturb the stable P $_8$ at $R=395.719$ . If we sample the $\tau$ sequence with $J=8$ , all subsequences $j=1,2,3,\dots ,8$ converge and all eight limits are different (recall figure 6). Applying the same $J=8$ sampling to a stable P $_4$ (or P $_2$ or P $_1$ ) will instead produce indistinguishable limits for $j$ and $j+4$ and, consequently, vanishing amplitudes $A_j^8=0$ .

To nail down PD $_3$ , a collection of DNS runs traversing the bifurcation point, i.e. between $R=395.711$ and $395.719$ , are required. The amplitude of P $_8$ drops fast as the parameter $R$ is decreased towards PD $_3$ , such that distinguishing it from P $_4$ becomes increasingly difficult. In addition, the dynamics become despairingly slow as the bifurcation point is approached from either side, rendering convergence with DNS downright impracticable. Figures 7(a) and 7(b) depict the torque subsequences at $R=395.7116$ and $R=395.7122$ , respectively, below and above PD $_3$ but moderately close to the bifurcation point. Every pair of sequences ( $\tau _j^8$ and $\tau _{j+4}^8$ , depicted in different shades of the same colour) appears to be converging to the same value at $R=395.7116$ , but not at $R=395.7122$ . At these values of  $R$ , the achievable degree of convergence is still reasonable but sufficiently slow to illustrate how the estimation of the amplitudes may be systematised to parameter values much closer to the bifurcation point. Well past initial transients, and as the sequences approach convergence, the dynamics is progressively determined by the governing equations linearised around the stable solution. Therefore, we can expect that a power law fit of the form

(3.4) \begin{equation} \tau _j^J(k) =(\tau _j^J)_{\infty }+ \left ((\tau _j^J)_0-(\tau _j^J)_{\infty }\right )(\lambda _j^J)^k, \end{equation}

with fitting parameters $(\tau _j^J)_{\infty }, (\tau _j^J)_0$ and $\lambda _j^J$ , conforms reasonably well to the final transients. The fit provides an estimate for both the asymptotic torque value $(\tau _j^J)_{\infty }$ and the dominant multiplier $\lambda _j^J$ . The dashed lines in the insets of the top panels of figure 7 are the fits to the sequences $\tau _j^J(k)$ , $k\geqslant 0$ . To expose the convergence rate and the accuracy of the fits, the same data have been plotted in the mid panels, now subtracting from each pair of branches $j$ and $j+4$ the mean value $((\tau _j^8)_\infty +(\tau _{j+4}^8)_\infty )/2$ at infinity. The fits provide an excellent approximation to the data points because the dynamics has already reached the linear regime. While all fitting curves decay to zero for $R=395.7116$ in this representation, four distinct amplitudes arise for $R=395.7122$ . This unequivocally identifies the former dynamics as converging on a P $_4$ and the latter on a P $_8$ . The bottom panels plot the same data yet again, but now in logarithmic scale. In this representation, all fits to the individual torque subsequences are seen to collapse onto a single straight line. The eight slopes $\lambda _j^8$ , indistinguishable from one another to the naked eye, provide an estimate to the dominant multiplier of the converging periodic solution, interpreted as a P $_8$ regardless of its actual period.

Figure 7. Analysis of torque sequences $\tau (\ell )$ , sampled with $N=8$ at (a) $R=395.7116$ and (b) $R=395.7122$ , below and above PD $_3$ ( $R=395.711841564$ , as we shall see later), respectively. The top panels show the eight sequences $\tau _j^8(k)$ (circles, magnified in the insets), alongside respective power law fits (dashed lines). To differentiate the eight subsequences, the points corresponding to each are coloured and represented with disks that are mainly blank except for a 1/8th sector, whose orientation ( $360^\circ \times j/8$ ) uniquely identifies the corresponding sequence $\tau _j^8$ (see the legend). The mid panels show the same data but with the mean value $\overline{\tau}_j^8 = [(\tau_j^8)_\infty + (\tau_{j+4}^8)_\infty] /2$ subtracted from every pair of branches to illustrate convergence. The bottom panels show again the same data but in logarithmic scale: $\tilde{\tau}_j^8 = [\tau_j^8 - (\tau_j^8)_\infty]/[(\tau_j^8)_0 - (\tau_j^8)_\infty]$ .

The four amplitudes $A_j^8$ and eight multipliers $\lambda _j^8$ have been computed, following the same procedure with the sampling $J=8$ , at several values of $R$ approaching PD $_3$ from either side, and plotted in figure 8. Figure 8(a) indicates that there is a well-defined critical value of $R$ , beyond which, all four $A_j^8$ start growing following a square root law, as expected for a period-doubling bifurcation. To accurately pinpoint the critical value of the parameter for PD $_n$ , a fit of the form

(3.5) \begin{equation} A_j^J(R) = a_j^J\sqrt {R-R_j^J} \end{equation}

with $J=2^n$ applied to the first few points after the period-doubling bifurcation. The fits provide, besides the $a_j^J$ values for the $J/2$ scaling factors, $J/2$ independent estimates $R_j^J$ for the critical value of the parameter. A unique estimate for the critical value is then obtained by taking the arithmetic mean of the $J/2$ individual estimates as $R_n=\langle R_j^J\rangle$ . Performing the $J/2=4$ fits to the data points of figure 8(a) (dashed lines), yields $R_3=\langle R_j^8\rangle =395.71184156\pm 5\times 10^{-8}$ (violet square) for PD $_3$ .

Figure 8. Analysis of PD $_3$ . (a) Torque amplitudes $A_j^8$ and (b) multipliers $\lambda _j^8$ (coloured, sampling $\tau (\ell )$ with $N=8$ ) and $\lambda _j^4$ (black, sampling with $N=4$ ), as a function of inner cylinder Reynolds number $R-\langle R_j^8\rangle$ . The square-root fits are indicated with dashed lines. The lower bounds for the uncertainty of each point, shown as error bars, are sufficiently small to be imperceptible. Following the graphical representation of the $\tau _j^8$ sequences of figure 7, each of the amplitudes is represented with a double-sectored disk that results from the superposition of the oppositely oriented sectored disks of the two torques that define the amplitude, hence the hourglass appearance of the symbols.

The eight $\lambda _j^8$ multipliers, shown in figure 8(b), behave as expected across a period-doubling bifurcation. The eight individual fits (3.4) to the eight individual $\tau _j^8$ subsequences produce nearly identical values of $\lambda _j^8$ at any given $R$ , hence the multicoloured filled circles, resulting from the superposition of the eight oriented sectored disks. The unique estimate of the multiplier at any given $R$ , obtained as the arithmetic mean $\lambda _8=\langle \lambda _j^8\rangle$ , increases towards unity as the period-doubling point is approached from either side. Since the solution above PD $_3$ is indeed a P $_8$ , $\lambda _8$ is directly its multiplier. This is not so for the P $_4$ solution below PD $_3$ . Plots analogous to those in figure 7(a) but sampling with $J=4$ show that convergence onto the four asymptotic branches follows an oscillatory pattern. The multiplier of the P $_4$ orbit in the vicinity of PD $_3$ must be either computed from fits to the torque subsequences sampled with $J=4$ or indirectly from $\lambda _8$ according to $\lambda _4=\langle \lambda _j^4\rangle =-\sqrt {\lambda _8}$ . Either way, the multiplier of the stable P $_4$ in the vicinity of PD $_3$ is real, negative and on the way of leaving the unit circle of the complex plane through $\lambda =-1$ as $R_3$ is approached from below (black disks in figure 8 b), consistent with the period-doubling nature of the bifurcation.

4. Verification of Feigenbaum universality

Signs that the period-doubling cascade may exhibit some sort of universality are already apparent from figure 9, where all orbits computed up to period $N=128$ (past PD $_7$ ) are shown as a function of $R$ . Successive magnifications, presented in figures 9(b) and 9(c), exhibit branch structures very similar to that of the full cascade of figure 9(a). Yet, at this stage, there is no guarantee that the universality observed here is of the same kind as discovered by Feigenbaum.

Figure 9. Successive magnifications of the period-doubling cascade. Same data as in figure 5(b), truncated at the accumulation point $R_{\infty }$ , the central branch indicated with black dots. (a) Overview and (b,c)two successive levels of magnification. The first few period-doubling bifurcations ( $R_n$ , dashed) and the accumulation point ( $R_\infty$ , solid) are indicated with grey lines and labels.

To ascertain Feigenbaum universality, we have computed the sequence of critical $R_n$ values along the period-doubling cascade up to order $n=7$ to very high accuracy, using the method described in § 3.2. The results are summarised in table 2. For every $n\geqslant 3$ , the table also lists the ratios $\delta _n$ , defined as

(4.1) \begin{equation} \delta _n = \dfrac {R_{n-1}-R_{n-2}}{R_{n}-R_{n-1}}, \end{equation}

which provide the best approximation to Feigenbaum’s first constant $\delta _{{F}}$ that may be obtained from data up to the $n$ th period doubling. The trend of $\delta _n$ exhibits evident signs of conforming to Feigenbaum universality, that is, $\delta _n \rightarrow \delta _{{F}}$ as $n\rightarrow \infty$ . It is instructive to compare these results with those in table 1. Limiting the analysis to $n=4$ , as customarily done in the fluid dynamics literature, does not guarantee convergence to even the second significant digit. Securing three digits precision is an arduous task that typically requires the accurate determination of at least $n=7$ period-doubling bifurcation points.

Table 2. Confirmation of Feigenbaum universality along the period-doubling cascade. The values of the critical inner cylinder Reynolds number $R_n$ at the $n$ th period doubling bifurcation point PD $_n$ are used to compute the $n$ th approximation, $\delta _n$ , to Feigenbaum’s first constant, according to (4.1). Approximations, $\alpha _n$ , to the second constant, are computed from central branch torque values following (4.8). Lower bounds to the uncertainties in the $R_n$ , $\delta _n$ , $\alpha _n$ and $R_\infty$ parameters have been rigorously estimated from the covariance matrices of the various fits involved in the process, combined with standard error propagation theory (note that we have not taken into account the error of the fit (3.4), as conducting a systematic study is difficult). The last row of the table corresponds to the accumulation point $R_{\infty }$ estimated by (4.7), and the actual values of Feigenbaum constants, $\delta _{{F}}$ and $\alpha _{{F}}$ .

The $R_n$ sequence is expected to converge onto the accumulation point, $R=R_\infty$ . This point may be computed by assuming that all period-doubling bifurcations of order higher than available are perfectly Feigenbaum-universal and, therefore, satisfy (4.1) with the left-hand side replaced by $\delta _{{F}}$ . However, this naive approach leads to systematic error. To obtain a better estimate of $R_{\infty }$ and, at the same time, bound the error, we first note that

(4.2) \begin{eqnarray} R_{\infty } = R_n+\sum _{m=n+1}^{\infty }{\Delta R_m}, \end{eqnarray}

where $\Delta R_m\equiv R_m-R_{m-1}$ denotes the interval between the $(m-1)$ and $m$ th period-doubling bifurcations. Since (4.1) can be written as $\delta _n=\Delta R_{n-1}/\Delta R_{n}$ , any successive $m$ th interval may be expressed recursively in terms of previous intervals as

(4.3) \begin{equation} \Delta R_{m}=\delta _m^{-1}\Delta R_{m-1}=\delta _m^{-1}\delta _{m-1}^{-1}\Delta R_{m-2}=\dots =\delta _m^{-1}\delta _{m-1}^{-1}\cdots \delta _{n+1}^{-1}\Delta R_{n}=\Delta R_{n}\prod _{\ell =n+1}^{m}\delta _\ell ^{-1}. \end{equation}

The accumulation point can then be formally expressed in terms of the n $th$ period doubling and the immediately preceding interval as

(4.4) \begin{eqnarray} R_{\infty } = R_n+\Delta R_n\left [\sum _{m=n+1}^{\infty }{\prod _{\ell =n+1}^{m}\delta _\ell ^{-1}}\right ]. \end{eqnarray}

If $\delta _\ell =\tilde {\delta }_n\gt 1$ were constant for all $ \ell \gt n$ , then the terms in the square brackets would form a geometric series of sum $1/(\tilde {\delta }_n-1)$ . We conveniently define the average ratio beyond PD $_n$ as

(4.5) \begin{equation} \tilde {\delta }_n = 1 + \left [ \sum _{m=n+1}^{\infty }\displaystyle {\prod _{\ell =n+1}^{m}{\delta _\ell }^{-1}} \right ]^{-1}, \end{equation}

so that (4.4) simplifies to

(4.6) \begin{equation} R_\infty = R_n + \dfrac {\Delta R_n}{\tilde {\delta }_n-1} = \dfrac {\tilde {\delta }_n R_n - R_{n-1}}{\tilde {\delta }_n-1}. \end{equation}

The constant $\tilde {\delta }_n$ is unknown and may only be estimated. For $n$ sufficiently large, the sequence $\delta _\ell$ for all $\ell \gt n$ is expected to approach $\delta _{{F}}$ monotonically from above. Under this assumption, inspection of (4.5) provides the inequalities $\delta _{{F}}\lt \delta _{n+1}\lt \tilde {\delta }_n \lt \delta _n$ . Accordingly, upper and lower bounds for $R_\infty$ can be obtained as

(4.7) \begin{equation} R_\infty \in \left [R_n+\dfrac {R_n-R_{n-1}}{\delta _n-1},R_n+\dfrac {R_n-R_{n-1}}{\delta _{{F}}-1}\right ], \end{equation}

by substituting $\tilde {\delta }_n=\delta _n$ and $\tilde {\delta }_n=\delta _{{F}}$ in (4.6), respectively. Using PD $_7$ , the upper and lower bounds for $R_\infty$ are already within the uncertainty range of one another (see table 2) and, combined, provide an accurate estimate for the accumulation point $R_\infty =395.721266624\pm 1.1\times 10^{-8}$ .

Figure 10. Self-similarity of the period-doubling cascade. Same data as figure 9 in log scale, the central branch represented with black dots. The coloured bottom panels are successive magnifications of the boxes in the top panel. The bifurcations are indicated with dashed vertical lines and labelled PD $_n$ according to their order.

Following the determination of the accumulation point, we can now plot the period-doubling cascade as a function of $R_\infty -R$ in logarithmic scale (figure 10). The asymptotic approach to universality and the accuracy with which the accumulation point has been obtained are evident from the similarities observed after two consecutive magnifications of the period-doubling cascade (red and blue bottom panels).

The occurrence of Feigenbaum universality implies that the dynamics of the Navier–Stokes system is governed by a (nearly) one-dimensional discrete map. Figure 11(a) depicts a three-dimensional phase map projection of the orbit at the accumulation point $R_\infty$ . If the manifold has dimension one, the coefficients $\mathbf {a}$ on the Poincaré section $\Sigma$ introduced in (2.10) can be uniquely parametrised by the torque, such that the first return map $f$ defined by $\tau (\ell +1)=f(\tau (\ell ))$ contains all the long-term properties of the Poincaré map. The aperiodicity of the orbit becomes apparent when plotting $\tau (\ell +1)$ as a function of $\tau (\ell )$ (green dots in figure 11 b). The structure of the map on $\Sigma$ is, however, not clear, as the $\Sigma$ -crossings of the orbit are sparsely distributed (see figure 11 a). Perturbing the orbit, dropping the first initial transients and plotting the remainder of the transients conveniently assists in populating the map more densely (black points). All the points thus obtained fall within a narrow, nearly one-dimensional band, seemingly connected by a smooth curve. This indicates that the discrete-time dynamics on $\Sigma$ in the vicinity of the accumulation point can be well approximated by a one-dimensional map $f$ . Note that, in the range of $\tau$ for which $f$ appears to be multivalued (see the black points), the lower branch is never visited by the attractor and can, therefore, be ignored as regards its dynamics. Accordingly, the function $f$ , with its leftmost branch duly trimmed, becomes unimodal and its smooth extremum (the minimum) admits a quadratic approximation. Note that unimodal maps with non-quadratic extrema provide universality constants that deviate from Feigenbaum’s well-known values (Briggs Reference Briggs1991).

Figure 11. Phase portrait at $R_{\infty }$ . (a) Three-dimensional phase map projection of the chaotic attractor. (b) First return iteration map on $\Sigma$ populated by the chaotic attractor (green dots) and transients leading back to it following a perturbation (black).

For period-doubling cascades in one-dimensional unimodal maps, approximations of the second Feigenbaum constant, $\alpha _{{F}}$ , are often obtained from the central branch of the bifurcation diagram, which is generated by the periodic point that is nearest to the extremum of the map. Since producing accurate bifurcation diagrams of cascades undergone by fluid flow systems is an arduous task, Feigenbaum (Reference Feigenbaum1979) used instead the frequency spectrum, recorded at a fixed value of the parameter, of a natural convection experiment. In the subsequent studies summarised in table 1, no evaluation of agreement with $\alpha _{{F}}$ was conducted.

However, we have a highly accurate bifurcation diagram available, rendering direct assessment of agreement with the second constant possible. In the torque sequences presented so far, we have indexed the $N$ distinct points visited by every periodic solution P $_N$ as $j=1,2,\dots N$ such that $j=1$ is the central branch point (the black dots in figure 10). Thus, using the torque values of the central branch at PD $_m$ , $\hat {\tau }_{m}=\tau _1^{2^m}(R_m)$ , we can compute the ratios

(4.8) \begin{eqnarray} \alpha _n = \dfrac {\hat {\tau }_{n-1}-\hat {\tau }_{n-2}}{\hat {\tau }_{n}-\hat {\tau }_{n-1}} \end{eqnarray}

to be compared against the $\alpha _{{F}}$ appearing in the Feigenbaum–Cvitanović functional equation. The values of $\alpha _n$ thus obtained from the bifurcation diagram, listed in table 2 alongside $\delta _n$ , evidence gradual convergence toward Feigenbaum universality.

5. Prediction of the period-doubling route to chaos

Following the seminal work of Feigenbaum (Reference Feigenbaum1979), who showcased the first comparison with experimental results, Feigenbaum (Reference Feigenbaum1980, Reference Feigenbaum1982, Reference Feigenbaum1983 Reference Feigenbaum) further extended renormalisation theory to include the universal scaling of all periodic points, not limited to the central branch, at different values of the parameter. The trajectory scaling function theory, which unravels the scaling between orbits featuring the same stability properties (i.e. having the same multiplier), was later applied to interpret experimental data on natural convection by Belmonte et al. (Reference Belmonte, Vinson, Glazier, Gunaratne and Kenny1988).

The occurrence of universality entails that, for any given unimodal map, one should be able to predict the entire structure of the bifurcation diagram near the accumulation point using data from only the initial few period-doubling bifurcations. Unfortunately, the trajectory scaling function method does not provide an easy way to predict the bifurcation cascade. Moreover, orbits become aperiodic beyond the accumulation point, hindering direct comparison with orbits at different parameter values. In this section, we devise a method for predicting the unfolding of the bifurcation cascade up to and beyond the accumulation point.

Consider a general one-dimensional discrete-time dynamical system described by

(5.1) \begin{eqnarray} x_{\ell +1}=f_{a}(x_\ell ). \end{eqnarray}

We assume that the function $f_{a}(x)$ is unimodal and that a fixed point of the system (5.1) undergoes a period-doubling cascade as the parameter $a$ is increased. In drawing parallels with our Navier–Stokes problem, we will later identify $a$ and $f_a$ with the Reynolds number $R$ and the function approximating the return map of the torque sequence, respectively. Our claim, motivated by Feigenbaum (Reference Feigenbaum1982), is that the approximation

(5.2) \begin{eqnarray} \frac {f^{2^{n}}_{a}(x)-X}{\gamma } \approx U^{2^n}_{M} \left ( \frac {x-X}{\gamma } \right ), \qquad M= \frac {a-a_{\infty }}{\mu } \end{eqnarray}

holds as long as $2^n$ is sufficiently large,

(5.3) \begin{eqnarray} \left |\frac {\delta ^{n} (a-a_{\infty })}{\mu } \right |\ll 1, \qquad \left | \frac {\alpha ^n(x-X)}{\gamma } \right | \lt 1, \end{eqnarray}

and provided the constants $\gamma$ , $X$ , $\mu$ and $a_{\infty }$ are suitably chosen. The derivation of (5.2) and (5.3) is given in Appendix. To simplify the notation, we use functional powers to denote repeated application of a function (e.g. $f^2_a(x)=f_a(f_a(x))$ ). Also, we have dropped the F subscript when referring to Feigenbaum’s first ( $\delta$ ) and second ( $\alpha$ ) constants.

Here, the function $U_M$ appearing in the approximation (5.2) is defined in terms of the Feigenbaum function $G$ and the Feigenfunction $\varPhi$ (with the normalisation $\varPhi (0)=1$ ) as

(5.4) \begin{eqnarray} U_M(\xi ) = G(\xi )+M\varPhi (\xi ). \end{eqnarray}

The $\delta$ and $\varPhi$ satisfy the eigenvalue problem $\mathcal {L}_G[\varPhi ]=\delta \varPhi$ . The precise form of the linear operator $\mathcal {L}_G$ is given in Appendix. The nature of $G$ and $\varPhi$ is well understood (Collet & Eckmann Reference Collet and Eckmann1980; Briggs Reference Briggs1991; Thurlby Reference Thurlby2021). For our purposes, it will suffice to know that they are both smooth even functions, and that numerical recipes for their computation abound in the literature.

Figure 12. Bifurcation diagram of the one-dimensional map (5.5). (a) Full period-doubling cascade. (b) Magnification around the accumulation point $M_\infty$ . The $2^n=4$ distinct sets of points arising for $M\gt M_*$ are labelled as $\xi _j^n$ , $j=1,2,3,4$ . Black dots are used for the central branch, grey for the rest. (c) Close-up of the central branch $\overline {\xi }=\xi _1^n$ . The red bullets at PD $_3$ and PD $_4$ are the points $(M,\xi )=(M_{n+1},\hat {\xi }_{n+1})$ and $(M_{n+2},\hat {\xi }_{n+2})$ selected for the matching with DNS data. (d) Full cascade transformed by (5.8) with $(\gamma ,X,\mu ,a_{\infty })$ obtained by (5.7), compared with the period-doubling cascade in Taylor–Couette flow (green dots).

If the four (map-dependent) parameters, $\gamma$ , $X$ , $\mu$ and $a_{\infty }$ , can be computed from the first few period doublings, then (5.2) provides a useful tool for predicting the behaviour of the dynamical system (5.1), subject to the constraints outlined by (5.3). In practice, fixing $n=2$ and using data up to PD $_{n+2}$ (i.e. PD $_4$ ) yields already reasonably accurate results. The procedure for prediction is straightforward, as we shall see shortly. All it takes is extracting template branches from the bifurcation diagram of the one-dimesnional dynamical system

(5.5) \begin{eqnarray} \xi _{\ell +1}=U_M(\xi _\ell ), \end{eqnarray}

and then stretch and translate them to align with corresponding branches of the bifurcation diagram of (5.1). The constants $X$ and $\gamma$ translate and scale the state $x$ , while $a_\infty$ and $\mu$ do the same for the parameter $a$ . The constant $a_{\infty }$ corresponds to the predicted accumulation point of the dynamical system (5.1). Thus, the first condition in (5.3) requires that the parameter $a$ be sufficiently close to the accumulation point. Likewise, $X$ corresponds to the predicted position of the extremum of the map $f_a(x)$ at $a=a_{\infty }$ , and hence the second condition in (5.3) limits the validity of the approximation to a close neighbourhood of the extremum of the map $f_a(x)$ , i.e. the central branch of the cascade. In brief, the conditions (5.3) imply that universality applies to the bifurcation diagram only locally. We will see later how the prediction can, nevertheless, be extended to the rest of the branches.

The bifurcation diagram of the dynamical system (5.5) is depicted in figure 12(a). It is easy to show that, for $M=-1$ , the system has a stable fixed point at the origin ( $\xi =0$ ) and that $M=M_\infty =0$ is the accumulation point of the period-doubling cascade. Assume that fixing the parameter $M$ to a small (in modulus) negative value produces a stable orbit of, say, period $N=2^n$ with $n\in \mathbb {N}$ . An example of $M=M_*$ producing a period-4 orbit ( $n=2$ ) is shown in figure 12(b), which contains a magnification of figure 12(a) around the accumulation point. Each of the $2^n$ branches at $M=M_*$ initiates its own period-doubling cascade as $M$ is increased from $M_*$ towards $M_\infty$ . After the accumulation point, the $2^n$ distinct clouds of points induced by the $2^n$ separate branches start merging in pairs following a reverse cascade until eventually forming one sole object. We have labelled the four cascades as $\xi _j^4$ ( $j=1,2,3,4$ ) in figure 12(b) and shown them up to the point where they first merge in two pairs ( $\xi _1^4$ with $\xi _3^4$ and $\xi _2^4$ with $\xi _4^4$ ). We shall see that each of these cascade branches can act as a template in our prediction, but first, let us focus on the cascade induced by the central branch, henceforth denoted as $\overline {\xi }=\xi _1^4$ (see figure 12 c).

To illustrate how the approximation (5.2) operates, we begin by fixing $M=M_*$ in figure 12(b), where we have an orbit of period $N=2^n=4$ , so $n=2$ . Each template $\xi _j^4$ corresponds to a single point at this value of the parameter, with the $j$ labels following the same sampling procedure of (3.1) with $J=2^2$ . Assuming that (5.2) is sufficiently accurate, we should be able to compare system (5.5) at $M=M_*$ with system (5.1) at $a_*=\mu M_*+a_{\infty }$ . At $M=M_*$ , the template solution $\overline {\xi }_{*}=\overline {\xi }(M_*)$ satisfies $U_{M_*}^{2^n}(\overline {\xi }_*)=\overline {\xi }_*$ . Substituting $x=\overline {x}(a_*)=\overline {x}_*=\gamma \overline {\xi }_*+X$ into (5.2) yields $f_{a_*}^{2^n}(\overline {x}_*)\approx \overline {x}_*$ , which implies that system (5.1) should exhibit a period- $2^n$ point well approximated by $x=\overline {x}_*$ at $a_*=\mu M_*+a_{\infty }$ . This orbit should be stable because differentiation of (5.2) using $\xi =(x-X)/\gamma$ demands that

(5.6) \begin{eqnarray} \left .\dfrac {{\textrm{d}}}{{\textrm{d}}x}\left (f^{2^{n}}_{a_*}\right )\right |_{x\,=\,\overline {x}_*}\approx \left .\dfrac {{\textrm{d}}}{{\textrm{d}}\xi }\left (U^{2^n}_{M_*}\right )\right |_{\xi \,=\,\overline {\xi }_*}, \end{eqnarray}

which tells that the multiplier of the $2^n$ orbit of the map $f_a$ must be reasonably well approximated by that of the universal map $U_M$ . This result is actually not surprising because (5.2) indicates that the dynamics induced by the map $f^{2^n}_a(x)$ around $x=X$ are quantitatively similar to those induced by $U_M^{2^n}(\xi )$ around $\xi =0$ . The same argument can be made for the correspondence between solutions $\overline {\xi }(M)$ and $\overline {x}(a)$ at any $a=\mu M+a_\infty$ , which need no longer be a single point as the central branch cascade unfolds. For $a\in (a_*,a_{\infty })$ , stable periodic orbits in the two systems can be related in the same way, and the approximate equality of multipliers follows from (5.6) analogously. Moreover, the dynamical similarity is not limited to periodic orbits, but applies also to aperiodic solutions for $a\gt a_\infty$ . In general, we must interpret that the template solution $\overline {\xi }=\overline {\xi }(M)$ consists of a set of points satisfying $U_{M}^{2^n}(\xi _i)\in \overline {\xi }$ for all $\xi _i \in \overline {\xi }$ (i.e. the set is invariant under $U_{M}^{2^n}$ ). Then, (5.2) implies that the set $\overline {x}(a)$ , defined as the collection of points $x_i=\gamma \xi _i+X$ for every $\xi _i \in \overline {\xi }(M)$ , is invariant under $f_{a}^{2^n}$ . Chaotic solutions related by $a=\mu M+a_\infty$ will have approximately the same Lyapunov exponents due to (5.6). All in all, the similarity of central branch dynamics between the two systems extends naturally all the way up to the accumulation point and beyond.

To find an adequate set of values for $\gamma , X, \mu ,$ and $a_{\infty }$ , one must pick two points in one of the bifurcation diagrams and find their analogues in the other. As a proof of concept, we select the first and second period-doubling bifurcation points that are encountered upon increasing $M$ beyond $M_*$ , where we have the stable period- $2^n$ orbit, namely PD $_{n+1}$ and PD $_{n+2}$ . We call the values of the parameter and of the central branch state variable $(M,\xi )=(M_{n+1},\hat {\xi }_{n+1})$ and $(M_{n+2},\hat {\xi }_{n+2})$ for PD $_{n+1}$ and PD $_{n+2}$ , respectively (see figure 12 c). These same period-doubling bifurcation points are found to occur at $(a,x)=(a_{n+1},\hat {x}_{n+1})$ and $(a_{n+2},\hat {x}_{n+2})$ for the dynamical system (5.1), which, in the context of our DNS computations, correspond to $(R,\tau )=(R_{n+1},\hat {\tau }_{n+1})$ and $(R_{n+2},\hat {\tau }_{n+2})$ of § 4. The matching of both points from one bifurcation diagram to the other is accomplished by enforcing $M_{m}= (a_{m}-a_{\infty })/\mu$ and $\hat {\xi }_{m}=(\hat {x}_{m}-X)/\gamma$ for $m=n+1$ and $n+2$ . Solving the four resulting algebraic equations, we get

(5.7) \begin{equation} \begin{array}{rcl} \gamma & = & \displaystyle {\frac {\hat {x}_{n+2}-\hat {x}_{n+1}}{\hat {\xi }_{n+2}-\hat {\xi }_{n+1}}},\\[1em] X & = & \displaystyle {\frac {\hat {\xi }_{n+2}\hat {x}_{n+1} -\hat {\xi }_{n+1} \hat {x}_{n+2}}{\hat {\xi }_{n+2}-\hat {\xi }_{n+1}}},\\[1em] \mu & = & \displaystyle {\frac {a_{n+2}-a_{n+1}}{M_{n+2}-M_{n+1}} }, \\[1em] a_{\infty } & = & \displaystyle {\frac {M_{n+2} a_{n+1} -M_{n+1}a_{n+2} }{M_{n+2}-M_{n+1}}}. \end{array} \end{equation}

If we set $n=2$ and exploit the DNS results of figure 10, the quadruplet of matching parameters takes the values $(\gamma , X,\mu , a_{\infty })\approx (-0.02518,0.2496,-0.4629,395.7213)$ . Finally, applying the transformation

(5.8) \begin{equation} \overline {x}(a)=\gamma \overline {\xi }(M)+X, \qquad a=\mu M+a_{\infty } \end{equation}

to the full cascade generated by system (5.5) results in the branch arrangement shown in figure 12(d). The transformed cascade exhibits a remarkably good agreement with DNS for the central branch (black dots). However, the quality of the match does not carry over to the rest of the template branches. It is precisely in this sense that the approximation (5.2) is to be considered only local.

The prediction can, however, be straightforwardly extended to all $2^n$ branches of the cascade undergone by the system (5.1), by simply allowing for a different set of scaling and translation parameters for each of the branches. The transformation of the $j$ th branch is performed according to

(5.9) \begin{equation} x_j^J(a)=\gamma _j \xi _j^J(M)+X_j, \qquad a=\mu M+a_{\infty }, \qquad j=1,2,3,\dots ,J, \end{equation}

where the pairs $(\gamma _j,X_j)$ are obtained by matching the bifurcation points in the respective bifurcation diagrams, following the same procedure as for the central branch. The prediction thus obtained for the complete unfolding of the bifurcation cascade undergone by the DNS is shown in figure 13. The agreement is now excellent for all four cascade branches. The success of the branch-by-branch matching process follows from the fact that the map $f^{2^n}_a(x)$ has $J=2^n$ extrema, the approximation (5.2) being applicable to each individually. Note that, in the chaotic regime after the accumulation point, the dynamics of the two systems is not expected to match in detail due to the sensitivity to initial conditions. Nonetheless, the statistical properties should, in principle, be in fair agreement, as evidenced by figure 13.

Figure 13. Branch-by-branch matching of period-doubling cascades. The black and grey points are the template branches shown in figure 12(b), re-adjusted according to the rescaling (5.9). The green points are the same DNS data as in figure 5(b). The matching process is based on the eight red points. The magenta points correspond to the P $_{12}$ orbit in the Taylor–Couette system.

For simplicity, we have limited our comparison across systems to orbits of the same period. It is also possible to relate stable orbits of different periods, as Feigenbaum (Reference Feigenbaum1982) did, to demonstrate the self-similarity of the bifurcation diagram (see Appendix). Furthermore, our theoretical result is consistent with the fact that self-similarity also holds beyond the accumulation point, and that the Lyapunov exponent exhibits the scaling law observed by Huberman & Rudnick (Reference Huberman and Rudnick1980). However, a detailed comparison of statistical quantities such as Lyapunov exponents is beyond the scope of our analysis, given that statistical convergence requires unaffordably long simulation times in the vicinity of the accumulation point (see Tirnakli, Tsallis & Beck Reference Tirnakli, Tsallis and Beck2009, for example). Instead, to gauge the accuracy of the prediction beyond the accumulation point, we have run DNS at a value of $R$ for which the occurrence of a periodic window with a P $_{12}$ orbit is anticipated. The time series of the Taylor–Couette system indeed appears to converge to a stable periodic orbit, which we have confirmed to be the Navier–Stokes counterpart of P $_{12}$ employing the PNK method (magenta dots in figure 13).

All self-similar features of the template period-doubling cascade should also show up in any cascade generated by a unimodal map. This implies that any other unimodal map, including the long-studied logistic map, could have been employed as a template. The advantage of the $U_M$ defined in (5.4) stems from its natural emergence in the derivation of the approximation theory (as seen in Appendix), which grants very rapid convergence to universal dynamics by construction. Extending the prediction method to more general cases is straightforward using the results by Briggs (Reference Briggs1991).

The core of our prediction method lies in completely separating the map-dependent scaling/translating parameters $(\mu , a_{\infty }, \gamma _j, X_j)$ from all universality-related considerations. A careful reading of Feigenbaum’s work suggests that he was likely aware of this possibility. Notably, the method of trajectory scaling functions was specifically designed to eliminate all ingredients that depend on any particular map. However, to the authors’ knowledge, there is no direct assertion in the literature pointing at the feasibility of such predictions, let alone their applicability beyond the accumulation point.