# Feigenbaum universality in subcritical Taylor-Couette flow

B. Wang1  R. Ayats1  K. Deguchi2<EMAIL_ADDRESS>A. Meseguer3  F.
Mellibovsky3<EMAIL_ADDRESS>1 Institute of Science and
Technology Austria (ISTA), 3400 Klosterneuburg, Austria 2 School of
Mathematics, Monash University, VIC 3800, Australia 3 Departament de Física,
Universitat Politècnica de Catalunya, 08034, Barcelona, Spain

###### Abstract

Feigenbaum universality is shown to occur in subcritical shear flows. Our
testing ground is the counter-rotation regime of the Taylor-Couette flow,
where numerical calculations are performed within a small periodic domain. The
accurate computation of up to the seventh period doubling bifurcation,
assisted by a purposely defined Poincaré section, has enabled the estimation
of the two Feigenbaum universal constants with unprecedented accuracy in a
fluid flow problem. We have further devised a method to predict the
bifurcation diagram up to the accumulation point of the cascade based on the
detailed inspection of just the first few period doubling bifurcations.
Remarkably, the method is applicable beyond the accumulation point, with
predictions remaining valid, in a statistical sense, for the chaotic dynamics
that follows.

##  1 Introduction

The counter-rotating regime of Taylor-Couette flow has long served as a
paradigm for subcritical transition to turbulence. The onset of fully
developed turbulence is preceded by intermittency, the flow typically
exhibiting alternating laminar and turbulent helical bands (Coles, 1965;
Andereck et al., 1986; Meseguer et al., 2009b; Dong, 2009) whose mean flow
topology (Wang et al., 2023) has recently been shown similar to the laminar-
turbulent stripe pattern observed for plane Couette flow (Barkley & Tuckerman,
2007). The multi-scale nature of the stripe structure complicates the detailed
understanding of subcritical transition with all the subtle features it
encompasses. Nevertheless, the aspects associated with the emergence of chaos
can often be approached in small periodic domains called minimal flow units
(Kawahara et al., 2012). A recent review paper (Feldmann et al., 2023)
hypothesised a connection between subcritical Taylor-Couette flow and plane
Couette flow, for which Kreilos & Eckhardt (2012) reported a period-doubling
cascade leading to chaos in a small periodic box. Here, we confirm that a
period-doubling route to chaos does indeed occur in Taylor-Couette flow,
within the parameter range studied by Meseguer et al. (2009b), Deguchi et al.
(2014), and Wang et al. (2022), hereafter abbreviated as W22. The minimal
periodic domain we use has an annular-parallelogram shape (Deguchi & Altmeyer,
2013) with two of its sides aligned along the helical band. As reported in
W22, the use of this domain was key to observing the first few period-doubling
bifurcations.

Our preliminary computations suggested that the subsequent period doubling
cascade may align exceptionally well with Feigenbaum universality (Feigenbaum,
1978, 1979, 1980, 1982), thus inspiring the detailed exploration we perform
here. This type of universality became widely recognised among fluid
dynamicists following the natural convection experiment conducted by Libchaber
et al. (1982), which led to the first estimation of Feigenbaum’s first
constant, $\delta_{\rm F}\approx 4.6692016$, corresponding to the limiting
ratio of a period-doubling bifurcation interval to the next (see also the
recent article by Libchaber, 2023). Their value, $\delta_{4}=4.4\pm 0.1$,
estimated from up to the fourth period-doubling bifurcation point (hence the
subindex $n=4$), was not too distant from the theoretical value $\delta_{\rm
F}$. Buzug et al. (1993) found an estimate close to $\delta_{\rm F}$ when
studying a Taylor-Couette apparatus with a tilted lid, while subsequent
experimental efforts, summarised in table 1, resulted in persistently larger
discrepancies.

System | Type |  $n$ | $\delta_{n}$ | Ref.
---|---|---|---|---
Rayleigh-Bénard (water) | E |  $4$ |  $4.3$ $\pm$ |  $0.8$ | Giglio et al. (1981)
Rayleigh-Bénard (mercury) | E |  $4$ |  $4.4$ $\pm$ |  $0.1$ | Libchaber et al. (1982)
Thermosolutal convection | N-2D |  $3$ |  $5$ $\pm$ |  $?$ | Moore et al. (1983)
Taylor-Couette (tilted lid) | E |  $3$ |  $4.67$ $\pm$ |  $?$ | Buzug et al. (1993)
Baffled channels | E & N-2D |  $3$ |  $4.65$ $\pm$ |  $?$ | Roberts & Mackley (1996)
Thermal convection | N-2D |  $2$ |  $4.9104$ $\pm$ |  $?$ | Lizée & Alexander (1997)
Triple periodic box | N-3D |  $2$ |  $4.88$ $\pm$ |  $0.5$ | Van Veen (2005)
Thermal convection | N-3D |  $4$ |  $4.321$ $\pm$ |  $?$ | Gao et al. (2015)
Forced oscillating cylinder | N-2D |  $5$ |  $4.52$ $\pm$ |  $?$ | Cheng et al. (2020)
MHD spherical Couette | N-3D |  $4$ |  $4.62$ $\pm$ |  $?$ | Garcia et al. (2021)
Table 1: Feigenbaum universality analyses in fluid systems. The table includes
the number of period doubling bifurcations analysed ($n$), and the
experimental (E) or numerical (N) nature of the study. Feigenbaum’s first
constant, as estimated from the last three period doubling bifurcations
analysed in each case is given in column $\delta_{n}$.

Estimates of Feigenbaum’s first universal constant have also been obtained
numerically for various fluid systems. Some of the values reported for
$\delta_{n}$, also listed in table 1, are reasonably close to $\delta_{\rm F}$
at first glance. However, it must be borne in mind that the reliability of the
estimates depends strongly on the number of consecutive period-doubling
bifurcations analysed and is extremely sensitive to the accuracy with which
they are located. Unfortunately, the estimates of $\delta_{\rm F}$ we have
found in the literature of fluid flows were not accompanied with a detailed
analysis as done, for example, for the Kuramoto-Sivashinsky equation (Smyrlis
& Papageorgiou, 1991). Accordingly, their reliability is debatable, as
concluding Feigenbaum universality from a handful of period-doubled solutions
is, at the very least, perilous.

All the studies presented in table 1 explore parameter ranges for which the
respective base flows are unstable, while our focus is set on the subcritical
regime. The aforementioned study by Kreilos & Eckhardt (2012) was the first to
show that period-doubling cascades may also occur in subcritical transition
problems. Their approach consisted in selecting the smallest possible periodic
domain that is capable of sustaining turbulence in plane Couette flow while,
at the same time, imposing specific discrete symmetries to further constrain
the dynamics. Symmetry restrictions were also employed by Moore et al. (1983)
and Van Veen (2005) in their respective works. A similar period-doubling
cascade in plane Couette flow has also been reported in a more recent paper by
Lustro et al. (2019). However, calculation of $\delta_{\rm F}$ in the
subcritical regime of fluid flow problems has not been conducted to date.

Although period-doubling cascades have been known to occur for over a century
(Collet, 2019), universality was not discovered until much later by Feigenbaum
(1978) and Coullet & Tresser (1978) independently, hence its being often
referred to as Feigenbaum-Coullet-Tresser universality by mathematicians. The
existence of universality swiftly disseminated throughout the mathematical
community, as vividly depicted by Khanin et al. (2021). Since its unveiling,
universality has been repeatedly observed in numerical simulations of low-
dimensional systems of ordinary differential equations. Checking universality
in low-dimensional models like the Rössler system or the Duffing equation has
become a common academic exercise. Feigenbaum universality was originally
established within the framework of discrete dynamical systems based on one-
dimensional iterated maps of the form $x_{\ell+1}=f(x_{\ell})$,
$\ell\in\mathbb{N}$, with $f(x)$ a unimodal function. Its necessary occurrence
in continuous time dynamical systems is typically justified by the use of
Poincaré sections. The phenomenon is often illustrated by initially replacing
the Smale horseshoe that occurs on the Poincaré section with a Hénon map,
which reduces to the logistic map in the limiting case of vanishing area
contraction rate after one mapping. Then, since the logistic map is a unimodal
map, universality of the period doubling cascade can be explained by means of
renormalisation theory, as done in many standard textbooks such as Collet &
Eckmann (1980); Schuster (1988); Glendinning (1994) or Strogatz (2024).
Universality manifests in both the parameter (coordinate) and state (ordinate)
axes of the bifurcation diagram, the scaling in the latter axis bearing close
relation with Feigenbaum’s second constant, $\alpha_{\rm F}\approx-2.5029079$.
In renormalisation theory, $\alpha_{\rm F}$ and the Feigenbaum function $G$
uniquely solve the so-called Feigenbaum-Cvitanović functional equation
$\mathcal{R}[G]=G$ (Feigenbaum, 1978; Cvitanović, 1989), with the condition
$G(0)=1$. Here, the action of the renormalisation operator, $\mathcal{R}$, on
a function $f$ consists merely in the twice repeated application of $f$,
mediated by a rescaling that involves a factor $\alpha_{\rm F}$:

