2.1. Flow configuration – governing equations

The flow past an axisymmetric rotating body is controlled by two parameters: the Reynolds number ($\mbox { {Re}}$) and the rotation rate (${{\varOmega }}$) which is defined as the ratio of the tangential velocity ${{\varOmega }^{\ast }}D^{*}/2$ on the sphere surface to the inflow velocity $W^{\ast }_\infty$. The fluid motion inside the domain is governed by the incompressible Navier–Stokes equations written in cylindrical coordinates $(r,\theta,z)$,

(2.1a)$$\begin{gather} \frac{\partial \boldsymbol{U}}{\partial t} + \boldsymbol{U}\boldsymbol{\cdot} \boldsymbol{\nabla}\boldsymbol{U} ={-}\boldsymbol{\nabla} P + \boldsymbol{\nabla} \boldsymbol{\cdot} \tau(\boldsymbol U ), \quad \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol U = 0, \end{gather}$$

(2.1b)$$\begin{gather}\text{with } \tau(\boldsymbol{U}) = \frac{1}{\mbox{{Re}}} (\boldsymbol{\nabla} \boldsymbol U + \boldsymbol{\nabla} \boldsymbol U^{\rm T}), \quad \mbox{{Re}} = \frac{W^{{\ast}}_\infty D^{{\ast}}}{\nu^{{\ast}}}, \quad \varOmega = \frac{{{\varOmega}^{{\ast}}}D^{{\ast}}}{2 W^{{\ast}}_\infty}, \end{gather}$$

(2.1c)$$\begin{gather}\boldsymbol{x} = \boldsymbol{x}^{{\ast}} \frac{1}{D^{{\ast}}}, \quad t = t^{{\ast}} \frac{W^{{\ast}}_\infty}{D^{{\ast}}}, \quad \boldsymbol{U} = \boldsymbol{U}^{{\ast}} \frac{1}{W^{{\ast}}_\infty} \quad P = P^{{\ast}} \left(\frac{1}{W^{{\ast}}_\infty}\right)^{2}. \end{gather}$$

Dimensional quantities are identified with the upperscript symbol $^{\ast }$. Reference scales are specified in (2.1c). The dimensionless velocity vector $\boldsymbol {U} = (U,V,W)$ is composed of the radial, azimuthal and axial components, $P$ is the dimensionless reduced pressure and the viscous stress tensor, $\tau (\boldsymbol U )$. For representation purposes, it is sometimes necessary to use the Cartesian coordinates ($x,y,z$), here $z$ denotes the streamwise direction, $y$ the vertical crosswise direction and $x$ the direction that forms a direct trihedral with $z$ and $y$. The incompressible Navier–Stokes equations (2.1) are complemented with the following boundary conditions:

(2.2)\begin{equation} \boldsymbol{U} = (0, \varOmega, 0) \text{ on } \varSigma_{b} \quad \boldsymbol{U} = (0, 0, 1) \text{ on } \varSigma_{i}. \end{equation}

No-slip boundary condition is set on the rotating sphere and a uniform boundary condition is set in the inlet, as shown in figure 1.

Figure 1. Sketch of the problem and geometric configuration.

In the sequel, Navier–Stokes equations (2.1) and the associated boundary conditions will be written symbolically under the form

(2.3)\begin{equation} \boldsymbol{B} \frac{\partial \boldsymbol{Q}}{\partial t} = \boldsymbol{F}(\boldsymbol{Q}, \boldsymbol{\eta}) \equiv \boldsymbol{L} \boldsymbol{Q} + \boldsymbol{N}(\boldsymbol{Q},\boldsymbol{Q}) + \boldsymbol{G} (\boldsymbol{Q},\boldsymbol{\eta}), \end{equation}

where $\boldsymbol {B}$ is the projection matrix onto the velocity field with the flow state vector $\boldsymbol {Q} = [ \boldsymbol {U},{P} ]^{\rm T}$, and the parameter vector $\boldsymbol {\eta } = [\mbox { {Re}}^{-1}, {{\varOmega }} ]^{\rm T}$. Such a form of the governing equations takes into account a linear dependency on the state variable $\boldsymbol {Q}$ through $\boldsymbol {L}$ and a quadratic dependency on parameters and the state variable through operators $\boldsymbol {N}(\cdot, \cdot )$ and $\boldsymbol {G}(\cdot, \cdot )$, which are detailed in Appendix A.

2.2. Nomenclature

Let us introduce some general concepts that will be employed throughout the study. Steady states, i.e. $\boldsymbol {Q}$ such that $\boldsymbol {F}(\boldsymbol {Q}, \boldsymbol {\eta }) = 0$, periodic orbits, i.e. $\boldsymbol {Q}(t) = \boldsymbol {Q}(t+T)$ for every $t \geq 0$, are the simplest invariants of (2.3). In general, an invariant set $V$ of the phase space of (2.3) is a set that is preserved under dynamics, i.e. for every initial solution $\boldsymbol {Q}(t_0) \in V$, we have $\boldsymbol {Q}(t) \in V$ for every $t \geq 0$. A $T^{n}$-quasiperiodic state, $n > 1,$ $n \in \mathbb {N}^{*}$, is an invariant of the system (2.3) that can be decomposed as a finite sum of $n$ incommensurate frequencies $\omega _n$, i.e.

(2.4)\begin{equation} \boldsymbol{Q} = \boldsymbol{Q}_0 + \sum_{\ell=1}^{n} \left( \hat{\boldsymbol{Q}}_{\ell} {\rm e}^{\mathrm{i} \omega_{\ell} t} + \text{c.c.} \right). \end{equation}

Incommensurate frequencies are those that are linearly independent, i.e. for $k_{\ell } \in \mathbb {Z}$, we have $\sum _{\ell =1}^{n} k_{\ell } \omega _{\ell } = 0$ if and only if every $k_{\ell } = 0$. Here, we determine the incommensurate frequencies as those corresponding to the fundamental modes (least stable eigenmodes) identified by LSA.

A second important property is the attractiveness of an invariant set. We denote as basin of attraction the set of initial conditions leading to long-time behaviour that approaches the attractor. The celebrated manuscript of Newhouse, Ruelle & Takens (Reference Newhouse, Ruelle and Takens1978) states that $T^{n}$-quasiperiodic states, with $n \geq 3$, are unusual attractors, in the sense that every $T^{n}$-quasiperiodic state can be perturbed by an arbitrarily small amount to a new vector field with a chaotic attractor. In other words, for any $T^{n}$-quasiperiodic state of (2.3), one may observe a chaotic Axiom A attractor by experimental or numerical means. Here, Axiom A attractor denotes a class of dynamical systems where the non-wandering set is hyperbolic and the attractor has a dense set of periodic orbits, more details about hyperbolicity may be found, for instance, in the recent article by Ni (Reference Ni2019).

2.3. Direct numerical simulation details

The flow governed by (2.1) is solved by means of DNS, following a time-stepping approach using the finite-volume library OpenFOAM$\circledR$. The domain shown in figure 1 consists of an upstream hemisphere of radius $r_{\infty }=15D$ and a downstream tube extending $z_{2}=50D$ downstream of the body.

