4.1. Correlation between bursters and topology of marginal stability curves

Based on the linear growth rate of axisymmetric disturbances, i.e. $m = 0$ $(k_c = 0)$ , Riahi et al. (Reference Riahi, Hayani Choujaa and Aniss2025) have noticed the presence of bursters near codimension-two bifurcation points where more than one instability mode co-exist and compete at the critical onset of the primary bifurcation. The presence of these multitude of instability modes gives rise to rough growth rate evolution versus the axial wavenumber, which is not the case when dealing with regular instability modes where unimodal and smooth growth rate evolutions are observed. Does the burster nature of the primary instability mode alter the topology of the marginal stability curves ? This is the main question that we aim to highlight in this subsection.

Figure 2.

Marginal stability curves computed using Floquet theory showing the apparition of new flow states and bicritical states under the effect of the frequency $\gamma _1$ in the case where $\omega = \sqrt {3}$ . On each curve, the vertical axis is the Taylor number $\mathit{Ta}$ while the horizontal one is the dimensionless axial wavenumber $q$ , ranging from 3 to 5.5.

For this purpose, marginal stability curves are depicted in figure 2 at $\omega = \sqrt {3} \approx 97/56$ and different values of the frequency $\gamma _1$ . As one notices, only smooth, regular and uni-modal marginal stability curves are obtained when $\gamma _1 \succcurlyeq 4$ in contrast to the case $\gamma _1 = 3.4$ where finger-like features appear on a low wavenumber branch. The latter becomes more evident and parabolic at $\gamma _1 = 3.155$ where a codimension-two bifurcation point is detected, corresponding to a bicritical state between smooth and rough branches. Indeed, the smooth branch corresponds to a regular non-reversing instability mode, whose numerical characterisation has been carried out recently by Riahi et al. (Reference Riahi, Hayani Choujaa and Aniss2025) at $\gamma _1 = 10$ , whereas the rough one corresponds to a family of bursting modes. Such roughness is a physical feature that indicates the existence of numerous instability modes close to each other, clustered in a specific marginal stability branch. Such roughness concerns also the growth rate of the perturbations, as was recently shown by Riahi et al. (Reference Riahi, Hayani Choujaa and Aniss2025) when dealing with the bursting mode discussed at $\gamma = 3.1$ , which belongs to the rough branch in the marginal stability curve presented in figure 2. This feature sheds light on the existence of a strong correlation between the topology of the marginal stability curves and the nature of the bifurcated solution, whether it corresponds to a regular or bursting mode: smooth marginal stability curves correspond to regular instability modes whereas bursting modes correspond to rough marginal curves. It is interesting to note here that such a bursting phenomenon requires a very small wavenumber increment to be adequately elucidated because even a small uncertainty in the selection of the critical wavenumber could lead to incorrect identification of the most dangerous mode on rough branches of marginal stability curves.

Figure 3.

Marginal stability curves computed using Floquet theory, showing the apparition of new flow states and bicritical states under the effect of the frequency $\gamma _1$ in the case where $\omega = 1/\sqrt {37}$ . On each curve, the vertical axis is the Taylor number $\mathit{Ta}$ while the horizontal one is the dimensionless axial wavenumber $q$ , ranging from 6 to 14.

A further decrease in the frequency parameter $\gamma _1$ suppresses the smooth instability branch, while other rough branches appear. The dynamics of the system becomes subsequently rich and complex in the sense that each rough instability branch is composed of numerous sub-branches and consequently several codimension-two bifurcation points arise. This results in discontinuities in the axial wavenumber evolution versus the frequency number. The frequency step in these numerical calculations has to be slowly and carefully varied owing to the strong competition between numerous harmonic and subharmonic instability modes.

For the case $\omega \prec 1$ and by considering $\omega = 1/\sqrt {37} \approx 12/37$ , strong mode competition is observed at sufficiently high frequencies. This feature is illustrated in figure 3, presenting examples of marginal stability curves obtained at $\omega = 1/\sqrt {37}$ and some representative values of the frequency $\gamma _1$ . As one can see, two smooth instability branches emerge at $\gamma _1 = 10$ . The most unstable mode is a regular non-reversing mode, as demonstrated by Riahi et al. (Reference Riahi, Hayani Choujaa and Aniss2025) at this frequency, reached by crossing $\mathit{Ta}_c = 511.868$ and $q_c = 10.3$ while the second one is a bursting mode. This bursting mode is represented by a smooth, narrow portion of a stability curve belonging to a rough stability branch. Such roughness is evident as the oscillation frequencies of the cylinders, i.e. $\gamma _1$ and $\gamma _2$ , diminish. In this case, the first high-frequency codimension-two bifurcation point occurs at $\gamma _1 = 8.75$ where the primary bifurcation switches from a harmonic regular non-reversing Taylor vortex flow emerging in the frequency interval $\gamma _1 \succ 8.75$ to a bursting harmonic reversing Taylor vortex flow emerging in the frequency interval $8.35 \prec \gamma _1 \prec 8.75$ . In this case, the non-reversing frequency range shrinks considerably in comparison with the case $\omega = \sqrt {3}$ , i.e. $\gamma _1 \succ 8.75$ against $\gamma _1 \succ 3.2$ and against $\gamma _1 \succ 1.3$ when $\omega = 1$ . The second high-frequency codimension-two bifurcation point occurs at $\gamma _1 = 8.35$ below which only rough marginal stability curves with both reversing harmonic and subharmonic bursting solutions are observed.