$\displaystyle\mathcal{R}[f](x)=\alpha_{\rm
F}f\left(f\left(\frac{x}{\alpha_{\rm F}}\right)\right).$ (1)

Despite its simplicity, this operator lies at the heart of the various
instances of self-similarity that arise in the bifurcation diagram of the map
generated by $f$. Likewise, Feigenbaum’s first constant $\delta_{\rm F}$ is
the leading eigenvalue of the linearisation around $G$ of the renormalisation
operator $\mathcal{R}$, and $\varPhi$ is the associated eigenfunction (Lanford
III, 1982; Epstein, 1986; Eckmann & Wittwer, 1987; Briggs, 1991; Sthephenson &
Wang, 1991), which has recently been dubbed Feigenfunction by Thurlby (2021),
as a tribute to Feigenbaum. The first mathematical proof of the existence of
the universal constants, $\delta_{\rm F}$ and $\alpha_{\rm F}$, and functions,
$G$ and $\varPhi$, is attributed to Lanford III (1982).

Ascertaining universality from the estimation of $\delta_{\rm F}$ and
$\alpha_{\rm F}$ in fluid flow problems presents a twofold challenge. First, a
numerical analysis analogous to that of low-dimensional models demands
significantly higher computational resources when dealing with the Navier-
Stokes equations. The first estimation of $\delta_{\rm F}$ in the field of
computational fluid mechanics was contributed by Moore et al. (1983), when
analysing two-dimensional thermosolutal convection. It took, however, until
Van Veen (2005) to expose universality in a three dimensional setup, and even
then, only for a fluid contained in the simplest computational domain, i.e.,
periodic in all three space directions. Second, the governing equations are
very high-dimensional (formally infinite dimensional before spatial
discretisation) and thus the link to one-dimensional maps is far less obvious
than for low-dimensional dynamical systems. Feigenbaum universality does not
apply unconditionally to high-dimensional systems exhibiting period-doubling
cascades. Cascades sequentially doubling the period at every bifurcation may
occur through mechanisms quite different from those leading to universality
(e.g., Yalim et al., 2019). Other routes to chaos are also common in fluid
systems; in the Taylor-Couette system, for example, chaos has been found to
arise following the Ruelle-Takens-Newhouse scenario (Swinney & Gollub, 1985),
or Shil’nikov-type bifurcations (Lopez & Marques, 2005). Therefore, detecting
the onset of chaos, which requires a great deal of effort, does not
necessarily indicate that a universal period-doubling cascade has taken place.
All these difficulties justify why demonstrating Feigenbaum universality in a
fluid system requires careful analysis of the phase space. Buzug et al. (1993)
reconstructed the dynamics in phase space via delay-time embedding, while
Kreilos & Eckhardt (2012) and Gao et al. (2015) investigated the relationship
between consecutive local maxima of the time series of some physical quantity.
Even though both ways of sampling a continuous system to produce a discrete
system are valid, the use of a purposely defined Poincaré section provides the
most reliable approach.

The stable orbits of ever increasing period that arise in succession as the
parameter is varied along the period-doubling cascade accumulate at the so-
called accumulation point, beyond which chaotic dynamics ensue. Our interest
extends also to phenomena occurring past this point. The existence of a
reverse cascade, whereby chaotic bands successively merge in pairs as one
moves away from the accumulation point, is a well-established property of
simple model maps (Grossmann & Thomae, 1977; Lorenz, 1980). It is noteworthy
that Huberman & Rudnick (1980) observed self-similarity in the numerical
analysis of Lyapunov exponents at the onset of chaos. Nevertheless, the extent
to which universality holds in the chaotic regime remains an elusive question.
Libchaber (1983) tried to detect universality beyond the accumulation point in
a natural convection experiment. To the authors’ knowledge, this is the first
and only endeavour to unveil universality past the accumulation point of a
period-doubling cascade in a fluid dynamics problem, but the results were
rather crude due to experimental limitations. In theory, if universality holds
beyond the accumulation point, the bifurcation diagram should be predictable,
not only up to but also past this point, from the initial few period-doubling
bifurcations. However, we have not been able to find any such attempts in the
literature, even for low-dimensional models.

The paper is structured as follows. The problem formulation is presented in
§2, alongside a brief account of the numerical methods employed. §3 introduces
the period doubling cascade that is our object of study and illustrates the
procedure by which period doubling bifurcation points are accurately computed.
The sequence of the first few such points is then exploited in §4 to estimate
the first and second Feigenbaum constants and to extrapolate the expected
occurrence of the accumulation point. A method for the detailed prediction of
the bifurcation diagram in the neighbourhood of the accumulation point from
just a few initial period doubling bifurcations is devised in §5. We further
show how meaningful predictions can be made, in a statistical sense, that
extend beyond the accumulation point into the chaotic regime. Finally, the
main findings are summarised and conclusions drawn in §6.