Regarding boundary conditions at the outlet, $\varSigma _{o}$, we impose an outflow condition that implements a Neumann condition for the velocity, $\boldsymbol {n}\boldsymbol {\cdot } \boldsymbol {\nabla }\boldsymbol {U} = 0$, where $\boldsymbol {n}$ is the outward normal, and a Dirichlet condition for the pressure, $P = 0$. The latter may be considered equivalent to setting a stress-free condition at the outlet for small values of the viscosity (as highlighted by Tomboulides & Orszag Reference Tomboulides and Orszag2000). Finally, at the outer radial boundary, $\varSigma _w$, we set a slip boundary condition, $\boldsymbol {n} \boldsymbol {\cdot } \boldsymbol {U}=0$. Note that such domain size and boundary conditions have been selected according to previous numerical works on rotating axisymmetric bodies (see, e.g. Jiménez-González et al. Reference Jiménez-González, Sanmiguel-Rojas, Sevilla and Martínez-Bazán2013; Lorite-Díez & Jiménez-González Reference Lorite-Díez and Jiménez-González2020). Additionally, second-order schemes have been employed for spatial and time integration. Nevertheless, for the sake of conciseness, the reader is referred to Appendix A in Lorite-Díez & Jiménez-González (Reference Lorite-Díez and Jiménez-González2020) for detailed information about the employed numerical schemes, convergence and validation studies. In the present simulations $\sim$2.6 millions of elements mesh, denoted $\#2$ in table 1 (Appendix A) therein, is used.

The three-dimensional time-stepping simulations were computed in parallel. In particular, the DNS are carried out, once converged, for $T\sim 500$ convective units for periodic regimes, and until $T\sim 1000$ convective units for quasiperiodic and most complex regimes. The employed time step is $\Delta t = 0.003$ for all simulations. In terms of computational cost, running on 16 Intel Xeon E5-2665 processors, a simulation lasting $T=1000$ convective time units corresponds to approximately 10 days.

3.1. Methodology

As a first step of the reduction procedure, we identify the base flow solution, which is defined as the steady solution $\boldsymbol {Q}_b$ of the (axisymmetric) Navier–Stokes equations, namely the solution of $\boldsymbol {F} ( \boldsymbol {Q}_b ) = \boldsymbol {0}$. We then characterize the dynamics of small-amplitude perturbations around this base flow by expanding them over the basis of linear eigenmodes

(3.1)\begin{equation} \boldsymbol{Q} = \boldsymbol{Q}_b + \varepsilon \sum_{\ell} \boldsymbol{q}_{(\varepsilon)}(t,\tau) = \boldsymbol{Q}_b + \varepsilon \sum_{\ell} \left( z_{\ell}(\tau) \hat{\boldsymbol{q}}_{(z_{\ell})}(r,z) \text{e}^{\mathrm{i}( m_{\ell} \theta + \omega_{\ell} t)} + \text{c.c.} \right), \quad \varepsilon \ll 1. \end{equation}

The eigenpairs $[\mathrm {i} \omega _{\ell }, \hat {\boldsymbol {q}}_{(z_{\ell })} ]$ are then determined as the solutions of the eigenvalue problem

(3.2)\begin{equation} \displaystyle \boldsymbol{J}_{(\omega_{\ell}, m_{\ell})} \hat{\boldsymbol{q}}_{(z_{\ell})} = \left({\rm i}\omega_{\ell} \boldsymbol{B} -\left. \frac{\partial \boldsymbol{F}}{\partial \boldsymbol{q} }\right|_{\boldsymbol{q} = \boldsymbol{Q}_b, \Delta \boldsymbol{\eta} = \boldsymbol{0}} \right) \hat{\boldsymbol{q}}_{(z_{\ell})}, \end{equation}

where $({\partial \boldsymbol {F} }/{\partial \boldsymbol {q} }|_{\boldsymbol {q} = \boldsymbol {Q}_b, \Delta \boldsymbol {\eta } = \boldsymbol {0}} ) \hat {\boldsymbol {q}}_{(z_{\ell })} \!=\! \boldsymbol {L}_{m_{\ell }} \hat {\boldsymbol {q}}_{(z_{\ell })} \!+\! \boldsymbol {N}_{m_{\ell }}(\boldsymbol {Q}_b, \hat {\boldsymbol {q}}_{(z_{\ell })} ) \!+\! \boldsymbol {N}_{m_{\ell }}(\hat {\boldsymbol {q}}_{(z_{\ell })}, \boldsymbol {Q}_b )$ $+ \boldsymbol {G}(\boldsymbol {Q}_b, \boldsymbol {\eta }_c)$, with $\boldsymbol {\eta }_c = [Re_c^{-1}, \varOmega _c]^{\rm T}$. The subscript $m_{\ell }$ indicates the azimuthal wavenumber used for the evaluation of the linearized Navier–Stokes operator $\boldsymbol {J}_{(\omega _{\ell }, m_{\ell })}$. Please note that here, the term $\Delta \boldsymbol {\eta } = [Re_c^{-1} - Re^{-1}, \varOmega _c - \varOmega ]^{\rm T}$ denotes the departure from the critical condition attained at $[Re_c^{-1}, \varOmega _c]^{\rm T}$. In the following, we consider that eigenmodes $\hat {\boldsymbol {q}}_{(z_{\ell })}(r,z)$ have been normalised in such a way that $\langle \boldsymbol {\hat q}_{(z_{\ell })}, \boldsymbol {\hat q}_{(z_{\ell })} \rangle _{\boldsymbol {B}} = \langle \boldsymbol {\hat u}_{(z_{\ell })}, \boldsymbol {\hat u}_{(z_{\ell })} \rangle = \int _{\varOmega } \boldsymbol {u}(\boldsymbol {x})^{\rm T} \boldsymbol {u}(\boldsymbol {x}) \,\text {d}\kern0.06em \boldsymbol {x} = 1$.

3.2. Neutral curves of stability

In the presence of supercritical self-sustained instabilities, rotating waves are predominant. These patterns prevail in axisymmetric flows, where the reflection symmetry regarding the azimuthal angle is broken. Here, the reflection symmetry is broken because of the rotation of the sphere, which induces a preferential direction of rotation. Consequently, bifurcations that lead to standing waves or to a symmetry breaking steady state do not occur generically. The existence of standing waves or a steady-state mode requires the matching between the phase speed of the helical pattern and the rotation of the body, which is another condition to be met. The global stability analysis of the flow past the sphere confirms that only rotating waves are linearly unstable for the range of Reynolds numbers $\mbox { {Re}} < 300$ and $\varOmega < 4$. The parametric linear stability study of the flow past the rotating sphere shows the existence of three neutral curves, which are associated to the three least stable modes identified by global stability analysis. These correspond to rotating waves, named $RW_1$, $RW_2$ and $RW_3$, which are depicted in figure 2. Linear stability results (figure 3a) reveal that the axisymmetric steady state, referred in the following as a trivial state, is stable in the white shaded region and unstable in the grey shaded region. The neutral curve of stability displays two regions in the parameter space $(\mbox { {Re}}, {{\varOmega }})$ for which the first primary bifurcations are rotating waves of low frequency where the wake past the sphere displays a single helix ($RW_1$), depicted in figure 2(a). In the second region, the flow pattern of the wake displays a double helix ($RW_3$) with a high frequency, depicted in figure 2(c). The onset of instability of the third branch ($RW_2$) displaying a flow pattern of the wake with a single helix with a medium frequency, depicted in figure 2(b), turns out to be linearly unstable for $\varOmega \leq 4$. Each pair of neutral curves intersects once, leading to three codimension-two points ($A$, $B$, $C$), identified in table 1. Another aspect of importance is the evolution of frequencies of the instability. Frequencies at critical parameters are reported in figure 3(b) as a function of ${{\varOmega }}$. The frequency evolution is divided into two regions, a first of rapid evolution for low rotation rates ${{\varOmega }} < 1$ and a second where the frequency of the three modes hardly depends on the rotation rate.

Figure 2. Cross-section view at $z=3.5$ of the three unstable modes. The streamwise component of the vorticity vector $\varpi _z$ is visualized by colours. Results are shown for (a) $RW_1$ at point $(\mbox { {Re}}_A, \varOmega _A) = (77,2.24)$; (b) $RW_2$ at point $(\mbox { {Re}}_B, \varOmega _B)=(188,1.01)$; (c) $RW_3$ at point $(\mbox { {Re}}_A, \varOmega _A)$.