Such roughness is a signature of strong mode competition between different and numerous harmonic and subharmonic modes appearing in very narrow wavenumber ranges. To clearly elucidate this feature, two examples of marginal stability curves are depicted in figure 4 for two subsequent codimension-two bifurcation points detected at $\gamma _1 = 8.32$ and $\gamma _1 = 8.26$ . The first codimension-two bifurcation point is observed by crossing $\gamma _1 = 8.32$ and $\mathit{Ta}_c = 417.444$ , where the system loses stability via two instability modes bifurcating simultaneously from the basic state at two critical wavenumbers, i.e $q_c = 7.02$ and $q_c = 7.25$ for the subharmonic (period-doubling) and (synchronous) harmonic modes, respectively. This feature is well illustrated in figure 5, that shows the calculated Floquet multipliers in the complex plane at the first codimension-two bifurcation point.

Figure 4.

Marginal stability curves computed using Floquet theory showing two bicritical states between harmonic and subharmonic responses in the case where $\omega = 1/\sqrt {37}$ at $\gamma _1 = 8.32$ and $\gamma _1 = 8.26$ . On each curve, the vertical axis is the Taylor number $\mathit{Ta}$ while the horizontal one is the dimensionless axial wavenumber $q$ , ranging from 6 to 8.5.

Figure 5.

Floquet multipliers at the bicritical state between harmonic and subharmonic instability modes reached at $\mathit{Ta}_c = 417.444$ , $q_c = 7.02$ , $q_c = 7.25$ when $\omega = 1/\sqrt {37}$ and $\gamma _1 = 8.32$ .

4.2. Burster identification strategy

Among the key characteristics of bursting behaviour is the slow and fast dynamics where the solution evolves slowly in some regions and suddenly jumps to a fast oscillatory state. This could be also observed from the phase portrait characterised by loops and clusters, i.e. a repetitive pattern where the system alternates between quiescent and oscillatory states.

Figure 6.