##  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. The
angular velocities of the inner and outer cylinders, of radii $r^{*}_{\rm{i}}$
and $r^{*}_{\rm{o}}$, are denoted as $\Omega_{\rm{i}}$ and $\Omega_{\rm{o}}$,
respectively. A complete set of independent dimensionless parameters
characterising the problem are the radius ratio
$\eta=r^{*}_{\rm{i}}/r^{*}_{\rm{o}}$ and the two Reynolds numbers
$R_{i}=dr^{*}_{\rm{i}}\Omega_{\rm{i}}/\nu$ and
$R_{o}=dr^{*}_{\rm{o}}\Omega_{\rm{o}}/\nu$, where
$d=r^{*}_{\rm{o}}-r^{*}_{\rm{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 $\mbox{\bf
v}=(v_{r},v_{\theta},v_{z})=v_{r}\,\hat{{\bm{r}}}+v_{\theta}\,\hat{{\bm{\theta}}}+v_{z}\,\hat{{\bm{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.,

$\displaystyle\partial_{t}\mbox{\bf v}+(\mbox{\bf v}\cdot\nabla)\mbox{\bf v}$
$\displaystyle=$ $\displaystyle-\nabla p+\nabla^{2}\mbox{\bf
v}+f\,\hat{{\bm{z}}},$ (2) $\displaystyle\quad\nabla\cdot\mbox{\bf v}$
$\displaystyle=$ $\displaystyle 0,$ (3) $\displaystyle Q(\mbox{\bf
v})=\int_{0}^{2\pi}{\int_{r_{\rm{i}}}^{r_{\rm{o}}}{(\mbox{\bf
v}\cdot\hat{{\bm{z}}})\,r\,dr\,d\theta}}$ $\displaystyle=$ $\displaystyle 0,$
(4)

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

$\mbox{\bf v}(r_{i},\theta,z)=(0,R_{i},0),\quad\mbox{\bf
v}(r_{o},\theta,z)=(0,R_{o},0),$ (5)

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

$\mbox{\bf v}_{\rm b}=V_{\rm
b}\hat{{\bm{\theta}}}=\displaystyle\left(Ar+\frac{B}{r}\right)\,\hat{{\bm{\theta}}},\;\;p_{\rm
b}(r)=\int\frac{V_{\rm b}^{2}}{r}\,{\rm d}r,\;\;f_{\rm b}=0,$ (6)

with $A=(R_{o}-\eta R_{i})/(1+\eta)$ and $B=\eta(R_{i}-\eta
R_{o})/\left[(1-\eta)(1-\eta^{2})\right]$. In what follows we express the
velocity and pressure fields as

$\mbox{\bf v}=\mbox{\bf v}_{\rm b}(r)+\mathbf{u}(r,\theta,z;t),\quad p=p_{\rm
b}(r)+q(r,\theta,z;t).$ (7)

The fields $q$ and
$\mathbf{u}=u_{r}\,\hat{{\bm{r}}}+u_{\theta}\,\hat{{\bm{\theta}}}+u_{z}\,\hat{{\bm{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 in Meseguer et al. (2007). 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\left[r_{\rm{i}},r_{\rm{o}}\right]\times\left[0,2\pi\right]\times\left[0,2\pi\right]$,
where the $2\pi$-periodicity is imposed in the transformed coordinates

$\xi=n_{1}\theta+k_{1}z,\quad\zeta=n_{2}\theta+k_{2}z,$ (8)

following Deguchi & Altmeyer (2013). Figure 1a 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)$.

(a) | (b)
---|---
|
Figure 1: (a) Annular-parallelogram computational domain defined by the
coordinates of (8) with wave numbers $(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.

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.

In order 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. (2020) 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_{\rm{i}}$ and $\tau_{\rm{o}}$,

$\displaystyle\kappa=\frac{E(\mbox{\bf u})}{E(\mbox{\bf
v}_{b})},\qquad\tau_{\rm i,\rm o}=\left.1+\frac{\partial_{r}(r^{-1}\langle
u_{\theta}\rangle_{\xi\zeta})}{\partial_{r}(r^{-1}V_{b})}\right|_{r=r_{\rm
i},r_{\rm o}},$ (9)

where

$E(\mbox{\bf v})=\frac{1}{2\mathcal{V}}\iiint_{\mathcal{V}}{\mbox{\bf
v}\cdot\mbox{\bf
v}\,d\mathcal{V}}=\frac{1}{2\mathcal{V}}\int_{0}^{2\pi}\int_{0}^{2\pi}\int_{r_{\rm{i}}}^{r_{\rm{o}}}\mbox{\bf
v}\cdot\mbox{\bf
v}\,r\,\mathrm{d}r\mathrm{d}\xi\mathrm{d}\zeta=\frac{1-\eta}{1+\eta}\int_{r_{\rm{i}}}^{r_{\rm{o}}}{\langle\mbox{\bf
v}\cdot\mbox{\bf v}\rangle_{\xi\zeta}\,r\,\mathrm{d}r}$ (10)

is the volume-averaged kinetic energy of some velocity field v. The volume of
the transformed computational domain is
$\mathcal{V}=2\pi^{2}(r_{\rm{o}}^{2}-r_{\rm{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_{\rm{i}}=\tau_{\rm{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
same velocity field represent the same solution. We have systematically
factored out the group orbit invariance employing the method of slices
(Budanur et al., 2015), using the same template as Wang et al. (2023). To
analyse the period doubling cascade within the framework of discrete-time
dynamical systems, we have devised a Poincaré section

$\Sigma=\left\\{{\bf a}\in\mathbb{X}\left|\;\tau_{i}({\bf a})=\tau_{o}({\bf
a}),\;\dfrac{d\tau_{i}}{dt}>\dfrac{d\tau_{o}}{dt}\right.\right\\},$ (11)

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_{\rm{i}}=\tau_{\rm{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, inner or outer
indistinctly, on $\Sigma$.

##  3 The period doubling cascade

Emulating the influential experiments by Andereck et al. (1986), 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 full-
fledged turbulence. A reduction of $R_{i}$ while still in the supercritical
region has the flow re-organise into a pattern of alternated laminar and
turbulent winding helical bands, commonly known as the spiral turbulence
regime (Coles, 1965; Andereck et al., 1986; Meseguer et al., 2009b; Dong,
2009). 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 persists below this threshold (Meseguer et al., 2009a; Deguchi et al.,
2014, 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 1a to seek for
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 1b) that coexists with the stable CCF. This drifting rotating wave,
DRW for short, embodies all the essential elements of the self-sustaining
mechanism (Wang et al., 2007; Hall & Sherwin, 2010) and plays a central role
in the transition process. In W22, two periodic solutions, P1, arising from a
Hopf bifurcation of DRW, and P2, originated at a period-doubling bifurcation
of P1, 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 The 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 subindex and simply
write $R$ to ease notation.

Figure 2a shows the DRW (square), P1 (dashed blue line) and P2 (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)$.

(a) (b) (c)

Figure 2: The Poincaré section $\Sigma$ and the periodic orbits P1 (dashed
blue) and P2 (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 P1 and P2. (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 P1 (empty blue) and P2 (filled
green) curves correspond to their representation on $\Sigma$.

The Poincaré section $\Sigma$ appears in this representation as the
transparent grey plane, which contains DRW and is pierced, in the direction
defined by (11), at a single point by P1 and at two different points by P2. In
the torque time series of figure 2b, 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 2a onto the plane $(\tau_{o},\tau_{i})$,
as depicted in figure 2c, further clarifies how the sequence of intersection
points of figure 2b collapses in phase space onto a single point for P1 (empty
blue disk) and as two distinct points for P2 (filled green circles), all
contained in $\Sigma$ (dashed grey straight line).

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 3ai
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 1a).

(ai) |  | (bi) |
---|---|---|---
(aii) |  | (bii) |
(ciii) |  | (biii) |
Figure 3: Space-time diagrams for (a) DRW, (b) P1 and (c) P2, 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.

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 P1, as
shown in figure 3bi. P1 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 dismissed by attaching the line probe to a moving
reference frame defined with the method of slices (Budanur et al., 2015; Wang
et al., 2023). In this reference frame, the space-time diagrams of DRW and P1
(figures 3aii and 3bii) are greatly simplified, the former appearing as time-
independent and the latter as purely time-periodic.

Further subtraction of the time average helps expose the true nature of the
time dependence of the solutions. For P1, the time-horizon considered in
figure 3biii allows for just over four repetitions, as indicated by the dashed
vertical lines. Solution P2, shown in figure 3ciii, has instead a natural
period about twice that of P1, 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 4a shows the bifurcation sequence that DRW undergoes when varying
the inner Reynolds number in the range $R\in[391.20,395.90]$.

(a) | (b)
---|---
|
Figure 4: Bifurcation scenario as recorded on the Poincaré section $\Sigma$.
(a) The initial steps of the bifurcation scenario. Shown are DRW (black), P1
(blue) and P2 (green), reported in W22. P4 (green) emerges at the second
period doubling bifurcation point PD2. 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.

DRW emerges from a saddle-node bifurcation (SN1) at $R\simeq 391.50$ and,
leaving aside subharmonic instabilities, the nodal (upper) branch (the one
represented in figure 1b and the square in figures 2a,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 subcritically stable non-trivial solution in 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, of course, be used to compute
unstable solution branches. A stable relative periodic orbit (P1, 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 P1 branch becomes unstable in a period doubling bifurcation
(PD1), whence a branch of stable period-doubled relative periodic orbits (P2,
pair of green lines) emerges. The P2 loses stability, in turn, in a second
period-doubling bifurcation (PD2) issuing a branch of period-4 solutions. The
$R$ value for which figures 2 and 3 were computed corresponds to just short of
PD2, hence the instability of DRW and P1, in contrast with P2, which is
stable.

The remainder of the period-doubling cascade and the onset of chaotic dynamics
at larger values of $R$ is depicted in figure 4b. Only stable solutions are
shown, so DNS is all that has been required to produce the diagram. Estimating
$\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 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. Both complications combined render the careful
estimation of 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 initial
transients have been overcome. 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 5a shows $\tau(\ell)$
for the P4 solution (red circles) at $R=395.711$ (between PD2 and PD3) and for
the P8 solution (blue) at $R=395.719$ (between PD3 and PD4).

(a) (b)
---
Figure 5: All orbits up to period 8 at (P1 and P2, unstable) and around (P4
and P8, at $R=395.711$ and $395.719$, respectively) PD3. (a) Torque
($\tau=\tau_{i}=\tau_{o}$) of P1 (cyan), P2 (green), P4 (blue) and P8 (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.

Both solutions are stable and well converged with DNS, such that the
respective series sequentially repeat, always in the same order, the same four
(P4) and eight (P8) distinct values, as indicated by the numbered horizontal
dashed lines. The effect of the period-doubling bifurcation PD3 can thus be
portrayed as doubling point $j$ of the P4 discrete-time orbit into points $j$
and $j+4$ of the P8 cycle, the latter two points having sprung from the former
and drifted away in opposite directions. The unstable P1 (cyan) and P2 (green)
orbits in the immediate vicinity of PD3 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 2c is shown in figure 5b, but with the
region where $\Sigma$ is pierced by the orbits duly magnified to further
clarify where the P4 and P8 cycles stand in relation to P1 and P2.

The $n$th period doubling bifurcation PDn can be readily analysed by splitting
the $\tau(\ell)$ sequences of interest in $J=2^{n}$ separate subsequences

$\tau_{j}^{J}(k)=\tau(kJ+j),$ (12)

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

$(\tau_{j}^{J})_{\infty}=\lim_{k\to\infty}\tau_{j}^{J}(k)\qquad
j\in\\{1,\ldots J\\}$ (13)

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

$A_{j}^{J}=\left\lvert(\tau_{j}^{J})_{\infty}-(\tau_{j+J/2}^{J})_{\infty}\right\rvert\qquad
j\in\\{1\ldots J/2\\},$ (14)

will all vanish for any orbit other than PJ. For instance, suppose we slightly
perturb the stable P8 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 5). Applying the same $J=8$ sampling to a stable P4
(or P2 or P1) will instead produce indistinguishable limits for $j$ and $j+4$
and, consequently, vanishing amplitudes $A_{j}^{8}=0$.

To nail down PD3, a collection of DNS runs traversing the bifurcation point,
i.e., between $R=395.711$ and $395.719$, are required. The amplitude of P8
drops fast as the parameter $R$ is decreased towards PD3, such that discerning
it from P4 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 6a and 6b
depict the torque subsequences at $R=395.7116$ and $R=395.7122$, respectively,
below and above PD3 but moderately close to the bifurcation point.

(a) | (b)
---|---
|
Figure 6: Analysis of torque sequences $\tau(\ell)$, sampled with $N=8$ at (a)
$R=395.7116$ and (b) $R=395.7122$, below and above PD3 ($R=395.711841564$, as
we shall see later), respectively. The top panels show the 8 sequences
$\tau_{j}^{8}(k)$ (circles, magnified in the insets), alongside respective
power law fits (dashed lines). To differentiate the 8 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
$((\tau_{j}^{8}(k))_{\infty}+(\tau_{j+4}^{8}(k))_{\infty})/2$ subtracted from
every pair of branches to illustrate convergence. The bottom panels show again
the same data but in logarithmic scale.

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

$\tau_{j}^{J}(k)=(\tau_{j}^{J})_{\infty}+\left((\tau_{j}^{J})_{0}-(\tau_{j}^{J})_{\infty}\right)(\lambda_{j}^{J})^{k},$
(15)

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 6 are the fits to the
sequences $\tau_{j}^{J}(k)$, $k\geq 0$. In order to expose the convergence
rate and the accuracy of the fits, the same data has 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 have 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 P4 and the latter on a P8. 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 P8 regardless
of its actual period.

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 PD3 from either side, and plotted in figure
7.

(a) (b) |
---|---
Figure 7: Analysis of PD3. (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 6, 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.

Figure 7a 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 PDn, a fit of the form

$A_{j}^{J}(R)=a_{j}^{J}\sqrt{R-R_{j}^{J}}$ (16)

with $J=2^{n}$ is 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 7a (dashed lines), yields $R_{3}=\langle
R_{j}^{8}\rangle=395.71184156\pm 5\times 10^{-8}$ (violet square) for PD3.

The eight $\lambda_{j}^{8}$ multipliers, shown in figure 7b, behave as
expected across a period doubling bifurcation. The eight individual fits (15)
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 PD3 is indeed a P8, $\lambda_{8}$ is directly its multiplier.
This is not so for the P4 solution below PD3. Plots analogous to those in
figure 6a but sampling with $J=4$ show that convergence onto the four
asymptotic branches follows an oscillatory pattern. The multiplier of the P4
orbit in the vicinity of PD3 must be either estimated 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 P4 in the vicinity of PD3 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 7b),
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 8, where all orbits computed up to period
$N=128$ (past PD7) are shown as a function of $R$.

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

Successive magnifications, presented in 8b and 8c, exhibit branch structures
very similar to that of the full cascade of 8a. Yet, at this stage, there is
no guarantee that the universality observed here is of the same kind as
discovered by Feigenbaum.

In order to ascertain Feigenbaum universality, we have estimated 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.

$n$ | $N=2^{n}$ | $R_{n}$ | $\delta_{n}$ | $\alpha_{n}$
---|---|---|---|---
1 | 2 | 395. | 4762796241073 |  $\pm 6.9$$\times$ | $10^{-13}$ | - | - | - | -
2 | 4 | 395. | 67542342 |  $\pm 5.3$$\times$ | $10^{-7}$ | - | - | - | -
3 | 8 | 395. | 711841564 |  $\pm 4.2$$\times$ | $10^{-8}$ | 5. | 46826 |  $\pm 1.0$$\times$ | $10^{-4}$ | -2. | 6668339 |  $\pm 8.1$$\times$ | $10^{-6}$
4 | 16 | 395. | 719257726 |  $\pm 8.7$$\times$ | $10^{-8}$ | 4. | 91065 |  $\pm 1.6$$\times$ | $10^{-4}$ | -2. | 611741 |  $\pm 1.6$$\times$ | $10^{-5}$
5 | 32 | 395. | 72083697363 |  $\pm 2.9$$\times$ | $10^{-10}$ | 4. | 69601 |  $\pm 3.4$$\times$ | $10^{-4}$ | -2. | 496646 |  $\pm 2.9$$\times$ | $10^{-5}$
6 | 64 | 395. | 7211746289 |  $\pm 4.7$$\times$ | $10^{-9}$ | 4. | 67710 |  $\pm 3.3$$\times$ | $10^{-4}$ | -2. | 510893 |  $\pm 2.8$$\times$ | $10^{-5}$
7 | 128 | 395. | 7212469242 |  $\pm 1.4$$\times$ | $10^{-9}$ | 4. | 67050 |  $\pm 4.6$$\times$ | $10^{-4}$ | -2. | 500577 |  $\pm 4.4$$\times$ | $10^{-5}$
$\dots$ | $\dots$ | $\dots$ | $\dots$ | $\dots$ |  | $\dots$ |  |
$\infty$ | $\infty$ | 395. | 721266624 |  $\pm 1.1$$\times$ | $10^{-8}$ | 4.66920160910299… | -2.50290787509589…
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 PDn are used to compute the $n$th
estimate to Feigenbaum’s first constant, $\delta_{n}$, according to (17). An
estimate to the second constant, $\alpha_{n}$, is computed from central branch
torque values following (24). 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 (15), as conducting a
systematic study is difficult). The last row of the table corresponds to the
accumulation point $R_{\infty}$ estimated by (23), and the actual values of
Feigenbaum constants, $\delta_{\rm F}$ and $\alpha_{\rm F}$.

For every $n\geq 3$, the table also lists the ratios $\delta_{n}$, defined as

$\delta_{n}=\dfrac{R_{n-1}-R_{n-2}}{R_{n}-R_{n-1}},$ (17)

which provide the best estimate of Feigenbaum’s first constant $\delta_{\rm
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_{\rm 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.

The $R_{n}$ sequence is expected to converge onto the accumulation point,
$R=R_{\infty}$. This point may be estimated by assuming that all period
doubling bifurcations of order higher than available are perfectly Feigenbaum-
universal and, therefore, satisfy (17) with the left hand side replaced by
$\delta_{\rm 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

$\displaystyle R_{\infty}=R_{n}+\sum_{m=n+1}^{\infty}{\Delta R_{m}},$ (18)

where $\Delta R_{m}\equiv R_{m}-R_{m-1}$ denotes the interval between the
$(m-1)$ and $m$th period-doubling bifurcations. Since (17) 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

$\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}.$ (19)

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

$\displaystyle R_{\infty}=R_{n}+\Delta
R_{n}\left[\sum_{m=n+1}^{\infty}{\prod_{\ell=n+1}^{m}\delta_{\ell}^{-1}}\right].$
(20)

If $\delta_{\ell}=\tilde{\delta}_{n}>1$ were constant for all $\ell>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
PDn as

$\tilde{\delta}_{n}=1+\left[\sum_{m=n+1}^{\infty}\displaystyle{\prod_{\ell=n+1}^{m}{\delta_{\ell}}^{-1}}\right]^{-1},$
(21)

so that (20) simplifies to

$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}.$
(22)

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

$R_{\infty}\in\left[R_{n}+\dfrac{R_{n}-R_{n-1}}{\delta_{n}-1},R_{n}+\dfrac{R_{n}-R_{n-1}}{\delta_{\rm
F}-1}\right],$ (23)

by substituting $\tilde{\delta}_{n}=\delta_{n}$ and
$\tilde{\delta}_{n}=\delta_{\rm F}$ in (22), respectively. Using PD7, 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 9: Self-similarity of the period doubling cascade. Same data as figure
8 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 labeled PDn
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 9). 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 10a depicts a three-dimensional phase map projection of the orbit at
the accumulation point $R_{\infty}$.

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

If the manifold has dimension one, the coefficients $\mathbf{a}$ on the
Poincaré section $\Sigma$ introduced in (11) can be uniquely parameterised 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 10b). The
structure of the map on $\Sigma$ is, however, not clear, as the
$\Sigma$-crossings of the orbit are sparsely distributed (see figure 10a).
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$. The
function $f$ is 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, 1991).

For period-doubling cascades in one-dimensional unimodal maps, the second
Feigenbaum constant, $\alpha_{\rm F}$, is often estimated 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 (1979) 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, the estimation of $\alpha_{\rm F}$
has not been conducted.

On the other hand, we have a highly accurate bifurcation diagram available,
making the direct estimation of the second constant possible. In the torque
sequences presented so far, we have indexed the $N$ distinct points visited by
every periodic solution PN as $j=1,2,\dots N$ such that $j=1$ is the central
branch point (the black dots in figure 9). Thus, using the torque values of
the central branch at PDm, $\hat{\tau}_{m}=\tau_{1}^{2^{m}}(R_{m})$, we can
compute the ratios

$\displaystyle\alpha_{n}=\dfrac{\hat{\tau}_{n-1}-\hat{\tau}_{n-2}}{\hat{\tau}_{n}-\hat{\tau}_{n-1}}$
(24)

to provide estimates for the $\alpha_{\rm 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}$, evince
gradual convergence toward Feigenbaum universality.

##  5 Prediction of the period doubling route to chaos

Following the seminal work of Feigenbaum (1979), who showcased the first
comparison with experimental results, Feigenbaum (1980, 1982, 1983) 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 so-called trajectory scaling function theory, which unravels
the scaling between orbits featuring the same stability properties (i.e.,
having the same value of the multiplier), was later validated against
experimental data on natural convection by Belmonte et al. (1988).

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 method of trajectory scaling function 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

$\displaystyle x_{\ell+1}=f_{a}(x_{\ell}).$ (25)

We assume that the function $f_{a}(x)$ is unimodal, and that a fixed point of
the system (25) 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 (1982), is that the approximation

$\displaystyle\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}$ (26)

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

$\displaystyle\left|\frac{\delta^{n}(a-a_{\infty})}{\mu}\right|\ll
1,\qquad\left|\frac{\alpha^{n}(x-X)}{\gamma}\right|<1,$ (27)

and provided the constants $\gamma$, $X$, $\mu$, and $a_{\infty}$ are suitably
chosen. The derivation of (26) and (27) is given in Appendix A. 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
subindex when referring to Feigenbaum’s first ($\delta$) and second ($\alpha$)
constants.

Here, the function $U_{M}$ appearing in the approximation (26) is defined in
terms of the Feigenbaum function $G$ and the Feigenfunction $\varPhi$ (with
the normalisation $\varPhi(0)=1$) as

$\displaystyle U_{M}(\xi)=G(\xi)+M\varPhi(\xi).$ (28)

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 A. The nature of $G$ and
$\varPhi$ is well understood (Collet & Eckmann, 1980; Briggs, 1991; Thurlby,
2021). 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.

If the four (map-dependent) parameters, $\gamma$, $X$, $\mu$, and
$a_{\infty}$, can be estimated from the first few period doublings, then (26)
provides a useful tool for predicting the behaviour of the dynamical system
(25), subject to the constraints outlined by (27). In practice, fixing $n=2$
and using data up to PDn+2 (i.e. PD4) 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

$\displaystyle\xi_{\ell+1}=U_{M}(\xi_{\ell}),$ (29)

and then stretch and translate them to align with corresponding branches of
the bifurcation diagram of (25). The constants $X$ and $\gamma$ translate and
scale the state $x$, while $a_{\infty}$ and $\mu$ do the same for the
parameter $a$, respectively. The constant $a_{\infty}$ corresponds to the
predicted accumulation point of the dynamical system (25). Thus the first
condition in (27) 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 (27) 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 (27) 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 branches.

The bifurcation diagram of the dynamical system (29) is depicted in figure
11a.

(a) | (b)
---|---
|
(c) | (d)
|
Figure 11:  Bifurcation diagram of the one-dimensional map (29). (a) The full
period doubling cascade. (b) Magnification around the accumulation point
$M_{\infty}$. The $2^{n}=4$ distinct sets of points arising for $M>M_{*}$ are
labeled 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 PD3 and PD4 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) The full cascade transformed by (32) with
$(\gamma,X,\mu,a_{\infty})$ obtained by (31), compared with the period
doubling cascade in Taylor-Couette flow (green dots).

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 11b, which contains a magnification of figure 11a
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 labeled the four cascades as $\xi_{j}^{4}$ ($j=1,2,3,4$) in figure 11b
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 11c).

To illustrate how the approximation (26) operates, we begin by fixing
$M=M_{*}$ in figure 11b, 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
(12) with $J=2^{2}$. Assuming that (26) is sufficiently accurate, we should be
able to compare system (29) at $M=M_{*}$ with system (25) 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 (26)
yields $f_{a_{*}}^{2^{n}}(\overline{x}_{*})\approx\overline{x}_{*}$, which
implies that system (25) 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 (26) using
$\xi=(x-X)/\gamma$ demands that

$\displaystyle\left.\dfrac{d}{dx}\left(f^{2^{n}}_{a_{*}}\right)\right|_{x=\overline{x}_{*}}\approx\left.\dfrac{d}{d\xi}\left(U^{2^{n}}_{M_{*}}\right)\right|_{\xi=\overline{\xi}_{*}},$
(30)

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 (26) 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 (30) analogously. Moreover, the dynamical similarity is not limited to
periodic orbits but applies also to aperiodic solutions for $a>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 (26) 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 (30). 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.

In order 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 PDn+1 and PDn+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 PDn+1 and PDn+2, respectively (see figure
11c). 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 (25), 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

$\begin{array}[]{rcl}\gamma&=&\displaystyle{\frac{\hat{x}_{n+2}-\hat{x}_{n+1}}{\hat{\xi}_{n+2}-\hat{\xi}_{n+1}}},\\\\[10.00002pt]
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}}},\\\\[10.00002pt]
\mu&=&\displaystyle{\frac{a_{n+2}-a_{n+1}}{M_{n+2}-M_{n+1}}},\\\\[10.00002pt]
a_{\infty}&=&\displaystyle{\frac{M_{n+2}a_{n+1}-M_{n+1}a_{n+2}}{M_{n+2}-M_{n+1}}}.\end{array}$
(31)

If we set $n=2$ and exploit the DNS results of figure 9, 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

$\overline{x}(a)=\gamma\overline{\xi}(M)+X,\qquad a=\mu M+a_{\infty}$ (32)

to the full cascade generated by system (29) results in the branch arrangement
shown in figure 11d. 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 template branches. It is
precisely in this sense that the approximation (26) 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 (25), 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

$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,$ (33)

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
12.

Figure 12: Branch-by-branch matching of period doubling cascades. The black
and grey points are the template branches shown in figure 11b, re-adjusted
according to the rescaling (33). The green points are the same DNS data as in
figure 4b. The matching process is based on the 8 red points. The magenta
points correspond to the P12 orbit in the Taylor-Couette system.

The agreement is now excellent for all 4 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 (26) being
applicable to each individually. Note that, in the chaotic regime after the
accumulation point, the dynamics of the two systems are 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 evinced
by figure 12.

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 (1982) did, to demonstrate the self-similarity of the
bifurcation diagram (see Appendix A). 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 (1980). 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 et
al. (2009), 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 P12 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
P12 employing the PNK method (magenta dots in figure 12).

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 template. The advantage of the $U_{M}$ defined in (28) stems from
its natural emergence in the derivation of the approximation theory (as seen
in Appendix A), 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, 1991).

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.

##  6 Conclusions

We have studied in detail a period-doubling cascade that arises in subcritical
counter-rotating Taylor-Couette flow employing small computational domains of
annular-parallelogram shape. The cascade is seeded from a family of drifting
rotating waves discovered by W22 and eventually leads to a chaotic regime
fueled by the self-sustained process of wall-bounded shear flows.

A pivotal outcome of our analysis is the confirmation of Feigenbaum
universality to an unprecedented degree of accuracy in the context of fluid
dynamics. Our success is of course largely ascribable to the fast development
of computational power over the past few decades, along with the outstanding
numerical accuracy of the methods we employ to solve Navier-Stokes flows,
which altogether has rendered thorough parameter sweeps feasible. That said,
establishing universality has required very long and costly numerical
simulations, a refined parametric exploration, and the deployment of very
accurate numerical techniques for the detailed analysis of the asymptotic
convergence of relative periodic orbits, the assessment of their stability and
the estimation of the period-doubling bifurcation points. Determining both the
first and second Feigenbaum constants to the third significant digit demands
the precise computation of up to the seventh period-doubling bifurcation along
the cascade. This has enabled an accurate estimation of the accumulation point
and, with it, the unambiguous exposure of the self-similar structure of the
period-doubling bifurcation cascade. A key element to make the analysis
systematic has been the deployment of a convenient Poincaré section based on
torque balance, which provides a robust sampling method to validate theories
of unimodal discrete-time dynamical systems.

Furthermore, our results constitute the first confirmation of universality in
a fluid flow problem subject to subcritical transition. It is also noteworthy
that all solutions along the period doubling cascade exhibit spatial drift, a
property that is unique to our system among the many previous studies on the
subject (see table 1). As expected, the conventional Feigenbaum theory becomes
applicable to the Taylor-Couette phase space once reduced by the method of
slices (Budanur et al., 2015).

Drawing from Feigenbaum’s theory, we have further developed a method to
predict the bifurcation diagram of the period-doubling cascade all the way up
to the accumulation point and, remarkably, also beyond. This approach has
broader applicability than the method of trajectory scaling functions
(Feigenbaum, 1983; Belmonte et al., 1988), which is limited to scrutinising
stable orbits with the same value of the multiplier. Our prediction method is
simple: it only requires stretching and shifting a template bifurcation
diagram produced by any smooth unimodal map at hand. Of course, the scaling
and translating parameters must be chosen such that the template bifurcation
diagram aligns with the period-doubling cascade of the target system (in our
case, Taylor-Couette flow). This can be accomplished, for example, by matching
two bifurcation points across systems. The remainder of the transformed
template diagram provides the prediction for the target system. The prediction
results are useful, among other things, to anticipate the location of periodic
windows in the bifurcation diagram of the target system beyond the
accumulation point. The emergence of a stable orbit in the Taylor-Couette
system precisely within the expected parameter range attests to the
exceptional predictive capability of the method.

Although the Navier-Stokes equations should formally be considered an
infinite-dimensional dynamical system, the dynamics are often confined to
finite-dimensional manifolds in phase space, sometimes designated as inertial
manifolds (Temam, 1989). This low dimensionality has recently attracted
significant attention (Ding et al., 2016; Haller et al., 2023), as it affords
valuable insight into the underlying dynamics of the flow, sometimes leading
to simplified models whose numerical simulation is computationally more
tractable. Our results can be interpreted as an extreme case of simplification
that ultimately leads to a one-dimensional map representations.

The observation of Feigenbaum universality implies that the dynamics around
the accumulation point of our Taylor-Couette setup can indeed be approximated
by a nearly-one-dimensional discrete map when analysed on an appropriate
Poincaré section. Hence, our finding provides a handy playground for testing
in fluid dynamics the application of iconic theorems from the early days of
chaos theory, such as Li-Yorke or Sharkovsky’s (Sharkovskii, 1964; Li & Yorke,
1975), as well as cycle expansion theory (Cvitanović et al., 2016). It is our
intention to explore the significance and implications of these theories to
the field of fluid dynamics in the near future.

BW’s and RA’s research have been both funded by the European Union’s Horizon
2020 research and innovation programme under the Marie Skłodowska-Curie Grant
Agreement No. 101034413. This research was supported by the Australian
Research Council Discovery Project DP230102188 and the Spanish MICINN grant
No. PID2020-114043GB-I00.

The authors report no conflict of interest.

##  Appendix A Derivation of a universal approximation to a unimodal map and
the conditions for validity

Consider the one-dimensional map (25), with $f_{a}(x)$ a unimodal function,
undergoing a period-doubling cascade as the parameter $a$ is increased. The
location where the function has a maximum for $a=a_{\infty}$, the accumulation
point, is denoted as $X=\text{argmax}f_{a_{\infty}}(x)$. It will be convenient
to shift the coordinate according to $y=x-X$ and define the shifted map as
$F_{a}(y)=f_{a}(y+X)-X$.

If $a$ is close to the accumulation point $a_{\infty}$, the function $F_{a}$
may be approximated by a Taylor expansion truncated at first order, i.e.

$\displaystyle F_{a}(y)\approx\varPsi(y)+(a-a_{\infty})\psi(y),$ (34)

where $\varPsi=F_{a}|_{a=a_{\infty}}$ and
$\psi=\left.\partial_{a}F_{a}\right|_{a=a_{\infty}}$. Unimodal maps are known
to be infinitely renormalisable at accumulation points (see Lanford III, 1982,
for example). Therefore, for large $n$ and $a=a_{\infty}$, one must be able to
find a function $g$, satisfying $\mathcal{R}[g]=g$, such that

$\displaystyle\mathcal{R}^{n}[\varPsi](y)\approx g(y).$ (35)

Note that $\gamma=g(0)$ is not necessarily unity and depends on the map
$f_{a}$. The function $g$ and the Feigenbaum function $G$ introduced in §1 are
related by

$\displaystyle G(\xi)=g(\gamma\xi)/\gamma.$ (36)

Next we consider the Fréchet derivative of $\mathcal{R}$, in the direction of
function $h(y)$, and evaluated at function $f(y)$:

$\displaystyle\mathcal{L}_{f}[h](y)$ $\displaystyle=$
$\displaystyle\lim_{\epsilon\rightarrow 0}\frac{\mathcal{R}[f+\epsilon
h](y)-\mathcal{R}[f](y)}{\epsilon}=\alpha
h\left(f\left(\frac{y}{\alpha}\right)\right)+\alpha
f^{\prime}\left(f\left(\frac{y}{\alpha}\right)\right)h\left(\frac{y}{\alpha}\right).$
(37)

The operative expression of $\mathcal{L}_{f}[h](y)$ in terms of $f$ and $h$ is
a well-known result (see Briggs, 1991; Thurlby, 2021, for example). The
definition (37), combined with the expansion (34), implies
$\mathcal{L}_{\varPsi}[\psi]\approx(\mathcal{R}[F_{a}]-\mathcal{R}[\varPsi])/(a-a_{\infty})$,
such that simple algebraic manipulation yields

$\displaystyle\mathcal{R}[F_{a}]\approx\mathcal{R}[\varPsi]+(a-a_{\infty})\mathcal{L}_{\varPsi}[\psi].$
(38)

Moreover, it can be shown by induction that the $n$th renormalisation of
$F_{a}$ is given by

$\displaystyle\mathcal{R}^{n}[F_{a}]\approx\mathcal{R}^{n}[\varPsi]+(a-a_{\infty})\mathcal{J}_{n}[\psi],$
(39)

where
$\mathcal{J}_{n}=\mathcal{L}_{\mathcal{R}^{n-1}[\varPsi]}\circ\mathcal{L}_{\mathcal{R}^{n-2}[\varPsi]}\circ\dots\circ\mathcal{L}_{\mathcal{R}^{0}[\varPsi]}$
and, of course, $\mathcal{R}^{0}[\varPsi]=\varPsi$. To deduce (39), let us
define $Q_{k}=\mathcal{R}^{k}[\varPsi]+(a-a_{\infty})\mathcal{J}_{k}[\psi]$.
For $n=1$ the statement (39) reduces to $\mathcal{R}[F_{a}]\approx Q_{1}$,
which is trivially satisfied in view of (38). Now, assuming that
$\mathcal{R}^{k}[F_{a}]\approx Q_{k}$ holds at step $k$, we just need to prove
that $\mathcal{R}^{k+1}[F_{a}]\approx Q_{k+1}$ also holds. Using (37), we get

$\displaystyle\mathcal{J}_{k+1}[\psi]=\mathcal{L}_{\mathcal{R}^{k}[\varPsi]}\left[\mathcal{J}_{k}[\psi]\right]\approx\frac{\mathcal{R}\left[\mathcal{R}^{k}[\varPsi]+(a-a_{\infty})\mathcal{J}_{k}[\psi]\right]-\mathcal{R}\left[\mathcal{R}^{k}[\varPsi]\right]}{a-a_{\infty}},$
(40)

which implies
$\mathcal{R}[Q_{k}]\approx\mathcal{R}^{k+1}[\varPsi]+(a-a_{\infty})\mathcal{J}_{k+1}[\psi]=Q_{k+1}$.
Operating $\mathcal{R}$ on the induction step assumption yields
$\mathcal{R}^{k+1}[F_{a}]\approx\mathcal{R}[Q_{k}]\approx Q_{k+1}$, which
completes the proof of (39).

The leading eigenvalue of the eigenvalue problem
$\mathcal{L}_{g}[\varphi]=\delta\varphi$ is $\delta$ (i.e. Feigenbaum’s first
constant) with associated eigenfunction $\varphi$. For sufficiently large $n$,
the relation

$\displaystyle\mathcal{J}_{n}[\psi]\approx\delta^{n}B\varphi,$ (41)

with $B$ some constant, is expected to hold, because $\mathcal{J}_{n}$ results
from the repeated application, many times, of $\mathcal{L}_{g}$, which
amplifies the most unstable mode. The eigenvalue problem can be recast in
terms of the universal functions as $\mathcal{L}_{G}[\varPhi]=\delta\varPhi$,
where

$\displaystyle\varPhi(\xi)=\mu B\varphi(\gamma\xi)/\gamma$ (42)

is normalised to $\varPhi(0)=1$ choosing the scaling factor appropriately as
$\mu=\gamma/B\varphi(0)$.

Substituting (35) and (41) in (39) yields

$\displaystyle\mathcal{R}^{n}[F_{a}](y)\approx
g(y)+(a-a_{\infty})\delta^{n}B\varphi(y).$ (43)

Further replacing $g$ and $\phi$ using (36) and (42), respectively, we obtain

$\displaystyle\mathcal{R}^{n}[F_{a}](y)=\alpha^{n}F^{2^{n}}_{a}\left(\frac{y}{\alpha^{n}}\right)\approx\gamma\left\\{G\left(\frac{y}{\gamma}\right)+\delta^{n}(a-a_{\infty})\mu^{-1}\varPhi\left(\frac{y}{\gamma}\right)\right\\},$
(44)

where we have made explicit the effect of the renormalisation operator upon
repeated application. In order to get back to the original map, we can use in
(44) the obvious proposition that repeated application of the original and
shifted maps is related by $F_{a}^{n}(y)=f_{a}^{n}(y+X)-X$ for any
$n\in\mathbb{N}$. Combined with the definition of the universal function (28),
the approximation

$\displaystyle f^{2^{n}}_{a}(x)\approx
X+\Gamma_{n}U_{M}\left(\frac{x-X}{\Gamma_{n}}\right),$ (45)

follows directly, where the rescaled parameter $M$ and magnification rate of
the variable $\Gamma_{n}$ have been defined as

$\displaystyle
M=\delta^{n}\mu^{-1}(a-a_{\infty})\;\;\text{and}\;\;\Gamma_{n}=\gamma/\alpha^{n}.$
(46)

Let us now consider what conditions are necessary for the approximation (45)
to be valid, besides the obvious condition that $2^{n}$ must be large. The
reliability of (45) depends on the degree of trust one can place on the
approximations (35) and (41). It is important to note that, in renormalisation
theory, comparisons of functions are conducted on magnified scales. That is,
while (35) may be a good approximation on the scale where
$y/\alpha^{n}=(x-X)/\alpha^{n}$ is not too large, it may not hold true on the
original scale where $y\sim O(1)$ (mathematically, this can be easily
understood by considering the case where $\varPsi$ is topologically conjugate
to $g$). In (45), it suffices to consider $U_{M}(\xi)$ with the argument in
the range $\xi\in[-1,1]$, which encompasses the entire period doubling cascade
of (29). Hence the restriction $|(x-X)/\Gamma_{n}|<1$, which makes (35) only
locally valid. The approximation (41) is, again, local. Moreover, in order to
use the approximation (40) at level $k=n$, $|\delta^{n}(a-a_{\infty})|$ must
necessarily be small, so we require $|\delta^{n}M|\ll 1$.

The approximation (45) must of course also hold for $U_{M}(\xi)$ itself, as it
is unimodal in the range $\xi\in[-1,1]$. We can therefore set $n=m$,
$x=\tilde{x}$, $a=\tilde{M}$, $X=0$, $a_{\infty}=0$, $\mu=\tilde{\mu}$,
$\gamma=\tilde{\gamma}$ and $f_{a}=U_{\tilde{M}}$ in (45)-(46) to obtain

$\displaystyle
U^{2^{m}}_{\tilde{M}}(\tilde{x})\approx\tilde{\Gamma}_{m}U_{M}\left(\frac{\tilde{x}}{\tilde{\Gamma}_{m}}\right),\qquad
M=\delta^{m}\tilde{M}/\tilde{\mu},\qquad\tilde{\Gamma}_{m}=\tilde{\gamma}/\alpha^{m},$
(47)

which should work as long as $2^{m}$ is large,
$|\tilde{x}/\tilde{\Gamma}_{m}|<1$, and $|\delta^{m}\tilde{M}|\ll 1$.

Finally, particularising (45) for $n=\ell+m$, we get

$\displaystyle\frac{f^{2^{\ell+m}}_{a}(x)-X}{\Gamma_{\ell+m}}\approx
U_{M}\left(\frac{x-X}{\Gamma_{\ell+m}}\right)\approx\frac{\alpha^{m}}{\tilde{\gamma}}U_{\tilde{M}}^{2^{m}}\left(\frac{x-X}{\gamma\tilde{\Gamma}_{\ell}}\right),$
(48)

where the rightmost approximation is found by setting
$\tilde{x}=(\tilde{\gamma}/\gamma)(x-X)\alpha^{-\ell}$ in (47). The parameters
$a$ and $\tilde{M}$ are linked by
$\tilde{M}=\tilde{\mu}\delta^{-m}M=(\tilde{\mu}/\mu)\delta^{\ell}(a-a_{\infty})$.
We can simplify the notation by renaming ($\gamma/\tilde{\gamma}$) and
($\mu/\tilde{\mu}$) as $\gamma$ and $\mu$, respectively. In summary, the
approximation

$\displaystyle\frac{f^{2^{\ell+m}}_{a}(x)-X}{\Gamma_{\ell}}\approx
U_{M}^{2^{m}}\left(\frac{x-X}{\Gamma_{\ell}}\right),$ (49)

is obtained if we redefine

$\displaystyle
M=\delta^{\ell}\mu^{-1}(a-a_{\infty}),\qquad\Gamma_{\ell}=\gamma/\alpha^{\ell}.$
(50)

This approximation is valid for $2^{m}$ large, irrespective of the value of
$\ell=0,1,2\dots$, as long as

$\displaystyle\left|\frac{x-X}{\Gamma_{\ell+m}}\right|<1,\qquad|\delta^{\ell+m}\mu^{-1}(a-a_{\infty})|\ll
1.$ (51)

The approximation (26) and conditions (27) quoted in the main text correspond
to the case $\ell=0$.

The self-similarity of the period-doubling cascade follows from the validity
of the approximation (49). Suppose we detect a period-$2^{m}$ stable orbit in
(29) for $M$ negative but close to zero. Then the relation (49) asserts the
presence of period-$2^{\ell+m}$ orbits in (25) at
$a=a_{\ell}=\delta^{-\ell}\mu M+a_{\infty}$ for every $\ell\in\mathbb{N}$.
Moreover, the orbits must share the same multiplier, because an argument
analogous to (30) can be applied.

The same self-similarity theory applies also to values of the parameter beyond
the accumulation point. For instance, suppose we fix $M$ at a small positive
number so that a period $2^{m}p$ window is observed in (29), with $p$ an odd
number. Then we can expect a period $2^{\ell+m}p$ window in (25) at
$a_{\ell}=\delta^{-\ell}\mu M+a_{\infty}$ for every $\ell\in\mathbb{N}$. Even
if the dynamics are predominantly chaotic, we can still expect similarity of
the dynamics in the statistical sense. In particular, if $\lambda(M_{*})$ is
the Lyapunov exponent of (29) at $M=M_{*}$, the Lyapunov exponent of (25) at
$a_{\ell}=\delta^{-\ell}\mu M_{*}+a_{\infty}$ should approximately be
$2^{-\ell}\lambda(M_{*})$. This is consistent with the observation by Huberman
& Rudnick (1980). Note that this implies that for all unimodal maps, the
distribution of Lyapunov exponents is identical immediately past the
accumulation point.

## References

  * Andereck et al. (1986) Andereck, C. D., Liu, S. S. & Swinney, H. L. 1986 Flow regimes in a circular Couette system with independently rotating cylinders. J. Fluid Mech. 164, 155–183.
  * Ayats et al. (2020) Ayats, R., Deguchi, K., Mellibovsky, F. & Meseguer, A. 2020 Fully nonlinear mode competition in magnetised Taylor-Couette flow. J. Fluid Mech. 897, A14.
  * Barkley & Tuckerman (2007) Barkley, D. & Tuckerman, L. S. 2007 Mean flow of turbulent-laminar patterns in plane Couette flow. J. Fluid Mech. 576, 109–137.
  * Belmonte et al. (1988) Belmonte, A. L., Vinson, M. J., Glazier, J. A., Gunaratne, G. H. & Kenny, B. G. 1988 Trajectory scaling functions at the onset of chaos: Experimental results. Phys. Rev. Lett. 61, 539–542.
  * Briggs (1991) Briggs, K. 1991 A precise calculation of the feigenbaum constants. Mathematics of computation 57 (195), 435–439.
  * Budanur et al. (2015) Budanur, N. B., Cvitanović, P., Davidchack, R. L. & Siminos, E. 2015 Reduction of SO(2) symmetry for spatially extended dynamical systems. Physical Review Letters 114, 084102\.
  * Buzug et al. (1993) Buzug, Th., von Stamm, J. & Pfister, G. 1993 Characterization of period-doubling scenarios in Taylor-Couette flow. Phys. Rev. E 47 (2), 1054.
  * Cheng et al. (2020) Cheng, L., Ju, X., Tong, F. & An, H. 2020 Transition to chaos through period doublings of a forced oscillating cylinder in steady current. J. Fluid Mech. 887, A5.
  * Coles (1965) Coles, D. 1965 Transition in circular Couette flow. J. Fluid Mech. 21, 385–425.
  * Collet (2019) Collet, P. 2019 A short historical account of period doublings in the pre-renormalization era. Comptes Rendus Mécanique 347, 287–293.
  * Collet & Eckmann (1980) Collet, P. & Eckmann, J-P 1980 Iterated maps on the interval as dynamical systems. Birkhäuser.
  * Coullet & Tresser (1978) Coullet, P. & Tresser, C. 1978 Iterations d’endomorphismes et groupe de renormalisation. J. Physique IV (Proceedings) Le Journal de Physique IV 39, C5–25–C5–28.
  * Cvitanović (1989) Cvitanović, P. 1989 Universality in Chaos, 2nd edn. Taylor & Francis.
  * Cvitanović et al. (2016) Cvitanović, P., Rrutuso, R., Mainieri, R., Tanner, G. & Vattay, G. 2016 Chaos: Classical and Quantum. Copenhagen: Niels Bohr Inst.
  * Deguchi & Altmeyer (2013) Deguchi, K. & Altmeyer, S. 2013 Fully nonlinear mode competitions of nearly bicritical spiral or Taylor vortices in Taylor-Couette flow. Phys. Rev. E 87, 043017.
  * Deguchi et al. (2014) Deguchi, K., Meseguer, A. & Mellibovsky, F. 2014 Subcritical equilibria in Taylor-Couette flow. Phys. Rev. Lett. 112, 184502.
  * Ding et al. (2016) Ding, X., Chaté, H., Cvitanović, P., Siminos, E. & Takeuchi, K. A. 2016 Estimating the dimension of an inertial manifold from unstable periodic orbits. Phys. Rev. Lett. 117, 024101.
  * Dong (2009) Dong, S. 2009 Evidence for internal structures of spiral turbulence. Phys. Rev. E 80, 067301.
  * Eckmann & Wittwer (1987) Eckmann, J. P. & Wittwer, P. 1987 A complete proof of the Feigenbaum conjectures. J. Stat. Phys. 46, 455–475.
  * Epstein (1986) Epstein, H. 1986 New proofs of the existence of the Feigenbaum functions. Commum. Math. Phys. 106, 395–426.
  * Feigenbaum (1978) Feigenbaum, M. J. 1978 Quantitative universality for a class of nonlinear transformations. Journal of Statistical Physics 19, 25–52.
  * Feigenbaum (1979) Feigenbaum, M. J. 1979 The onset spectrum of turbulence. Physics Letters A 74 (6), 375–378.
  * Feigenbaum (1980) Feigenbaum, M. J. 1980 The transition to aperiodic behavior in turbulent systems. Communications in mathematical physics 77, 65–86.
  * Feigenbaum (1982) Feigenbaum, M. J. 1982 Some formalism and predictions of the period-doubling onset of chaos. In North-Holland Mathematics Studies, , vol. 61, pp. 379–394. Elsevier.
  * Feigenbaum (1983) Feigenbaum, M. J. 1983 Universal behavior in nonlinear systems. Physica D: Nonlinear Phenomena 7 (1), 16–39.
  * Feldmann et al. (2023) Feldmann, D., Borrero-Echeverry, D., Burin, M. J., Avila, K. & Avila, M. 2023 Routes to turbulence in Taylor-Couette flow. Phil. Trans. R. Soc. A 381, 20220114\.
  * Gao et al. (2015) Gao, Z., Podvin, B., Sergent, A. & Xin, S. 2015 Chaotic dynamics of a convection roll in a highly confined, vertical, differentially heated fluid layer. Phys. Rev. E 91, 013006\.
  * Garcia et al. (2021) Garcia, F., Seilmayer, M., Giesecke, A. & Stefani, F. 2021 Long term time dependent frequency analysis of chaotic waves in the weakly magnetized spherical Couette system. Physica D: Nonlinear Phenomena 418, 132836.
  * Giglio et al. (1981) Giglio, M., Musazzi, S. & Perini, U. 1981 Transition to chaotic behavior via a reproducible sequence of period-doubling bifurcations. Phys. Rev. Lett. 47, 243–246.
  * Glendinning (1994) Glendinning, P. 1994 Stability, Instability and Chaos: An Introduction to the Theory of Nonlinear Differential Equations. Cambridge University Press.
  * Grossmann & Thomae (1977) Grossmann, S. & Thomae, S. 1977 Invariant distributions and stationary correlation functions of one-dimensional discrete processes. Zeitschrift für Naturforschung A 32.
  * Hall & Sherwin (2010) Hall, P. & Sherwin, S. 2010 Streamwise vortices in shear flows: harbingers of transition and the skeleton of coherent structures. J. Fluid Mech. 661, 178–205.
  * Haller et al. (2023) Haller, G., Kaszás, B., Liu, A. & Axås, J. 2023 Nonlinear model reduction to fractional and mixed-mode spectral submanifolds. Chaos: An Interdisciplinary Journal of Nonlinear Science 33 (6), 063138.
  * Huberman & Rudnick (1980) Huberman, B. A. & Rudnick, J. 1980 Scaling behavior of chaotic flows. Phys. Rev. Lett. 45, 154–156.
  * Kawahara et al. (2012) Kawahara, G., Uhlmann, M. & van Veen, L. 2012 The significance of simple invariant solutions in turbulent flows. Ann. Rev. Fluid Mech. 44 (1), 203–225.
  * Khanin et al. (2021) Khanin, K., Lyubich, M., Siggia, E. D. & Sinai, Y. 2021 Mitchell feigenbaum. Notices of the American Mathematical Society 68 (5), 757–767.
  * Kreilos & Eckhardt (2012) Kreilos, T. & Eckhardt, B. 2012 Periodic orbits near onset of chaos in plane Couette flow. CHAOS 22, 047505\.
  * Lanford III (1982) Lanford III, O. E. 1982 A computer-assisted proof of the feigenbaum conjectures. Bull. American Math. Soc. 6, 427–435.
  * Li & Yorke (1975) Li, T. Y. & Yorke, J. A. 1975 Period three implies chaos. Am. Math. Month. 82, 985–992.
  * Libchaber (1983) Libchaber, A. 1983 Experimental aspects of the period doubling scenario. In Dynamical System and Chaos (ed. Luis Garrido), pp. 157–164. Berlin, Heidelberg: Springer Berlin Heidelberg.
  * Libchaber (2023) Libchaber, A. 2023 A journey through nonlinear dynamics: The case of temperature gradients. Ann. Rev. Cond. Matt. Phys. 14, 1–19.
  * Libchaber et al. (1982) Libchaber, A., Laroche, C. & Fauve, S. 1982 Period doubling cascade in mercury, a quantitative measurement. J. Physique Lett. 43 (7), 211–216.
  * Lizée & Alexander (1997) Lizée, A. & Alexander, J. Iwan D. 1997 Chaotic thermovibrational flow in a laterally heated cavity. Phys. Rev. E 56, 4152–4156.
  * Lopez & Marques (2005) Lopez, Juan M & Marques, Francisco 2005 Finite aspect ratio Taylor–Couette flow: Shil’nikov dynamics of 2-tori. Physica D: Nonlinear Phenomena 211 (1-2), 168–191.
  * Lorenz (1980) Lorenz, E. N. 1980 Noisy periodicity and reverse bifurcation. Annals of The Lyceum of Natural History of New York 357, 282–291.
  * Lustro et al. (2019) Lustro, J.R.T., Kawahara, G., van Veen, L., Shimizu, M. & Kokubu, H. 2019 The onset of transient turbulence in minimal plane Couette flow. J. Fluid Mech. 862, R2.
  * Meseguer et al. (2007) Meseguer, A., Avila, M., Mellibovsky, F. & Marques, F. 2007 Solenoidal spectral formulations for the computation of secondary flows in cylindrical and annular geometries. Eur. Phys. J. Special Topics 146, 249–259.
  * Meseguer et al. (2009a) Meseguer, A., Mellibovsky, F., Avila, M. & Marques, F. 2009a Families of subcritical spirals in highly counter-rotating Taylor-Couette flow. Phys. Rev. E 79, 036309.
  * Meseguer et al. (2009b) Meseguer, A., Mellibovsky, F., Avila, M. & Marques, F. 2009b Instability mechanisms and transition scenarios of spiral turbulence in Taylor-Couette flow. Phys. Rev. E 80, 046315.
  * Moore et al. (1983) Moore, D. R., Toomre, J., Knobloch, E. & Weiss, N. O. 1983 Period doubling and chaos in partial differential equations for thermosolutal convection. Nature 303, 663–667.
  * Roberts & Mackley (1996) Roberts, E. P. L. & Mackley, M. R. 1996 The development of asymmetry and period doubling for oscillatory flow in baffled channels. J. Fluid Mech. 328, 19.
  * Schuster (1988) Schuster, H. G 1988 Deterministic Chaos, An Introduction, 2nd Ed.. Weinheim: VCH Publishers.
  * Sharkovskii (1964) Sharkovskii, A. N. 1964 Co-existence of cycles of a continuous mapping of the line into itself. Ukrainian Math. J. 16, 61–71.
  * Smyrlis & Papageorgiou (1991) Smyrlis, Y. S. & Papageorgiou, D. T. 1991 Predicting chaos for infinite dimensional dynamical systems: the kuramoto-sivashinsky equation, a case study. Proceedings of the National Academy of Sciences 88 (24), 11129–11132.
  * Sthephenson & Wang (1991) Sthephenson, J. & Wang, Y. 1991 Relationships between eigenfunctions associated with solutions Feigenbaum’s equation. Appl. Math. Lett. 4 (3), 53–56.
  * Strogatz (2024) Strogatz, S. H. 2024 Nonlinear Dynamics and Chaos - With Applications to Physics, Biology, Chemistry, and Engineering, 3rd Ed.. CRC Press.
  * Swinney & Gollub (1985) Swinney, H. L. & Gollub, J. P. 1985 Hydrodynamic Instabilities and the Transition to Turbulence. Topics in Applied Physics 45\. Springer-Verlag, Berlin, Heidelberg.
  * Temam (1989) Temam, R. 1989 Do inertial manifolds apply to turbulence? Physica D: Nonlinear Phenomena 37, 146–152.
  * Thurlby (2021) Thurlby, J. A. 2021 Rigorous calculations of renormalisation fixed points and attractors. PhD Thesis, University of Portsmouth, UK, Portsmouth.
  * Tirnakli et al. (2009) Tirnakli, U., Tsallis, C. & Beck, C. 2009 Closer look at time averages of the logistic map at the edge of chaos. Phys. Rev. E 79 (5), 056209.
  * Van Veen (2005) Van Veen, L. 2005 The quasi-periodic doubling cascade in the transition to weak turbulence. Physica D: Nonlinear Phenomena 210 (3-4), 249–261.
  * Wang et al. (2022) Wang, B., Ayats, R., Deguchi, K., Mellibovsky, F. & Meseguer, A. 2022 Self-sustainment of coherent structures in counter-rotating Taylor-Couette flow. J. Fluid Mech. 951, A21.
  * Wang et al. (2023) Wang, B., Mellibovsky, F., Ayats, R., Deguchi, K. & Meseguer, A. 2023 Mean structure of the supercritical turbulent spiral in Taylor–Couette flow. Phil. Trans. R. Soc. A 381 (2246), 20220112.
  * Wang et al. (2007) Wang, J., Gibson, J. & Waleffe, F. 2007 Lower branch coherent states in shear flows: transition and control. Phys. Rev. Lett. 98 (20), 204501.
  * Yalim et al. (2019) Yalim, J., Welfert, Bruno D. & Lopez, J.M. 2019 Parametrically forced stably stratified cavity flow: complicated nonlinear dynamics near the onset of instability. J. Fluid Mech. 871, 1067–1096.