Figure 3. Linear stability properties of the rotating sphere configuration. (a) Neutral curve of stability: the onset of the primary instability is portrayed with a solid black line ($\textbf {---}$), whereas the continuation of the neutral curves is depicted with dashed black lines ($\textbf {- - -}$). (b) Frequency evolution with respect to ${{\varOmega }}$ of linear modes at the critical Reynolds number ($\mbox { {Re}}_c({{\varOmega }}$)).

Table 1. Location in the parameter space $(\mbox { {Re}},{{\varOmega }}$) and the pair of modes involved at the codimension-two points. It also lists the main properties of the primary bifurcation of the flow past the static sphere. The last two columns are related to non-normality effects and are defined in § 3.3.

The neutral curve of stability reveals that the static configuration ($\varOmega = 0$) exhibits the largest critical Reynolds number. Then, the critical value of the Reynolds number is hardly modified by weak rotating speeds, in the range ${{\varOmega }} < 0.3$. However, there is a clear threshold around ${{\varOmega }} \approx 0.4$ where the critical Reynolds number passes from around $\mbox { {Re}}_c \approx 200$ to $\mbox { {Re}}_c \approx 100$ in a narrow interval ${{\varOmega }} \in [0.4, 1.2]$. The critical Reynolds number remains approximately constant up to the point $A$, the point which divides the boundary of stability. Below the point $A$, that is, for ${{\varOmega }} < {{\varOmega }}_A$, the steady-state flow transits supercritically to a single helix rotating wave $RW_{1}$; above the point $A$, i.e. ${{\varOmega }} > {{\varOmega }}_A$, the steady-state flow transits supercritically to the double helix rotating wave, $RW_{3}$. Such a point corresponds to a double-Hopf bifurcation between modes 1 and 3, and its analysis is left to §§ 4 and 5. Other two double-Hopf bifurcation points exist, denoted $B$ and $C$, which characterize the interaction between modes 2 and 3, and 1 and 2, respectively. Yet, at points $B$ and $C$ the trivial state is already unstable, thus, instabilities associated with these points are not directly observed in experiments or numerical simulations. Instead, these organizing centres play a role in the pattern formation of secondary instabilities, which is left to §§ 4 and 5, where we interpret the subtle implications of these points in dynamics. In addition, authors have looked for the presence of a primary bifurcation that leads to the $RW_{2}$ state. For the studied configuration, there does not exist such a region in the range $0 < {{\varOmega }} < 6$.

3.3. Properties of the axisymmetric steady state

The analysis presented in this section studies the linear stability of the axisymmetric steady-state solution in the range $\mbox { {Re}} \leq 250$ and $\varOmega \leq 4$. Typical axisymmetric steady-state solutions (TS) at codimension-two points are portrayed in figure 4, which shows the neutrally stable trivial state state at $(Re_A,\varOmega _A)=(77,2.24)$ and the two other unstable trivial states at $(Re_B,\varOmega _B)=(188,1.01)$ and $(Re_C,\varOmega _C)=(73,3.95)$, respectively. The flow visualization illustrates the recirculation region behind the sphere, delimited by the separatrix, which divides the recirculation bubble and the unperturbed flow field. Such a line, depicted with a thick solid line in figure 4 connects the separation point on the sphere surface and the stagnation point on the $r=0$ axis. The development of the recirculation bubble can be measured using the maximum extent of the region

(3.3)\begin{equation} L_r = \max \left\{ z - \frac{D}{2} \text{ } |\text{ } W(r=0,z) \leq 0 \right\}, \end{equation}

where $D$ is the diameter of the sphere. Figure 5(a) displays the evolution of the length of the recirculation bubble by varying $\varOmega$ and $\mbox { {Re}}$. The length of the bubble increases monotonically with the angular velocity $\varOmega$ of the sphere as well as the largest negative values of the streamwise velocity behind the sphere, from around $40\,\%$ for $\varOmega = 0$ to around $60\,\%$ for the largest values of $\varOmega$ explored. A similar trend was identified by Kim & Choi (Reference Kim and Choi2002) at $Re$=100; however, we should consider that the trends observed in figure 5(a) are only valid before bifurcation. After that, $L_r$ does not have to increase with $\varOmega$, as seen by Lorite-Díez & Jiménez-González (Reference Lorite-Díez and Jiménez-González2020) and Kim & Choi (Reference Kim and Choi2002). The results at the onset of stability of the steady state are synthesized in figure 5(b), with a domain of existence of a stable steady state (white shaded) and another of an unstable steady state (grey shaded). In § 3.4 we identify the core of the $RW_1$ and $RW_3$ instabilities, which are found within the recirculation region. In particular, a passive control that shortens the recirculation region is an efficient technique to stabilize the flow. Therefore, it is not surprising that the neutrally stable flow is characterized by a shorter recirculation region with respect to the unstable steady state.

Figure 4. Spatial distribution of the streamwise velocity (contours) at the steady state along with flow streamlines (solid lines) and recirculation region separatrix (thick solid lines). Results are shown for (a) $(\mbox { {Re}} = 212, \varOmega = 0)$ in the upper half and $(\mbox { {Re}}_A, \varOmega _A)$ in the lower half; (b) $(\mbox { {Re}}_B, \varOmega _B)$ in the upper half and $(\mbox { {Re}}_C, \varOmega _C)$ in the bottom half. Points $A$, $B$ and $C$ are defined in table 1.

Figure 5. Evolution of the recirculating length $L_r$ in the plane $(\mbox { {Re}}, {{\varOmega }})$. The unstable region is the grey-shaded area, delimited by a thick line. (a) The recirculating length of the stable solution is painted with solid lines, dashed lines are employed for the unstable steady state. (b) Length of the recirculating region at the onset of stability $(\mbox { {Re}}_c(\varOmega ))$.

Finally, we briefly discuss the influence of non-normality mechanisms, lift-up and convective non-normality as they are partly related to recirculation region length. The main results are included in table 1, where we can see a lower influence of non-normality effects through the obtained values for $\gamma$ and $\theta _{N}$, with respect to the static sphere configuration. The estimator $\theta _N$ measures the importance of non-normality, the lower $\theta _N$ the more important non-normal effects are. On the other hand, the estimator $\gamma$ characterizes the relative contribution between the lift-up and the convective non-normality mechanisms to the total non-normality effects. A $\gamma$ value close to $0$ indicates the dominance of the lift-up effect. The largest non-normal effects have been measured at point $B$ (lowest values of $\theta _N$), which corresponds to the point with the largest critical Reynolds number among the codimension-two points. The values of $\theta _N$ obtained at point $B$ are associated with a larger non-normality than the stationary mode (the case of $RW_1$ with $O(2)$ symmetry) and $RW_2$ at the threshold for $(\varOmega = 0, \mbox { {Re}} = 281)$, which was found by Meliga, Chomaz & Sipp (Reference Meliga, Chomaz and Sipp2009b) to be $1^{\circ }$. Thus, one may conclude that the rotation of the sphere increases the effect of non-normality, however, it induces an earlier transition with regard to the Reynolds number, which turns out to globally reduce the effect of non-normality. This previous statement can be also indirectly verified from the satisfactory comparison between normal form estimations and DNS results in § 5.1. Furthermore, the analysis of the direct global mode shows a dominant effect of the convective non-normality, which is responsible at most of around $90\,\%$ (mode $RW_1$) and $98\,\%$ (mode $RW_3$) at point $A$ and around $80\,\%$ for the remainder modes at points $B$ and $C$. In comparison, the stationary and oscillating modes of static configuration ($\varOmega = 0$) displayed $\gamma = 0.76$ and $\gamma = 0.94$. More details about the non-normality study such as the definition of $\theta _N$ and $\gamma$ can be found in Appendix B.

3.4. Identification of the physical mechanisms from a control perspective