(a) Time series of the radial perturbation velocity and (b) its corresponding PSD recorded at the point $(x, \theta , z) = (1/2, 0, 0)$ over $2$ oscillation periods ( $\tau ' \in [4T, 6T]$ ) in the case where $\omega = 1/\sqrt {37}$ , $\gamma _1 = 10$ . The critical instability parameters are $\mathit{Ta}_c = 584.66$ and $q_c = 6.1$ .

In order to confirm the bursting nature of the left-hand axisymmetric instability mode observed in the marginal stability curve at the sufficiently high frequency $\gamma _1 = 10$ (figure 3), time series of the radial velocity and its corresponding power spectral density (PSD) at the point $(x, \theta , z) = (1/2, 0, 0)$ over $2$ oscillation periods ( $\Delta \tau ' = 2\pi$ ) are plotted in figure 6. It is evident that this time series is $\pi \hbox{-}$ periodic and exhibits a change of sign between negative and positive values. This behaviour means that we are dealing with a harmonic reversing flow characterised by sharp peaks, indicating fast oscillations superimposed on a slower modulation. There is no clear alternation between quiescent (low-amplitude) phases and sudden bursts as in classical bursting systems. However, the waveform suggests the presence of eigenmode interactions leading to sharp spikes. The PSD indicates that we are not dealing with a purely $\pi $ -periodic response, owing to the strong presence of additional frequency components. Also, the existence of sharp peaks suggests that the system is likely undergoing modulated periodic motion rather than fully chaotic bursting.

A set of phase portraits providing more insights into this instability mode is depicted in figure 7. Figure 7(a) represents a three-dimensional phase space using the $(u_r, u_\theta , u_z)$ velocities where it is noticed that the trajectory evolves in a confined, structured manner and exhibits oscillatory behaviour with slow and fast dynamics, a hallmark of bursting instabilities. Projection in the $(u_r, u_\theta )$ phase space is depicted in figure 7(b), which highlights the complex nature of the trajectory, suggesting a slow evolution around one region, followed by a rapid excursion (for a detailed version of this trajectory coloured by time along with time markers, readers should refer to figure 8). This is typical in slow–fast systems, reinforcing the idea of a bursting mode. A projection in the phase space $ (u_r, ( {{\rm d} u_r}/{{\rm d}\tau '}) )$ is illustrated in figure 7(c) showing a butterfly-shaped structure describing thus a classical bursting phase portrait. The latter sheds light on the trajectory that expands and contracts, forming two large lobes on both sides of $u_r = 0$ . Near this point, i.e. $u_r = 0$ , there are dense, chaotic-like oscillations, indicating a metastable state. This suggests that the system spends some time in a relatively slow phase before a burst occurs. In addition, this phase portrait shows fast excursions to large amplitudes where the system abruptly jumps to large values of $ {{\rm d} u_r}/{{\rm d}\tau '}$ , indicating thus the strong presence of rapid bursting events. These are characteristics of a slow–fast dynamics, where the system alternates between slow motion and rapid jumps (bursts). The projected trajectory in the $ (u_\theta , ( {{\rm d} u_\theta}/{{\rm d}\tau '}))$ space (figure 7 d) aligns with these conclusions and shows a rapid transition from a low-amplitude regime to a high-amplitude one.

Figure 7.

Phase portraits illustrating the bursting instability mode obtained at the mid-gap $(x, \theta , z) = (1/2, 0, 0)$ when $\omega = 1/\sqrt {37}$ , $\gamma _1 = 10$ , $\mathit{Ta}_c = 584.66$ and $q_c = 6.1$ during the time interval $\tau ' \in [4T, 5T]$ . (a) Three-dimensional trajectory in the $(u_r, u_\theta , u_z)$ space, highlighting the oscillatory dynamics. (b) Projection onto the $(u_r, u_\theta )$ plane, showing the complex interplay between the radial and azimuthal velocity components. (c) Phase portrait in the $(u_r, ({{\rm d}u_r}/{{\rm d}\tau '}))$ space, revealing characteristic slow–fast transitions indicative of bursting behaviour. (d) Projection onto the $(u_\theta , ( {{\rm d}u_\theta }/{{\rm d}\tau '}))$ space, emphasising the role of azimuthal velocity variations. Parameters as in figure 6. The red and green markers represent respectively the final and initial instants.

Figure 8.

(a) Time series of the radial perturbation velocity recorded at the point $(x, \theta , z) = (1/2, 0, 0)$ over one oscillation period ( $\tau ' \in [4T, 5T]$ ). (b) Time-coloured projection of the phase-space trajectory with equidistant time markers, highlighting the fast (sparser) and slow (denser) dynamics consistent with the temporal signal during the time interval $\tau ' \in [13.08, 13.25]$ . Parameters as in figure 6. The denser region (region 1) of the trajectory corresponds to slow evolution where the system remains temporarily trapped in an oscillatory regime in contrast to the sparser region (region 2) that reflects a fast transition between two oscillatory states as indicated by the steep slope in the time series signal.

From a physical point of view, bursting modes occur when interactions between two oscillatory Stokes layer flows generated near each cylinder are strongly involved (Riahi et al. Reference Riahi, Hayani Choujaa and Aniss2025). It is interesting to recall here that, at these frequencies, i.e. $\gamma _1 = 10$ and $\gamma _2 = \gamma _1/\sqrt {\sqrt {37}} \approx 0.405\gamma _1 \approx 4.05$ , a purely oscillatory Stokes layer flow sets in near the inner cylinder, meaning that this oscillatory flow generates a damped viscous wave whose penetration depth does not exceed the Stokes layer thickness $\delta _1 = 0.1d$ . Meanwhile, the penetration depth of the damped viscous wave generated by the outer cylinder, i.e. $\delta _2 = 2.46d$ , is greater than $\delta _1$ and requires more than the gap width to be attenuated. This feature suggests that we are dealing with a flow configuration mainly controlled by the outer cylinder drive. The veracity of these speculations is substantiated and clearly illustrated in snapshots of two-dimensional (2-D) and 3-D contour plots of the radial component of the perturbed velocity field depicted in figure 9 over one oscillation period. As one can see, the gap between the cylinders is filled by pairs of large Taylor vortices rotating in a privileged direction (clockwise say for vortices located at a representative vertical position, e.g. $z = 0.5$ ) at $\tau ' = 0$ with small vortices, having opposite polarity (anti-clockwise say at $z = 0.5$ ), near the outer cylinder. These small outer-cylinder-driven vortices undergo a growth process over time and transport, towards the inner cylinder, the large vortices that become smaller at $\tau ' = T/8$ . Other outer-cylinder-driven vortices, rotating in the clockwise direction, emerge and become larger just after $\tau ' = T/8$ , suppressing thus the existing vortices. This scenario is repeated several times during one period and highlights the reversing flow process in this parameter space. The burster is generated by interactions, taking place near the inner cylinder, between the rapid and repetitive small outer-cylinder-driven vortices and the inner cylinder drive. The latter pushes these vortices inwards (see the instants $\tau ' = 15/16T$ and $\tau ' = 4/16T$ ) or outwards depending on the instantaneous order of magnitude of the drives.

Figure 9.

Snapshots of the radial component of the perturbed velocity field relative to the bursting instability mode using two visualisation modes. The left panel displays a 2-D contour plot, while the right panel presents a 3-D isosurface rendering. These snapshots have been taken after transients have vanished. In particular, a time shift of four periods $t \to t - 4T$ is applied. Parameters as in figure 6.

It is interesting to note here that this bursting mode exists at $\gamma _1 = 10$ but it is considered as the primary bifurcation of the system just within the frequency interval $8.35 \prec \gamma _1 \prec 8.75$ where the Stokes layer thicknesses $\delta _1$ and $\delta _2$ become larger, enhancing thus the bursting phenomenon. The strong contribution of the outer cylinder drive is evident in this bursting mode in contrast to the most unstable mode observed at this frequency, i.e. $\gamma _1 = 10$ , where the Taylor vortex flow becomes mainly driven by the inner cylinder. At $\gamma _1 = 10$ , the resulting flow corresponds to a regular oscillatory Stokes layer flow adjacent to the inner cylinder, as is well depicted in figure 10 using the same snapshots as in figure 9. This means that no interactions between the Stokes layer flows generated by the cylinders are expected when dealing with regular instability modes.

Phase portraits equivalent to those used in figure 7 are illustrated in figure 11 and are related to the regular instability mode observed at $\gamma _1 = 10$ . These phase portraits suggest a different type of dynamic behaviour compared with the previous bursting one. Indeed, this figure shows more structured, layered oscillations, especially in panels (a) and (b), which might indicate transitions between different periodic states. In panel (c), it is evident that we are dealing with a superposition of numerous and different limit cycles characterised by different frequencies. This behaviour aligns quite well with the multi-frequencies observed in the PSD of the time series of the local radial velocity, illustrated by Riahi et al. (Reference Riahi, Hayani Choujaa and Aniss2025) for this instability mode in this parameter space.

Figure 10.

Snapshots of the radial component of the perturbed velocity field relative to the regular instability mode using two visualisation modes. The left panel displays a 2-D contour plot, while the right panel presents a 3-D isosurface rendering. These snapshots have been taken after transients have vanished at $\omega = 1/\sqrt {37}$ , $\gamma _1 = 10$ , $\mathit{Ta}_c = 511.8685$ and $q_c = 10.3$ . In particular, a time shift of four periods $t \to t - 4T$ is applied.

Figure 11.

Phase portraits illustrating the regular instability mode obtained at the mid-gap $(x, \theta , z) = (1/2, 0, 0)$ when $\omega = 1/\sqrt {37}$ , $\gamma _1 = 10$ , $\mathit{Ta}_c = 511.8685$ and $q_c = 10.3$ during the time interval $\tau ' \in [4T, 5T]$ . Panels are similar to those used in figure 7.

4.3. Bursters at different approached values of the frequency ratio $\omega$ : the case $\omega = \sqrt {3}$

A strong relationship exists between the computational costs of the numerical calculations and the approached values assigned to the frequency ratio parameter $\omega$ in the sense that, the higher the integers values of $\tilde {p}$ and $\tilde {q}$ in the stability problem, the longer calculation time required. The subsequent investigation seeks to ascertain the impact of the designated values for the integer parameters $\tilde {p}$ and $\tilde {q}$ on both the critical instability parameters and the burster phenomenon. That is to say, what are the lowest values of $\tilde {p}$ and $\tilde {q}$ for which the critical instability parameters and the burster phenomenon remain unaffected?

The frequency ratio $\omega = \sqrt {3}$ , for instance, can be approached by different ratios. If one considers the decimal number $7/4 = 1.75$ where $\tilde {p} = 7$ and $\tilde {q} = 4$ as an underestimated value of $\sqrt {3} \approx 1.73205$ with an absolute error of $0.01795$ , we are dealing with a periodic forcing system whose quasiperiodic character is not involved. This is not the case if one approaches the number $\sqrt {3}$ by the rational number $26/15 \approx 1.73333\ldots$ , where $\tilde {p} = 26$ and $\tilde {q} = 15$ , or by $97/56 \approx 1,7321428\ldots$ , where $\tilde {p} = 97$ and $\tilde {q} = 56$ . In these cases, the quasiperiodicity is strongly involved.

Figure 12 illustrates marginal stability curves when different approached values are assigned to the frequency ratio $\omega = \sqrt {3}$ when $\gamma _1 = 3.155$ . We consider especially the ratios $7/4, 26/15$ and $97/56$ . As one can notice, when $\omega = 7/4$ , only two smooth instability branches, corresponding to two regular modes, are observed, meaning that no bursting modes occur in this case. However, the left-hand stability branch becomes rough when $\omega = 26/15$ , indicating the existence of numerous instability modes, a hallmark of the bursting phenomenon. Such bursting becomes more pronounced and intense when $\omega = 97/56$ while the critical instability parameters remain almost unchanged. This feature means that the bursting phenomenon is a consequence of the non-periodic character in the forcing system that injects two different time scales into the stability problem. It is also interesting to point out that the critical instability parameters are weakly dependent on the nature of the approached value to the frequency ratio such as $\mathit{Ta}_c = 135.02$ at $q_c = 3.73$ at $\omega = 7/4$ , $\mathit{Ta}_c = 138.22$ at $q_c = 3.55$ at $\omega = 26/15$ against $\mathit{Ta}_c = 138.8$ at $q_c = 3.65$ at $\omega = 97/56$ . This means that approaching $\sqrt {3}$ by $362/209$ or $1351/780$ leads to a great increase in both the computational cost and the number of the stability modes constituting the rough branch in the stability curve in contrast to the critical instability parameters that remain unchanged.

By comparing these critical instability parameters with the benchmark values obtained by direct numerical simulations (performed recently by Riahi et al. Reference Riahi, Hayani Choujaa and Aniss2025) where the ratio $\sqrt {3}$ is fully taken into account, one can conclude that the Floquet predictions, based on rational approximations of $\sqrt {3}$ , exhibit rapid convergence toward the direct numerical simulation (DNS) results. The relative error between Floquet and DNS computations falls below 0.5 $\,\%$ for ratios higher than 7/4, confirming thus the accuracy and reliability of the rational-approximation method (see table 1 for more details).

Table 1.

Comparison between Floquet predictions obtained using different rational approximations of the irrational frequency ratio $\omega = \sqrt {3}$ and DNS results adopted by Riahi, Hayani Choujaa and Aniss (Reference Riahi, Hayani Choujaa and Aniss2025). The relative error is defined as $|\mathit{Ta}_c^{(\textit{Floquet})} - \mathit{Ta}_c^{(\textit{DNS})}|/\mathit{Ta}_c^{(\textit{DNS})}$ .

By limiting the frequency number $\sqrt {3}$ to the first number after the decimal point $\sqrt {3} \approx 1.7$ , which is equivalent to the decimal ratio number $17/10$ , only four instability modes are detected in the marginal stability curve computed at the same frequency $\gamma _1 = 3.155$ , as depicted in figure 13. However, by considering $\sqrt {3} \approx 1.73\,(173/100)$ we reproduce quite well the same marginal stability curve obtained at $\sqrt {3} \approx 97/56$ but the required calculation time becomes significantly larger.

Figure 12.

Marginal stability curves computed using Floquet theory at $\gamma _1 = 3.155$ when the frequency $\omega = \sqrt {3}$ is approached by different ratios: $\omega \approx 7/4$ , $\omega \approx 26/15$ and $\omega \approx 97/56$ . The vertical axis is the Taylor number $\mathit{Ta}$ while the horizontal one is the dimensionless axial wavenumber $q$ , ranging from 3 to 5.5.

Figure 13.

Marginal stability curves computed using Floquet theory at $\gamma _1 = 3.155$ when the frequency $\omega = \sqrt {3}$ is approached by different rational ratios: $\omega \approx 1.7$ and $\omega \approx 1.73$ . The vertical axis is the Taylor number $\mathit{Ta}$ while the horizontal one is the dimensionless axial wavenumber $q$ , ranging from 3 to 5.5.

4.4. Effect of the quasiperiodic frequency ratio $\omega$ on the stability diagrams of axisymmetric disturbances

4.4.1. The case $\omega \succ 1$ : $\omega = \sqrt {3}$

Figure 14.

Stability diagrams obtained by Floquet analysis in the case where $\omega = \sqrt {3}$ . Panels show (a, b) $\mathit{Ta}_c(\gamma _1)$ and (c, d) $q_c(\gamma _1)$ .

By considering axisymmetric disturbances, i.e. $m = 0$ , we illustrate in figure 14 a stability diagram showing the evolution of the critical Taylor number (figure 14 a,b) and wavenumber (figure 14 c,d) versus the inner cylinder oscillation frequency when the quasiperiodic frequency parameter $\omega = \sqrt {3}$ is approached by $97/56$ . Due to the strong presence of numerous instability modes, we plot these critical instability parameters in figure 14(a,c) in the frequency intervals $\gamma _1 \in [0, 3.5]$ and in figure 14(b,d) when $\gamma _1 \in [0, 10]$ . As one can see, these stability diagrams are composed of thirteen branches, resulting in twelve cusp points/discontinuities in the critical Taylor/wavenumber curve. This feature is due to the multiplicity of the flow states emerging under the effect of the irrational frequency number $\omega$ . The instability exchange between these states occurs via codimension-two bifurcation points, giving rise to twelve bicritical states in the computed marginal stability curves. This bicriticality sets in between either two harmonic, two subharmonic instability modes, or harmonic and subharmonic instability modes. For instance, the first bicritical state between two harmonic modes occurs at $\gamma _1 = 3.155$ , where a drop in $q_c$ from $5$ to $3.65$ is detected. This means that, at this frequency, two harmonic Taylor vortex flows with different critical wavenumbers bifurcate simultaneously from the laminar basic state. The flow reversal nature of these modes could be identified on the basis of time series of the local radial velocity of the disturbances whose change of sign indicates opposite polarity of the Taylor vortices. Recently, Riahi et al. (Reference Riahi, Hayani Choujaa and Aniss2025) have demonstrated that at $\gamma _1 = 10$ we are dealing with a regular non-reversing flow whereas at $\gamma _1 = 3.1$ a bursting reversing flow emerges. This means that the quasiperiodic nature of the frequency ratio alters the flow reversal of the system where non-reversing/reversing flows set in within the frequency interval $\gamma _1 \succ 3.2/ \gamma _1 \prec 3.2$ when $\omega = \sqrt {3}$ against $\gamma _1 \succ 1.3/\gamma _1 \prec 1.3$ when $\omega = 1$ (Riahi et al. Reference Riahi, Aniss and Ouazzani Touhami2019).

The second bicritical state occurs between two harmonic reversing flows at $\gamma _1 = 2.9$ , where a drop in $q_c$ from $4.5$ to $3.1$ is detected, while others are observed between harmonic and subharmonic modes at $\gamma _1 = 2.5, 2.3, 1.45, 1.27, 0.78$ . It is important to mention here that subharmonic instability modes do not appear when $\omega = 1$ and this is one of the important features induced by the presence of a two-frequency forcing.

In addition, one should highlight the flow stabilisation observed in the high-frequency limit where the critical Taylor number reproduces quite well the asymptotic behaviour obtained in the case $\omega = 1$ where: $\mathit{Ta}_c = 15.28\gamma ^{3/2}$ (Aouidef et al. Reference Aouidef, Normand, Stegner and Wesfreid1994; Riahi et al. Reference Riahi, Aniss and Ouazzani Touhami2019). Nevertheless, the low-frequency behaviour of the system is significantly altered by the presence of the quasiperiodicity $\omega$ . Indeed, the overall trend of the system is towards more destabilisation, especially when $\gamma _1 \prec 0.8$ , instead of a stabilisation in the same frequency range when $\omega = 1$ . This feature means that when $\omega = \sqrt {3}$ , the system is no longer potentially unstable at an intermediate-frequency value and therefore the frequencies $\gamma _{1p}$ and $\gamma _{2p}$ lose their physical relevance in the current context.

To understand this behaviour, one can recall the low-frequency behaviours of the system under different boundary conditions. Indeed, when the outer cylinder is slowly oscillating and the inner one is fixed, Riahi, Hayani Choujaa & Aniss (Reference Riahi, Hayani Choujaa and Aniss2024) have demonstrated using Floquet theory and two quasi-steady approaches that the system remains linearly stable, similarly to the steady rotation case of the cylinder. This result is in agreement with those obtained by Carmi & Tustaniwskyi (Reference Carmi and Tustaniwskyj1981). This means that the oscillation of the outer cylinder alone does not lead to any kind of instability, in contrast to the case where the only oscillating cylinder is the inner one, where the system loses stability at $\mathit{Ta}_c = 41.31$ and $q_c \to 3.12$ (Krueger et al. Reference Krueger, Gross and Di Prima1966; Riahi et al. Reference Riahi, Aniss, Ouazzani Touhami and Skali Lami2015) similarly to the steady rotation case. However, when both the inner and outer cylinders are slowly co-oscillating with the same frequency, the system becomes equivalent to a solid body rotation or to inviscid fluid flow where the shear strain between the cylinders vanishes and consequently $\mathit{Ta}_c$ tends towards infinity ( $\mathit{Ta}_c \to \infty$ ) (Aouidef et al. Reference Aouidef, Normand, Stegner and Wesfreid1994; Riahi et al. Reference Riahi, Aniss and Ouazzani Touhami2019). When the cylinders are slowly oscillating with two incommensurate frequencies, two irrational time scales are involved and a phase lag between the flows generated by each cylinder is created. A shear strain between the cylinders is then produced and the dynamics of the system becomes governed by the inner cylinder forcing, especially when the outer cylinder forcing does not lead to flow destabilisation, as demonstrated in Carmi & Tustaniwskyi (Reference Carmi and Tustaniwskyj1981), Riahi et al. (Reference Riahi, Hayani Choujaa and Aniss2024). In this framework, the low-frequency behaviour $\mathit{Ta}_c = 41.31$ and $q_c \to 3.12$ could be expected. This question will be discussed and studied in detail in § 5, where a quasi-steady analysis assuming low-frequency oscillations of the cylinders with an incommensurate frequency ratio is performed.

Figure 15.

Stability diagram obtained by Floquet theory over the frequency interval $\gamma _1 \in [7.6, 9]$ in the case where $\omega = 1/\sqrt {37}$ : (a) $\mathit{Ta}_c(\gamma _1)$ and (b) $q_c(\gamma _1)$ . Here, ‘s’ and ‘h’ denote subharmonic and harmonic modes, respectively.

4.4.2. The case $\omega \lt 1$ : $\omega = 1/ \sqrt {37}$

In figure 15, stability diagrams related to the case $\omega = 1/ \sqrt {37}$ when it is approached by $12/73$ are depicted for $\gamma _1 \in [7.6, 9]$ . As one can notice, the frequency range when the regular modes (non-reversing) set in shrinks, i.e. $\gamma _1 \succ 8.75$ in contrast to $\gamma _1 \succ 3.155$ when $\omega = \sqrt {3}$ . This means that, in comparison with the case $\omega = \sqrt {3}$ , the bursting frequency range becomes larger where we are dealing with a complex behaviour owing to the strong interplay between families of harmonic and subharmonic modes (reversing). Consequently, this stability diagram is composed of several instability branches featuring several codimension-two bifurcation points besides discontinuities in the critical wavenumber and cusp points in the Taylor number evolution, as shown in figure 16. With a decrease in the frequency $\gamma _1$ , the overall trend is towards destabilisation, similarly to the case $\omega = \sqrt {3}$ . In the low-frequency limit, one can envisage a constant critical Taylor number and wavenumber close to those calculated for a fixed outer cylinders, i.e. $\mathit{Ta}_c \to 41.31$ and $q_c \to 3.12$ (Krueger et al. Reference Krueger, Gross and Di Prima1966; Riahi et al. Reference Riahi, Aniss, Ouazzani Touhami and Skali Lami2015).

It is worth noting that a small frequency step is required, in this case, in order to capture all the eventually involved instability modes whose nature is very sensitive to the frequency number $\gamma _1$ .

The strong mode competition renders the experimental reproduction of the present results difficult. Pattern formation and abrupt changes in the critical wavenumber that may be obviously determined experimentally from the spatial structure of the Taylor cells could, however, be used as an efficient tool to distinguish between different primary instability modes arising when the frequency parameter $\gamma _1$ changes.

It is important to point out that the uncovered burster phenomenon in this study is a consequence of applying two incommensurate frequencies to the inner and outer drives of a Taylor–Couette system. These bursts are linked to resonance phenomena that could be induced or suppressed via an appropriate choice of input frequencies. Such behaviour has numerous real-world practical implications ranging from energy-efficient flow control and enhanced mixing to signal-based modulation in rotating systems. The dependence of bursting on the frequency forcing ratios suggests efficient strategies for regulating flow states in applications where hydrodynamic instabilities are undesired due to geometric constraints or high energy costs. By a proper tuning of the incommensurate frequency ratio, one can operate at regular, weakly and strongly bursting regimes and effectively the flow intensity or mixing characteristics are controlled without changing system geometry or fluid properties.

Figure 16.

Stability diagram obtained by Floquet theory over the frequency interval $\gamma _1 \in [0, 10]$ in the case where $\omega = 1/\sqrt {37}$ : (a) $\mathit{Ta}_c(\gamma _1)$ and (b) $q_c(\gamma _1)$ .

4.5. Effect of the quasiperiodic frequency ratio $\omega$ on the axisymmetric stability diagrams

In figure 17, we report stability diagrams computed at different values of the quasiperiodic frequency ratio in order to shed light on the effect of this latter on the onset of axisymmetric instabilities. This feature will be better elucidated by comparing the results related to the case $\omega \neq 1$ with those for $\omega = 1$ . It is important to mention here that two other values of $\omega$ are considered, namely $\omega = \sqrt {37}$ and $\omega = \sqrt {11}$ when they are approached respectively by $37/12$ and $63/19$ . As one can see, the quasiperiodic frequency parameter $\omega$ does not alter the stability of the system when the frequency of the cylinders is sufficiently high provided that $\omega \succ 1$ ( $\omega _2 \succ \omega _1$ ). In this case, we are dealing with an unstable oscillatory Stokes layer flow occurring especially in narrow regions near the inner cylinder. Therefore, the critical Taylor number in this frequency range is reminiscent of that obtained for the case of an oscillating inner cylinder with a fixed outer one and the cases of co- and counter-oscillating cylinders with $\omega = 1$ (Riahi et al. Reference Riahi, Aniss and Ouazzani Touhami2019, Reference Riahi, Hayani-Choujaa and Aniss2022): $\mathit{Ta}_c = 15.28 \gamma _1^{3/2}$ .

Figure 17.

Stability diagrams obtained by Floquet analysis for different values of the parameter $\omega$ .

The bifurcated solutions in this high-frequency limit correspond to regular non-reversing Taylor vortex flows that emerge when $\gamma _1 \succ 3.85$ for $\omega = \sqrt {1.1}$ , $\gamma _1 \succ 3.7$ for $\omega = \sqrt {1.2}$ , $\gamma _1 \succ 3.155$ for $\omega = \sqrt {3}$ , $\gamma _1 \succ 2.4$ for $\omega = \sqrt {11}$ , $\gamma _1 \succ 1.85$ for $\omega = \sqrt {37}$ and $\gamma _1 \succ 1.3$ for $\omega = 1$ . Therefore, one can conclude that the irrationality of the frequency ratio enhances/inhibits the reversing/non-reversing Taylor vortex flow. In comparison with the case $\omega = 1$ , the quasiperiodic frequency ratio in this case has a destabilising effect, which becomes more pronounced in the low-frequency limit. The parameter $\omega$ does not considerably alter the critical Taylor number in this frequency limit, where it is expected that the oscillating inner cylinder becomes the main destabilising drive of the system.

Particular behaviour is observed in the case $\omega \prec 1$ , i.e. $\omega _1 \succ \omega _2$ , especially at $\omega = 1/\sqrt {37}$ where non-reversing solutions emerge only when $\gamma _1 \succ 8.75$ . In this case, and in comparison with the case $\omega = 1$ , the quasiperiodic parameter ratio $\omega$ has a stabilising effect which becomes less pronounced by decreasing the frequency $\omega _1$ . A cross-over point is reached at $\gamma _1 \approx 1.5$ below which the effect of $\omega$ becomes destabilising. In the low-frequency limit, the critical Taylor number becomes close to the critical values calculated when $\omega \succ 1$ .

In previous studies dealing with time-quasiperiodic forcing, Boulal et al. (Reference Boulal, Aniss, Belhaq and Rand2007, Reference Boulal, Aniss, Belhaq and Azouani2008) and Boulal et al. (Reference Boulal, Riahi and Aniss2024), who were interested in studying the effect of quasiperiodic gravitational modulation on the stability of a heated fluid layer, have shown that the quasiperiodic frequency $\omega$ alters the stability characteristics of the system if it belongs to the interval $(( {1}/{\sqrt {37}}), \sqrt {37})$ . Similar behaviour has been observed by Yagoubi & Aniss (Reference Yagoubi and Aniss2017), who studied the effect of vertical quasiperiodic vibrations on the Faraday instability. When dealing with time-quasiperiodically modulated Taylor–Couette flow, this feature is conserved only when $\gamma _1 \succcurlyeq 1.5$ .

4.6. Effect of the quasiperiodic frequency ratio $\omega$ on the primary bifurcation taking into account non-axisymmetric disturbances

In this subsection, we shed light on the ability of the quasiperiodic frequency ratio $\omega$ in changing the axisymmetric nature of the primary bifurcation with respect to the case $\omega = 1$ studied experimentally by Aouidef et al. (Reference Aouidef, Normand, Stegner and Wesfreid1994). This purpose can be reached by a numerical resolution of the stability equations under small and continuous variations of the wavenumber $k$ . The resulting marginal stability curves, showing evolution of the Taylor number versus the axial wavenumber, allow the identification of the most unstable mode. The latter is three-dimensional and non-axisymmetric if the critical Taylor number calculated at $k = 0$ is higher than that at $k \neq 0$ . The critical value of the wavenumber $k$ , denoted $k_c$ , corresponds to that of the lowest value of the critical Taylor number. In such cases, the critical azimuthal wavenumber is calculated from the relation $k_c = m_c \varepsilon ^{1/2}$ for a given value of $\varepsilon \ll 1$ , where only positive integer values of $m_c$ are physically significant. This procedure requires a wide enough range of real values of $k$ to accurately identify the nature of the most unstable mode and determine its critical wavenumber $m_c$ . However, one can achieve the same objective in a simpler manner by choosing a definite value of $\varepsilon$ and numerically resolving the stability equations when only integer values are assigned to the azimuthal wavenumber $m$ . In the present paper, $\varepsilon = 0.14$ is adopted in our numerical calculations since it is a reasonable small-gap value in typical experimental geometries. We recall here that Feugaing et al. (Reference Feugaing, Crumeyrolle, Yang and Mutabazi2014) carried out experiments for $\varepsilon = 0.139$ in the case where the outer cylinder is fixed while the inner one is rotating with an angular velocity modulated in time about a non-zero mean value. At sufficiently small modulation amplitudes, their experimental instability parameters are in excellent agreement with the classical ones calculated using the narrow-gap approximation (Chandrasekhar Reference Chandrasekhar1961). Experimental results obtained by Aouidef et al. (Reference Aouidef, Normand, Stegner and Wesfreid1994) with $\varepsilon = 0.11$ reproduce quite well their numerical findings in the co-oscillating case with $\omega = 1$ using the narrow-gap approximation. Due to the rich and complex dynamics observed at $m = 0$ and owing to the wide frequency range considered in this study, our attention is restricted only to some values of the frequency $\gamma _1$ .

For $\omega = \sqrt {3}$ , we found that the most unstable mode is no longer axisymmetric, especially at the frequencies $\gamma _1 = 1.5$ , $\gamma _1 = 2$ and $\gamma _1 = 3.5$ where $m_c = 1$ and $\gamma _1 = 2.5$ and at $\gamma _1 = 3$ where $m_c = 2$ . These results are presented in table 2. It is also important to mention that the quasiperiodic parameter alters also the type of the most unstable mode where non-axisymmetric subharmonic responses set in at some frequencies, e.g. $\gamma _1 = 2$ , where the axisymmetric modes are harmonic.

Table 2.

Numerical values of the critical instability parameters, i.e. $(\mathit{Ta}_c, q_c, m_c)$ , at different values of the frequency $\gamma _1$ when the quasiperiodic frequency parameter $\omega = \sqrt 3$ .

The same dynamics is observed when $\omega = 1/\sqrt {37}$ , as shown in table 3. One can mention in this case that non-axisymmetric instability modes are more suitable than axisymmetric modes.

Table 3.

Numerical values of the critical instability parameters, i.e. $(\mathit{Ta}_c, q_c, m_c)$ , at different values of the frequency $\gamma _1$ when the quasiperiodic frequency parameter $\omega = 1/\sqrt {37}$ .