In this section we analyse the physical mechanisms leading to the $RW_1$ and $RW_3$ states at the point $A$. However, we do not discuss the $RW_2$ state as it will be seen in § 5, this state is not expected to be observed. First, we consider what is the effect of a steady axisymmetric forcing term, which represents the presence of a small obstacle, wall suction/blowing (as the control applied in Niazmand & Renksizbulut Reference Niazmand and Renksizbulut2005), etc. In this case the governing equations of the resulting flow are the same as (2.3) with the addition of a forcing term $\boldsymbol {H}_0 \equiv \boldsymbol {\hat {H}}_0$,

(3.4)\begin{equation} \boldsymbol{B} \frac{\partial \boldsymbol{Q}}{\partial t} = \boldsymbol{F}(\boldsymbol{Q}, \boldsymbol{\eta}) \equiv \boldsymbol{L} \boldsymbol{Q} + \boldsymbol{N}(\boldsymbol{Q},\boldsymbol{Q}) + \boldsymbol{G} (\boldsymbol{Q},\boldsymbol{\eta}) + \boldsymbol{\hat{H}}_0. \end{equation}

This case has been treated in the past by Marquet, Sipp & Jacquin (Reference Marquet, Sipp and Jacquin2008) in the case of the flow past a circular cylinder and by Sipp (Reference Sipp2012) in the case of the open cavity flow. The introduction of the forcing induces a modification of the eigenvalue $\mathrm {i} \omega _{\ell } \mapsto \mathrm {i} \omega _{\ell } + \Delta \mathrm {i} \omega _{\ell }^{0}$, where $\Delta \mathrm {i} \omega _{\ell }^{0} = \langle {\boldsymbol {\nabla }_{\boldsymbol {H}_{0}} \mathrm {i} \omega _{\ell } } , {\boldsymbol {\hat {H}}} \rangle$. Therefore, the control that induces the largest deviation of the growth rate (respectively frequency) of the mode $\ell$ is in the direction of $\boldsymbol {\nabla }_{\boldsymbol {H}_{0}} \mathrm {i} \omega _{\ell }$, which is defined as

(3.5)\begin{equation} \boldsymbol{\nabla}_{\boldsymbol{H}_{0}} \mathrm{i} \omega_{\ell} = \boldsymbol{B}^{\rm T} \boldsymbol{J}_{(0,0)} \boldsymbol{B} \boldsymbol{\nabla}_{\boldsymbol{U}_b} \mathrm{i} \omega_{\ell}\quad \text{for } \ell = 1,2,3. \end{equation}

Here $\boldsymbol {\nabla }_{\boldsymbol {U}_b} \mathrm {i} \omega _{\ell }$ is the sensitivity of the eigenvalue of the mode $\ell$ ($\ell = 1,2,3$) with respect to variations in the axisymmetric steady state, cf. (Marquet et al. Reference Marquet, Sipp and Jacquin2008). The sensitivity of the $\ell ^{th}$ eigenvalue $\boldsymbol {\nabla }_{\boldsymbol {H}_{0}} \lambda _{\ell }$ to the introduction of a steady axisymmetric forcing is represented in figure 6 for the two modes present in the codimension point $A$. The low-frequency mode ($RW_1$) is most sensitive to a steady axisymmetric forcing at the leftmost end of the recirculation region (see figure 6a,b). This forcing corresponds to one that accelerates the streamwise motion at the end of the recirculation region, thus reducing the counterclockwise motion of the recirculation zone, which would induce an effective decrease of the growth rate (respectively frequency). This is in accordance with the fact that the recirculation motion will be weaker, and the convective motion will be slower (note this is also the case for the sensitivity of the frequency $RW_3$ to steady forcing figure 6d). On the other hand, the high-frequency mode is most sensitive in a near wake region behind the sphere, close to the recirculation bubble (see figure 6c,d). In this case, a forcing that decelerates the clockwise motion within the recirculation region would cause the largest stabilization effect.

Figure 6. Steady forcing at codimension point $A$. Sensitivity of amplification rate to the steady axisymmetric forcing $\boldsymbol {\nabla }_{\boldsymbol {H}_{0}} \lambda _{\ell }$ for (a) the low-frequency mode, also known as $RW_1$, and (c) the high-frequency mode, also known as $RW_3$. Sensitivity of the frequency to the steady axisymmetric forcing $\boldsymbol {\nabla }_{\boldsymbol {H}_{0}} \lambda _{\ell }$ for (b) the low-frequency mode, $RW_1$, and (d) the high-frequency mode, $RW_3$. The magnitude of the growth rate and frequency sensitivities is pictured by colours and their orientation by arrows.

Second, let us consider the receptivity of the flow to the presence of localized feedbacks, as in Giannetti & Luchini (Reference Giannetti and Luchini2007). The harmonic forcing $\boldsymbol {H} \equiv \boldsymbol {H}_{(z_{\ell })} \exp ({\mathrm {i} (\omega _{\ell } t + m_{\ell } \theta )})$ is defined as

(3.6)\begin{equation} \boldsymbol{H}_{(z_{\ell})} = \delta(\boldsymbol{x} - \boldsymbol{x}_0) \boldsymbol{C}_{(z_{\ell})} \boldsymbol{\cdot} \hat{\boldsymbol{u}}_{(z_{\ell})}, \quad \ell = 1,2,3, \end{equation}

where $\boldsymbol {C}_{(z_{\ell })}$ is a generic feedback matrix and $\delta (\boldsymbol {x} - \boldsymbol {x}_0)$ is the Dirac distribution centred at the point $\boldsymbol {x}_0 = (z_0, r_0, \theta _0)$. Thus, the variation of the eigenvalue due to the introduction of the localized feedback is

(3.7)\begin{equation} \Delta^{u} \mathrm{i} \omega_{\ell} = \langle {\hat{\boldsymbol{q}}_{(z_{\ell})}^{{{{\dagger}}}}} , {\delta \boldsymbol{H}_{(z_{\ell})} } \rangle = \boldsymbol{C}_{(z_{\ell})} : {\boldsymbol{\mathsf{S}}}_s^{({\ell})} (\boldsymbol{x}_0) , \quad \ell \in \textrm{I} \end{equation}

The rank two tensor of (3.7) is commonly designated as the structural sensitivity tensor, here denoted as ${\boldsymbol{\mathsf{S}}}_s^{({\ell })}$,

(3.8)\begin{equation} {\boldsymbol{\mathsf{S}}}_s^{({\ell})} \equiv \hat{\boldsymbol{u}}_{(z_{\ell})}^{{{{\dagger}}}} \otimes \bar{\hat{\boldsymbol{u}}}_{(z_{\ell})}, \quad \ell = 1,2,3. \end{equation}

The spectral norm of the structural sensitivity tensor for low- and high-frequency modes is depicted in figure 7. Similar to the receptivity to axisymmetric steady forcing, the recirculation bubble ($RW_3$) and the leftmost end of the recirculation region ($RW_1$) are the most sensitive regions of the flow.

Figure 7. Spectral norm of the structural sensitivity tensor. (a) Low-frequency mode, (b) high-frequency mode.

4.1. Classification of solutions

In the following, the right-hand side of (4.6) is designated $\boldsymbol {f}(\boldsymbol {z})$ where $\boldsymbol {z} = (z_1,z_2,z_3)$. The reduced vector $\boldsymbol {f}$ is equivariant under the action of the group $\varGamma \equiv SO(2) \times \mathbb {T}^{3}$, with the following action representation:

(4.8)\begin{equation} \left. \begin{aligned} \theta \boldsymbol{\cdot} \boldsymbol{z} & \equiv (z_1 {\rm e}^{{\rm i} l\theta},z_2 {\rm e}^{{\rm i} r\theta},z_3 {\rm e}^{{\rm i} s\theta}),\\ (\psi_1,\psi_2,\psi_3) \boldsymbol{\cdot} \boldsymbol{z} & \equiv (z_1 {\rm e}^{{\rm i} \psi_1},z_2 {\rm e}^{{\rm i} \psi_2},z_3 {\rm e}^{{\rm i} \psi_3}). \end{aligned} \right\} \end{equation}

Here $l,r,s \in \mathbb {Z}$, $\theta \in [0, 2{\rm \pi} )$ and $\psi _i \in [0,2{\rm \pi} )$ for $i = 1,2,3$; $(\psi _1,\psi _2,\psi _3)$ and $\theta$ are the representations in $\mathbb {C}^{3}$ of the actions of the group $\varGamma$, which correspond to the time shift and rotational invariance, respectively. The substitution of the polar decomposition of $\boldsymbol {z} = \boldsymbol {r} \textrm {e}^{\textrm {i} \boldsymbol {\varPhi } }$, with $\boldsymbol {r} = (r_1,r_2,r_3 )$ and $\boldsymbol {\varPhi } = ( \phi _1, \phi _2, \phi _3 )$, into (4.6) yields the following decoupled phase-amplitude system:

(4.9)\begin{equation} \left. \begin{aligned} \dot{{r}}_{\ell} & = r_{\ell} \left[ \varLambda^{R}_{\ell} + \mathcal{V}^{R}_{\ell k} {r}^{2}_{k} \right], \quad k,\ell = 1,2,3,\\ \dot{\phi}_{\ell} & = \varLambda^{I}_{\ell} + \mathcal{V}^{I}_{\ell k} {r}^{2}_k, \quad k,\ell = 1,2,3. \end{aligned} \right\} \end{equation}

Here $\varLambda = \varLambda ^{R} + \varLambda ^{I} \equiv (\lambda _1,\lambda _2,\lambda _3)^\textrm {T}$ and the matrix $\mathcal {V} = \mathcal {V}^{R} + \textrm {i} \mathcal {V}^{I}$ is

(4.10)\begin{equation} \mathcal{V} \equiv \begin{pmatrix} \nu_{11} & \nu_{12} & \nu_{13} \\ \nu_{21} & \nu_{22} & \nu_{23} \\ \nu_{31} & \nu_{32} & \nu_{33} \end{pmatrix}. \end{equation}

To ease the presentation of the fixed-point solutions of (4.9), let us introduce the inverse of the linear operator $\mathcal {V}$, which can be written as

(4.11)\begin{equation} \mathcal{V}^{{-}1} = \dfrac{1}{\det{\mathcal{V}}} \begin{pmatrix} \det{\mathcal{V}}_{11} & \det{\mathcal{V}}_{21} & \det{\mathcal{V}}_{31} \\ \det{\mathcal{V}}_{12} & \det{\mathcal{V}}_{22} & \det{\mathcal{V}}_{32} \\ \det{\mathcal{V}}_{13} & \det{\mathcal{V}}_{23} & \det{\mathcal{V}}_{33} \end{pmatrix}, \end{equation}

where $\det {\mathcal {V}}_{k\ell }$ denotes the minor of the matrix $\mathcal {V}$, obtained by eliminating the line $k$ and the column $\ell$.

In the following, the notation $\dot {\boldsymbol {r}} = \boldsymbol {f}^{R}( \boldsymbol {r} )$ will be adopted to denote the amplitude equation of the nonlinear system (4.9). The remainder of this subsection will be devoted to the study of the three fixed-point solutions of (4.9).

The classification of the solutions of the generic triple-Hopf bifurcation interaction with $SO(2)$ symmetry is based on maximal isotropy subgroups of the group $\varGamma$. This technique predicts the existence up to tertiary bifurcations of fixed points of the complex normal form (4.6). These isotropy subgroups correspond to the symmetries of the solutions within the fixed-point subspace of each isotropy group (cf. table 2).

Table 2. Nomenclature and symmetry group of fixed-point solutions of the system (4.9).

In our discussion, we identify the subgroups of $SO(2) \times \mathbb {T}^{3}$. Each element in the group has the form

(4.12)\begin{equation} (\theta, \psi_1, \psi_2, \psi_3) \in SO(2) \times \mathbb{T}^{3}. \end{equation}

Using this notation, the subgroup $S(k,l,r,s)$ of $SO(2) \times \mathbb {T}^{3}$ is defined as

(4.13)\begin{equation} S(k,l,r,s) = \{ (k \theta, l \theta, r \theta, s \theta ) | \theta \in S^{1} \}. \end{equation}

The conjugacy classes of isotropy subgroups of $SO(2) \times \mathbb {T}^{3}$ are documented with the representative of the fixed-point subspace in polar coordinates and the number of incommensurate frequencies in table 2. Additionally, a graphical representation of the isotropy lattice is displayed in figure 8 in terms of the class representative of the fixed-point subspace.

Figure 8. Isotropy lattice of the triple-Hopf bifurcation.

Rotating waves correspond to the simplest non-trivial fixed point of (4.9), which in the original set of equations is a periodic solution. They arise as the result of a supercritical Hopf bifurcation of the steady state (named trivial state in table 2) and they may eventually bifurcate into mixed modes; the eigenvalues of rotating waves may be found in the first row of table 3. Mixed modes, defined in table 3, are the result of the interaction between two rotating waves. A mixed mode has a representative in the normal form with two non-zero amplitude terms, thus, they correspond to a $T^{2}$-quasiperiodic state in the original system of equations. These states may experience two kinds of bifurcations. They may lose stability in the transversal direction or within their own subspace, these two conditions are listed in table 3. Eventually, a bifurcation in the transversal direction of a mixed mode may be associated with the appearance of an interacting mixed mode ($IMM_{123}$) attractor. An interacting mixed mode corresponds to a $T^{3}$-quasiperiodic state in the original system of equations, and it is represented by three non-zero amplitude terms. However, $T^{3}$-quasiperiodic states are hardly observed in numerical simulations of dissipative systems, as it is the case of Navier–Stokes equations (2.1), instead a chaotic attractor is usually detected. A more exhaustive analysis of the unfolding of the triple-Hopf bifurcation is left to Appendix C.

Table 3. Defining equations and eigenvalues of the solutions of the polar third-order normal form (4.9).

4.2. Illustration of the procedure

Let us detail the procedure followed to compute the bifurcation scenario, a procedure that is also followed in § 5 for the determination of the parametric portrait. For the sake of simplicity, we first discuss the bifurcation diagram for a constant rotation rate $\varOmega = 1.75$ in terms of the amplitudes $(r_1, r_2, r_3)$. We would like to remind the reader that the amplitudes $(r_1, r_2, r_3)$ are representative of the kinetic energy of the velocity fluctuations, based on the normalization choice of § 3.1. First, we need to determine the coefficients of the normal form, listed in tables 4 and 5, following the procedure of Appendix A. Then, one may evaluate the linear and cubic coefficients of the normal form at $\varOmega = 1.75$ for a variable Reynolds number from the evaluation of (4.7). Please note that, for the evaluation of cross-diagonal cubic coefficients, the expansion parameter $\varepsilon _\varOmega ^{2} = \varOmega _c - \varOmega$ depends on the location of the critical rotation rate $\varOmega _c$, that is, to evaluate $\nu _{13}$ one evaluates $\varepsilon _{\varOmega,A}^{2} = \varOmega _{A} - \varOmega$ whereas to evaluate $\nu _{23}$ one evaluates $\varepsilon _{\varOmega,B}^{2} = \varOmega _{B} - \varOmega$. The diagonal cubic coefficients may be evaluated directly at the bifurcation point for every rotation rate $\varOmega$ as a function of $\varepsilon _\nu ^{2}$ or by considering the cubic coefficient of the nearest codimension-two point. In our procedure, we found good agreement with time-stepping simulations in the range $1 \leq \varOmega \leq 3$ if we consider $\nu _{11} = \nu _{11}^{A}$, $\nu _{22} = \nu _{22}^{B}$ and $\nu _{33} = \nu _{33}^{B}$; the consideration of $\nu _{33}^{A}$ induces a small deviation in the transition from $MM_{23}$ to $IMM_{123}$ of few units of the Reynolds number. The corresponding coefficients for $\varOmega = 1.75$ are listed in table 6 (Appendix A). Please note that the procedure illustrated herein corresponds to a method to determine the coefficients of the normal form; nonetheless, these coefficients can be estimated from numerical simulations as in Fabre et al. (Reference Fabre, Auguste and Magnaudet2008) or following a data-driven approach, cf. Callaham, Brunton & Loiseau (Reference Callaham, Brunton and Loiseau2022), Loiseau & Brunton (Reference Loiseau and Brunton2018) and Loiseau, Noack & Brunton (Reference Loiseau, Noack and Brunton2018). Bifurcation events are designated by their corresponding value of the Reynolds number $Re_{{state}_a}^{{state}_b}$, where $\textrm {state}_a$ stands for the simplest state that exists before the bifurcation and ${state}_b$ stands for the resulting state after the bifurcation. In addition, the notation $Re_{{state}_a}^{\sigma _k,s}$ indicates a bifurcation of the $\textrm {state}_a$ where the eigenvalue $\sigma _k$ ($k=1,2,3$) has changed sign, $s$ indicates stabilization and $u$ indicates the change from stable to unstable of the referring eigenvalue/eigenmode pair. In the following, there is only a bifurcation of this kind, the one associated to the mixed mode $MM_{12}$ that is stabilized/destabilized because of a change of sign of the eigenvalue in the transversal direction ($r_3$). Thus, we simplify the notation to $Re_{{state}_a}^{s}$ or $Re_{{state}_a}^{u}$.

Table 4. Linear coefficients of the normal form (4.9) evaluated at codimension-two points.

Table 5. Cubic coefficients of the normal form (4.9) evaluated at codimension-two points.

Table 6. Cubic coefficients of the normal form (4.9) evaluated at $\varOmega = 1.75$.

Figure 9 displays the bifurcation diagram, with Reynolds number as the control parameter for $\varOmega = 1.75$. There exist three primary bifurcations, i.e. bifurcations from the axisymmetric steady state, located at $Re_{TS}^{RW_1}$, $Re_{TS}^{RW_2}$ and $Re_{TS}^{RW_3}$, respectively. However, the $RW_2$ branch remains unstable all along the analysed interval. The first transition to occur is a supercritical Hopf bifurcation leading to the $RW_1$ solution, which is then followed by another supercritical Hopf bifurcation leading to $RW_3$. For the range of Reynolds numbers $Re_{RW_1} < Re < Re_{RW_3}^{MM_{13}}$, there exists a single stable attractor, which corresponds to the limit cycle associated with the solution $RW_1$. At $Re_{RW_3}^{MM_{13}}$ the $RW_3$ branch experiences a Neimark–Sacker bifurcation that results in the appearance of the mixed-mode solution $MM_{13}$. In the interval $Re_{RW_3}^{MM_{13}} < Re < Re_{RW_1}^{MM_{13}}$ both primary solutions ($RW_1$ and $RW_3$) are stable under any arbitrary perturbation and in addition they are connected by the unstable mixed mode $MM_{13}$, which is located on the separatrix of the basin of attraction of the two primary solutions, the phase portrait of this scenario is sketched in figure 9(b i). Eventually, the solution branch $MM_{13}$ terminates at $Re = Re_{RW_1}^{MM_{13}}$, which makes $RW_3$ the single attractor of the system for the interval $Re_{RW_1}^{MM_{13}} < Re < Re_{RW_3}^{MM_{23}}$. The $RW_3$ branch eventually bifurcates into the mixed-mode branch $MM_{23}$, which is a stable attractor within the interval $Re_{RW_3}^{MM_{23}} < Re < Re_{MM_{23}}^{IMM_{123}}$. The other primary branch, the unstable $RW_1$, undergoes another Neimark–Sacker bifurcation at $Re_{RW_1}^{MM_{12}}$ which results in the existence of the $MM_{12}$ branch, yet unstable for perturbations in the transversal direction of the mixed mode (in the $r_3$ direction). The $MM_{12}$ mixed-mode branch appears to be stable only within a small interval $Re_{MM_{12}}^{s} < Re < Re_{MM_{12}}^{u}$, where two bifurcations, which are associated to an instability in the transversal direction $r_3$, occur at the two limit values. We have employed $s$ and $u$ to denote the stable or unstable nature of the $MM_{12}$ regime. Thus, for $Re_{MM_{12}}^{s} < Re < Re_{MM_{12}}^{u}$, there is a second region with multiple stable attractors, which is schematically displayed in figure 9(b ii). The last bifurcation accounted by the normal form is the destabilization of the $MM_{23}$ branch at $Re = Re_{MM_{23}}^{IMM_{123}}$ that leads to the appearance of the $IMM_{123}$ branch, whose phase portrait is sketched in figure 9(b iii). Please note that despite the fact that $IMM_{123}$ is a fixed-point solution of the normal form, the Newhouse–Takens–Ruelle theorem indicates that the original system of equations may exhibit a chaotic attractors shadowing the $IMM_{123}$ solution.

Figure 9. (a) Transition scenario at ${{\varOmega }} = 1.75$. Attractors are depicted with solid lines, whereas unstable invariant states are represented with dashed lines.(b) Schematic representation of phase portraits. (i) Two stable rotating waves separated by a mixed-mode solution. (ii) Two stable mixed modes. (iii) An interacting mixed-mode attractor, the chaotic attractor that shadows the $IMM_{123}$ is sketched in a lighter blue colour.

5.1. Comparison with DNS

In this section we assess the validity of the normal form to characterize the bifurcation scenario, as well as its capability to predict accurately the frequencies and force coefficients of the flow. The estimations of the normal form are compared with DNS results, which are performed at constant rotation rate $\varOmega =1.75$, the scenario analysed in § 4.2. As a first guess, we show the accurate prediction of the fundamental frequencies of each of the invariant states from normal form analysis in figure 10(a) in comparison with DNS results (markers), which will be discussed below.

Figure 10. Frequency characterization of the flow at $\varOmega = 1.75$. (a) Frequency evolution estimated from the normal form (continuous lines), where attractors are represented with continuous lines, whereas unstable invariant solutions are depicted with dashed lines. The markers of figure (a) denote the resulting pattern obtained from a time-stepping simulation; axisymmetric steady state $\circ$, $RW_1$ $\vartriangleleft$, $RW_3$ $\vartriangleright$, $MM_{12}$ $\square$, $MM_{23}$ $\lozenge$, IMM ₁₂₃ $\star$, the family of initial conditions is visualized by colours (family I: grey, family II: blue). Figures (b–g) display the FFT fluctuating velocity spectra for the different regimes obtained by means of DNS. Two velocity components: $W'$ (red solid line), $U'$ (black solid line) and locations (0, 0, 3.5) and (0, 0.5, 3.5), respectively, are selected to characterize all the frequencies in the wake. Unstable mode frequencies are included: low frequency (light red solid line), medium frequency (light green dashed line) and high frequency (light blue dotted line). Results are shown for (b) $Re = 110$, (c) $Re = 110$, (d) $Re = 150$, (e) $Re = 150$, ( f) $Re = 181$, (g) $Re = 210$.

Direct numerical simulations have been carried out to confirm the existence of the bi-stability of full governing equations (2.1) at $Re=110$. In particular, two families of time-stepping simulations have been performed with two distinct initial conditions, in such a way that, after a transient period, each of them converged towards different time-periodic solutions. On the one hand, we have used a static, axisymmetric base flow initially obtained at $Re = 70$, $\varOmega = 0$ which is able to develop the $RW_{1}$ state (family I, grey markers). On the other hand, we have used the solution obtained by Lorite-Díez & Jiménez-González (Reference Lorite-Díez and Jiménez-González2020) at $Re = 250$, $\varOmega = 2.2$ as initial seed to find the $RW_{3}$ state (family II, blue markers), which in turn confirms the existence of multiple stable attractors at $Re = 110$.

To determine the observed regimes, we have computed the frequency components corresponding to $St_{RW_1}$ and $St_{RW_3}$, displayed in figure 10(b,c), using fast Fourier transform (FFT) spectra of pointwise streamwise and radial velocities in the near wake. The spectra are calculated using the oscillatory part of the velocity-time evolution, i.e. $U' = U-\bar {U}$ for the radial velocity component, where $\bar {\cdot }$ stands for the temporal averaging operator. The three-dimensional topology of the two rotating wave patterns are displayed by means of iso-surfaces of $Q$ quantity in figure 11, where single and double helices topologies are shown.

Figure 11. Three-dimensional structures of $RW_1$ (a) and $RW_3$ (b) at $Re=110$ and $\varOmega =1.75$. We have used isosurfaces of Q-criterion, $Q=0.001$, coloured by streamwise vorticity, $\varpi _{z} \in [-1,1]$ (blue to red) to depict the flow structure.

Similarly, the existence of two stable mixed-mode attractors ($MM_{12}$ and $MM_{23}$) is confirmed at $Re=150$. Two time-stepping simulations of the full governing equations at $\mbox { {Re}} = 150$ were performed, using $RW_{1}$ and $RW_{3}$ solutions as initial seeds, the resulting frequency spectra of these patterns are displayed in figure 10(d,e), where the different frequencies associated with $MM_{12}$ and $MM_{23}$ are identified. In particular, we can see that the appearance of the medium-frequency component originates quasiperiodic regimes. The mixed mode $MM_{12}$ has been detected only in a small interval of Reynolds numbers, which is faithfully captured by the normal form; however, the value of $Re^{s}_{MM_{12}} \approx 154$ slightly differs from the results of the DNS, which show a stable $MM_{12}$ for $Re=150$. The corresponding patterns are displayed in figure 12(a,b).

Figure 12. Three-dimensional structure of $MM_{12}$ (a), $MM_{23}$ (b) and $IMM_{123}$ (c) at $Re=150$ (a,b) and $Re=210$ (c) and $\varOmega =1.75$. We have used isosurfaces of Q-criterion, $Q=0.001$, coloured by streamwise vorticity, $\varpi _{z} \in [-1,1]$ (blue to red) to depict the flow structure.

The $T^{3}$-quasiperiodic state $IMM_{123}$ has been detected with DNS for Reynolds numbers $Re \approx 181$. Such a state seems to be the single stable attractor in the analysed range $181 < Re < 210$. A series of DNS were carried out with two families of initial conditions: the mixed modes $MM_{12}$ and $MM_{23}$, both obtained at lower Reynolds numbers. Eventually, every DNS converged to the $IMM_{123}$ state, which seems to confirm the claim that it is the single stable attractor. Their associated spectra is depicted in figure 10( f,g) and its complex topology can be seen in figure 12(c). The identification of the three main frequencies (low, medium and high) in the spectra, figure 10( f), along with the multiple nonlinear interactions between them is possible for Reynolds number values near the bifurcation value $Re^{IMM_{123}}_{MM_{23}}\simeq 181$. However, it rapidly departs from the $T^{3}$-quasiperiodic state towards a more irregular state ($Re=210$) with a nearly continuous velocity spectrum, depicted in figure 10(g).

Globally, the good agreement between normal form analysis and DNS frequencies shows the predictive capability of the normal form within the range of Reynolds numbers studied (see figure 10a).

A further investigation of the dynamics of this attractor has been carried out by means of the $0-1$ test. Such a test was introduced by Gottwald & Melbourne (Reference Gottwald and Melbourne2004, Reference Gottwald and Melbourne2009) to distinguish between regular and chaotic dynamics. More precisely, it corresponds to a dichotomy test where an estimate $K$, associated to an asymptotic growth of the dynamics, takes discrete values $[0,1]$ which are associated with non-chaotic ($0$) and chaotic ($1$). Further details are given in Appendix F. The results corresponding to the application of the test to the local radial velocity $U'(0,0.5,3.5)$, obtained for the two families of computations, indicate a rapid departure towards chaotic dynamics, displayed in the Appendix (figure 16). These results confirm that the transition scenario is eventually ended by the Newhouse–Takens–Ruelle route to chaos.

Furthermore, DNS results can be used to illustrate the spatial pattern associated with each fundamental frequency (low, medium, high). To that aim, a high-order dynamic mode decomposition (HODMD) technique (see Vega & Le Clainche Reference Vega and Le Clainche2020, and references therein) has been applied to instantaneous fields of streamwise vorticity $\varpi _z$ located at $z = 2.5$ to isolate the spatial distributions of the main frequency components. More details about the HODMD technique and its application to present data can be found in Appendix E. In particular, the application of the technique to flow patterns $MM_{12}$, $MM_{23}$ and $IMM_{123}$ allowed the spatial characterization of the fundamental frequencies and their interactions. Apart from these frequencies, the methodology also provides the approximate contribution of each fundamental mode in the nonlinear state. The spatial patterns identified by HODMD are depicted using contours of spanwise vorticity without normalization. Thus, for the $MM_{12}$ state, the HODMD decomposition identifies four energetic modal contributions, corresponding to frequencies $St=0, 0.103, 0.174$ and $0.072$, which are depicted in figure 13(a). The use of instantaneous snapshots of $\varpi _z$ allows identifying mean flow ($St=0$) given by the constant streamwise rotation. Likewise, the low-frequency component ($St=0.103$) displays a dipole $m=-1$ topology. Similarly, the medium-frequency component ($St=0.174$), also features a $m=-1$ structure, although their topology resembles the Yin-Yang mode (Auguste, Fabre & Magnaudet Reference Auguste, Fabre and Magnaudet2010). Finally, an axisymmetric $m=0$ topology is identified at $St=0.072$, as the product of the interaction between low-frequency (LF) and medium-frequency (MF) modes, $St_{m_0} \simeq St_{MF}-St_{LF}$. Such a mode is specially energetic close to the longitudinal axis and even dominant in some locations. In Lorite-Díez & Jiménez-González (Reference Lorite-Díez and Jiménez-González2020) the authors identified such a mode as a fundamental frequency of the flow, named $f_b$ therein. However, given the results from linear stability and normal form analysis, we have now identified this mode as a subproduct.

Figure 13. Nonlinear patterns identified through HODMD analysis. The patterns are depicted using streamwise vorticity contours, $\varpi _z$. Results are shown for (a) $MM_{12}$, (b) $MM_{23}$ and (c) $IMM_{123}$ regimes. Tags values inside the contours refer to the corresponding frequency of each mode, $St_i$.

The same analysis for the $MM_{23}$ regime allows the identification of three main frequency components, depicted in figure 13(b). Apart from the mean flow ($St=0$), the flow decomposition pinpoints a mode with high frequency (HF) $St_{HF}=0.304$, that displays a $m=-2$ structure. In addition, the Yin-Yang mode is also retrieved but now with a frequency, $St_{MF}=0.238$. Note that, the amplitude associated with this medium-frequency component is very small. The weak energy associated with such a medium-frequency mode in the $MM_{23}$ regime is also observable in the corresponding spectra in figure 10(e).

Furthermore, the HODMD analysis of the complex regime $IMM_{123}$, present at $Re= 210$ and $\varOmega =1.75$, is depicted in figure 13(c). The decomposition identifies the axisymmetric mean flow, the three fundamental frequencies and many interactions between them, although only five energetic components have been selected for depiction. For instance, the low-frequency and medium-frequency modes display $m=-1$ symmetries and respective frequencies $St_{LF}=0.095$ and $St_{MF}=0.231$, which are similar to those corresponding to the $MM_{12}$ and $MM_{23}$ regimes. Similarly, the frequency value of the high-frequency mode ($m=-2$) remains nearly constant with respect to previous regimes, $St_{HF}=0.299$, indicating that the dependence of frequencies with $Re$ is small (as in figure 3b). Among the main subproducts, the axisymmetric pattern ($m=0$) produced by the interaction of fundamental frequencies $St_{m0}\simeq St_{LF}+St_{HF}-St_{MF}=0.163$ stands out. It should be noted that there is a small mismatch between the identified peaks in FFT and the frequencies obtained by HODMD that is produced by the different sampling period and the corresponding recording frequency during the simulation.

To complement the previous analysis, we next focus on the effect of the different flow states, shown in figure 11 and figure 12, on the sphere's aerodynamic forces. Thus, we present in figure 14(a) the evolution with respect to Reynolds number at $\varOmega =1.75$ of the time-averaged drag coefficient, $\overline {C_{D}}=\overline {C_{z}}$, and on figure 14(b) the total transverse force coefficient $\overline {C_L}=\overline {\sqrt {C^{2}_{x}+C^{2}_{y}}}$ for normal form analysis and DNS.

Figure 14. Evolution of the time-averaged forces with respect to Reynolds number for a constant rotation rate $\varOmega =1.75$. The force coefficients determined from the normal form analysis are represented by continuous lines. The force coefficients determined from DNS are depicted with markers. Same legend as in figure 10.

The comparison between force coefficients obtained by means of normal form and DNS approaches is indeed fairly satisfactory. The method captures the main trends in the forces and, for most of the states, the prediction is reasonably similar, as it was with the different fundamental frequencies. More precisely, for simpler regimes, such as rotating waves, the normal form analysis and DNS provide the same results. At $Re = 150$, for mixed modes, both methods display a small discrepancy (below five percent). This seems to be related to a small nonlinear contribution of the medium-frequency mode identified by the DNS, as it can be seen in the corresponding FFT spectra (see figure 10d,e). That said, the general trend of the mean drag $\overline {C_D}$ displays a general reduction with $Re$, which may be partly due to a smaller viscous drag contribution. Besides, the pressure drag component is likely to decrease as well, since $L_r$ is shown to increase with $Re$ (see figure 5b) for a given rotation rate, $\varOmega$. In particular, as discussed by Roshko (Reference Roshko1993) for three-dimensional bluff body wakes, an increase of the recirculating length leads to a decrease in the drag values on account of a pressure recovery associated to changes in the curvature of the separatrix line. Additionally, major changes in the trends are reported between different flow regimes, as expected from strong modifications in the near wake topology and flow separation (Lorite-Díez & Jiménez-González Reference Lorite-Díez and Jiménez-González2020).

Similarly, the agreement is also good for the mean total transverse coefficient, $\overline {C_L}$, although it displays small deviations for states $RW_1$ and $MM_{23}$ (see figure 14b). The value of $\overline {C_L}$ is strongly affected by the wake regime and the corresponding azimuthal symmetry. Therefore, the two families of simulations display quite a different evolution. Moreover, it should be noted that the high-frequency mode does not create a net component of transverse force due their symmetric wake topology, as it is seen for $RW_{3}$ and $MM_{23}$ regimes, where this mode is dominant, causing $\overline {C_L}\simeq 0$. In those cases, the wake structures net eccentricity is small, inducing a negligible transverse force. In view of such results, such regimes should be favoured in case of control if a stabilization of the trajectory is wished, e.g. for freely rising or falling rotating spheres, as the transverse displacement of the body might be limited. Conversely, $RW_1$, $MM_{12}$ and $IMM_{123}$ are likely to cause lateral shift and destabilization of the trajectory for freely moving bodies due to their greater mean lateral force and their corresponding eccentric wake structures.

5.2. Parametric exploration

Let us discuss the influence of the rotation rate in the dynamics of the flow past the rotating sphere. For that purpose, we determine the stable attractors of the normal form in the range $\mbox { {Re}} < 300$ and ${{\varOmega }} < 4$. The cubic normal form coefficients are determined following the same procedure as in § 4.2. However, for low rotation rates ($\varOmega < 0.8$), we found that the linear coefficients of the normal form were not correctly estimated. So, for these values of the rotation rate, the linear coefficients are determined exactly at the threshold of instability of each codimension-one bifurcation.

The flow past the static sphere ($\varOmega = 0$), analysed by Fabre et al. (Reference Fabre, Auguste and Magnaudet2008), experiences a symmetry breaking bifurcation that leads to the reflection-symmetric bifid wake (steady-state mode) and eventually a Hopf bifurcation that leads to the RSP. The phase diagram depicted in figure 15 shows that both the rotating wave $RW_1$ and the mixed mode $MM_{12}$ are the continuation in the parameter space of the steady-state and RSP mode, respectively. Thus, dynamics for low values of the rotation rate $(\varOmega < 1)$ are qualitatively similar to the flow past the static sphere. However, for rotation rate values slightly larger than $\varOmega _B$, one starts detecting a wake with a double helix. This fact has been evidenced by the numerical simulations of Lorite-Díez & Jiménez-González (Reference Lorite-Díez and Jiménez-González2020), Pier (Reference Pier2013) and the experimental work of Skarysz et al. (Reference Skarysz, Rokicki, Goujon-Durand and Wesfreid2018) who reported that at around ${{\varOmega }} \approx 1.5$ for large Reynolds numbers ($Re > 250$) the quasiperiodic motion of the wake changes from a single to double helix pattern, which is consistent with the phase diagram in figure 15. However, if we look in detail at how the flow transits from a single helix to the double helix wake in terms of the rotation of the sphere, one can observe three bi-stable regions, between $RW_{1}$ and $RW_{3}$, between $RW_{3}$ and $MM_{12}$, and between $MM_{12}$ and $MM_{23}$. These bi-stable sections connect two regions of the parameter space with distinct attractors, the rotating waves $RW_{1}$ and $RW_{3}$ and the mixed modes $MM_{12}$ and $MM_{23}$. It is of interest that the formation of these new states occurs near the codimension-two point $B$, which exhibits the importance of this bifurcation as an organizing centre, even though it occurs as a bifurcation for an already unstable trivial state. Furthermore, the importance of the codimension point $A$ is clearly evidenced by figure 15, which also acts as an organizing centre of dynamics. Around point $A$, one finds four distinct regions with inequivalent dynamics. If we move counterclockwise from point $A$, we have the left region where the trivial axisymmetric state is stable, the lower region where the $RW_1$ state is the single attractor, a region with two stable attractors ($RW_1$ and $RW_3$) and the upper region where the $RW_3$ is the single stable state. Finally, the significance of the other codimension-two bifurcation, the point $C$, is rather more subtle. This point is located in the only region where one can observe the mixed mode $MM_{13}$, which is only observed for very large rotation rates, and it connects the $RW_3$ state and the $T^{3}$-quasiperiodic state $IMM_{123}$.

Figure 15. Phase diagram of the nonlinear patterns in the range $\varOmega < 4$ and $\mbox { {Re}} < 300$, as predicted by the normal form (4.9). The axisymmetric state (TS) persists in the white region. Shaded regions indicate the existence of a stable pattern (respectively, stable patterns). Dashed lines illustrate unstable rotating waves neutral curves obtained by LSA.

Therefore, one may conclude that the rotation of the sphere has a mild effect on the bifurcation scenario for low rotation rates ($\varOmega < 1$). Rotation rates between $1 < \varOmega < 2$ favourise the appearance of a double helix wake and hysteresic behaviour, whereas large rotation rates ($\varOmega > 2$) have a destabilizing effect which rapidly triggers the emergence of chaotic dynamics via a Ruelle–Takens–Newhouse route.