2.1. Studied configuration and governing equations

We consider an axisymmetric acoustic streaming jet flowing in a closed cylindrical cavity (figure 1). The problem is defined using a cylindrical coordinate system in which $\boldsymbol{e_x}$ , $\boldsymbol{e_r}$ and $\boldsymbol{e_\theta }$ refer to the axial, radial and azimuthal directions, respectively. The spatial coordinates in the corresponding directions are $x$ , $r$ and $\theta$ ; the origin is centred on the upstream wall (figure 1) so that the cavity axis is $r=0$ . Finally, we shall define the hydrodynamics and acoustics problems using the same dimensionless parameters and variables as in our previous work (Vincent et al. Reference Vincent, Henry, Kumar, Botton, Pothérat and Miralles2025). In particular, all lengths are normalised by the dimensional source diameter $D_s$ , and we shall use $\rho D_s^3$ as mass scale (where $\rho$ is the dimensional fluid density) and $D_s^2 / \nu$ as time scale (where $\nu$ stands for the dimensional kinematic viscosity). From now on, we shall only rely on non-dimensional parameters based on these scales.

Figure 1. Sketch of the investigated set-up in the $( x, r )$ plane, where $x$ and $r$ are the axial and radial coordinates, respectively. The circular transducer of unit diameter (grey rectangle) at $x=0$ emits a beam of linear sound waves (green shaded area) within an elongated cylindrical cavity of length $L$ and diameter $D$ . As the acoustic pressure waves travel across the cavity, their amplitude decays at a rate $N/L$ , where $N$ is the ratio between the cavity length and the acoustic pressure attenuation length. Attenuation of these waves forces a mean flow of velocity $\boldsymbol{u}$ . All the boundaries are no-slip walls and completely absorb the acoustic waves.

The cavity, of diameter $D$ and length $L$ , is entirely filled with a Newtonian fluid. The cavity is fitted with a plane circular transducer of unit diameter at $x=0$ . The transducer, whose axis is aligned with $\boldsymbol{e_x}$ , radiates an axisymmetric beam of linear acoustic pressure waves. As these waves travel towards the opposite cavity end, their amplitude decays for two reasons: diffraction, which also causes the beam radius to increase at a rate $S$ (figure 1), and sound attenuation. The latter gives rise to an axisymmetric force driving a streaming flow of velocity $\boldsymbol{u}$ . The injected momentum is scaled by $M$ , the maximum dimensionless momentum flux to be acquired by the jet if $L \rightarrow \infty$ . The exponential decay rate of the acoustic pressure wave amplitudes is defined by the sound attenuation coefficient $N/L$ , where $N$ is the ratio between the cavity length and the attenuation length of the acoustic pressure wave amplitude. All cavity walls absorb incident sound waves thus preventing the development of standing sound waves. The studied configuration is thus defined by five non-dimensional parameters that are defined in terms of dimensional parameters as follows (Vincent et al. Reference Vincent, Henry, Kumar, Botton, Pothérat and Miralles2025):

(2.1) \begin{equation} L = \frac {L_c}{D_s} , \quad D = \frac {D_c}{D_s} , \quad N = \alpha L_c , \quad S = 1.22 \frac {\lambda }{D_s} , \quad M = \frac {P_{\textit{ac}}}{\rho c \nu ^2} , \end{equation}

where $L_c$ and $D_c$ are the dimensional cavity length and diameter, respectively. The dimensional acoustic pressure attenuation coefficient is denoted by $\alpha$ and the dimensional speed of sound by $c$ . The dimensional acoustic power radiated by the transducer is denoted by $P_{\textit{ac}}$ and $\lambda$ is the dimensional wavelength.

Whilst the acoustic field is characterised by small length scales and fast time scales, the streaming flow associated with them is significantly slower than the speed of sound and involves length scales that greatly exceed the acoustic wavelength (Lighthill Reference Lighthill1978; Orosco & Friend Reference Orosco and Friend2022). These scale differences allow us to model the mean motion as that of an incompressible flow due to an external momentum source (Moudjed et al. Reference Moudjed, Botton, Henry, Ben Hadid and Garandet2014a ):

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

(2.3) \begin{align} \frac { \partial \boldsymbol{u} }{\partial t} + \left ( \boldsymbol{u} \boldsymbol{\cdot } \boldsymbol{\nabla } \right ) \boldsymbol{u} &= - \boldsymbol{\nabla } p + {\nabla }^2 \boldsymbol{u} + \boldsymbol{F_{\textit{ac}}} , \end{align}

(2.4) \begin{align} \boldsymbol{F_{\textit{ac}}} &= {\textit{Gr}}_{\textit{ac}} \boldsymbol{ \widetilde {I} }. \end{align}

Here, $p$ and $\boldsymbol{u}$ respectively refer to the operating pressure and velocity fields averaged over one acoustic period of time. The operating pressure, $p$ , accounts for the weight of the fluid averaged over one acoustic period. The volume momentum source for the mean flow, $\boldsymbol{F_{\textit{ac}}}$ , couples the acoustics and the hydrodynamics problems. The spatial structure of $\boldsymbol{F_{\textit{ac}}}$ is set by the normalised acoustic intensity $\boldsymbol{\widetilde {I}}$ , which is bounded by unity for an unattenuated sound field. The magnitude of the forcing is defined by the Grashof number:

(2.5) \begin{equation} {\textit{Gr}}_{\textit{ac}} = \frac {32 \alpha P_{\textit{ac}} D_s}{\pi \rho c \nu ^2} = \frac {32}{\pi } \frac {\textit{N M}}{L}, \end{equation}

comparing the magnitude of the acoustic forcing $\boldsymbol{F_{\textit{ac}}}$ with the magnitude of the viscous forces at the scale of the transducer (El Ghani et al. Reference El Ghani, Miralles, Botton, Henry, Ben Hadid, Ter-Ovanessian and Marcelin2021).

The cavity walls impose the no-slip boundary condition:

(2.6) \begin{equation} \boldsymbol{u} = \boldsymbol{0}, \end{equation}

at $x \in \left \{ 0, L \right \}$ and $r = D/2$ . Additionally, the flow is axisymmetric and swirl-free, so that for all $x$ and $r$ :

(2.7) \begin{equation} u_{\theta } = \frac {\partial u_x}{\partial \theta } = \frac {\partial u_r}{\partial \theta } = \frac {\partial p}{\partial \theta } = 0, \end{equation}

where $u_x$ , $u_r$ and $u_{\theta }$ refer to the axial, radial and azimuthal velocity components, respectively. Axisymmetry of both $\boldsymbol{u}$ and $p$ further imposes

(2.8) \begin{equation} u_r = \frac {\partial u_x}{\partial r} = \frac {\partial p}{\partial r} = 0 \end{equation}

on the axis ( $r=0$ ).

To model the forcing field $\boldsymbol{F_{\textit{ac}}}$ , we assume the following.

1. (i) The length and time scales associated with the acoustic problem are significantly smaller than those of the streaming flow, as verified previously in experiments (Kamakura et al. Reference Kamakura, Sudo, Matsuda and Kumamoto1996; Mitome Reference Mitome1998; Moudjed et al. Reference Moudjed, Botton, Henry, Ben Hadid and Garandet2014a ). As a result, both the streaming flow velocity levels and gradients do not alter the propagation of sound (Pierce Reference Pierce1990). Therefore, $\boldsymbol{F_{\textit{ac}}}$ is time-independent.

2. (ii) The acoustic waves are irrotational. This common assumption in acoustics is known to hold in the bulk of the fluid (Coulouvrat Reference Coulouvrat1992; Kinsler et al. Reference Kinsler, Frey, Coppens and Sanders2000; Blackstock Reference Blackstock2000); however, sound waves are no longer irrotational in a thin viscous layer called the acoustic boundary layer (ABL). The ABL develops near solid boundaries, including at the transducer surface, and forces the acoustic velocity fluctuation to match the boundary velocity. For the experiment of Moudjed et al. (Reference Moudjed, Botton, Henry, Ben Hadid and Garandet2014a ), on which we shall base the structure of $\boldsymbol{F_{\textit{ac}}}$ , the ABL thickness is $\delta \approx 1.3 \times 10^{-5}$ . As $\delta$ is several orders of magnitude smaller than both our cavity sizes and flow length scales, the supplementary attenuation and possible streaming flow within the ABL can be safely neglected (Dridi et al. Reference Dridi, Henry and Hadid2010).

3. (iii) Both $L$ and the sound attenuation distance $L / N$ are significantly smaller than the acoustic shock formation distance $L_s$ . Nonlinear acoustic effects can thus be discarded: the acoustic waves are linear and the superposition principle applies. Accounting for acoustic nonlinearities would require two additional dimensionless parameters: a frequency parameter, as in Moudjed et al. (Reference Moudjed, Botton, Henry, Ben Hadid and Garandet2014a ), and the acoustic nonlinearity parameter (Hamilton & Blackstock Reference Hamilton and Blackstock2024). The reduced number of fluid-dependent parameters allowed by the linear acoustics framework thus makes our stability study as generic as possible. In practice, our forcing model is appropriate for $\min ( L, L/N ) \lt L_s$ .

4. (iv) The waves are time-harmonic monochromatic signals.

5. (v) Sound attenuation is due to the viscosity of the medium only. Coulouvrat (Reference Coulouvrat1992) showed that thermally-induced attenuation is of the order of $ ( \gamma - 1 )$ , where $\gamma$ stands for the adiabatic index. For common liquids such as water, $\gamma \approx 1$ ; hence, thermally-induced attenuation can be dropped (Riaud et al. Reference Riaud, Baudoin, Bou Matar, Thomas and Brunet2017) (this approximation however does not hold for liquid metals). We also discard any sound attenuation caused by relaxation.

6. (vi) The sound attenuation length $L/N$ is large compared with the wavelength (i.e. $NS / (1.22 L ) \ll 1$ ).

Under these assumptions, we adapted in a previous work the classical equations for a plane baffled piston in a semi-infinite domain to account for viscous attenuation of sound (Vincent et al. Reference Vincent, Henry, Kumar, Botton, Pothérat and Miralles2025). We showed that $\boldsymbol{ \widetilde {I} }$ can be obtained from a set of integrals defining the acoustic pressure $\widetilde {p}_{\textit{ac}}$ and velocity $\boldsymbol{ \widetilde {u}_{\textit{ac}} }$ at any point $ ( x, r )$ in the fluid domain:

(2.9a) \begin{align} \boldsymbol{ \widetilde {I} } &= \textrm{Re} \left \{ \widetilde {p}_{\textit{ac}} \, \left ( \boldsymbol{ \widetilde {u}_{\textit{ac}} }^* \right ) \right \} , \end{align}

(2.9b) \begin{align} \widetilde {p}_{\textit{ac}} &= \mathrm{i} \frac { k^2 S }{8 \times 1.22 \pi ^2} \int _{0}^{1/2} \int _{0}^{2 \pi } \frac { e^{- \mathrm{i} k \Vert \boldsymbol{d} \Vert } }{ \Vert \boldsymbol{d} \Vert } r'\, \mathrm{d}r'\, \mathrm{d}\theta ' , \end{align}

(2.9c) \begin{align} \boldsymbol{ \widetilde {u}_{\textit{ac}} } &= \frac {1}{4 \pi } \int _{0}^{1/2} \int _{0}^{2\pi } \left ( 1 + \mathrm{i} k \Vert \boldsymbol{d} \Vert \right ) \frac { e^{- \mathrm{i} k \Vert \boldsymbol{d} \Vert } }{ \Vert \boldsymbol{d} \Vert ^2} \frac { \boldsymbol{d} }{ \Vert \boldsymbol{d} \Vert } r' \,\mathrm{d}r'\, \mathrm{d}\theta ' , \end{align}

(2.9d) \begin{align} k &= \frac {2.44 \pi }{S} - \mathrm{i} \frac {N}{L} , \end{align}

where $^{*}$ refers to the complex conjugate and $\textrm{Re}$ to the real part. In (2.9), $\boldsymbol{d} = x \, \boldsymbol{e_x} + [ r - r' \cos ( \theta ' ) ] \boldsymbol{e_r} - r' \sin ( \theta ' ) \boldsymbol{e_{\theta }}$ is the distance between any point source on the transducer surface (located by the set of radial and azimuthal coordinates $r'$ and $\theta '$ ) and any point at $ ( x, r )$ . Finally, $\boldsymbol{\widetilde {I}}$ is mainly shaped by the complex wave number $k$ : the real part of $k$ sets the beam diffraction, whilst the imaginary part of $k$ involves the sound attenuation coefficient $N/L$ and thus adds attenuation to $\boldsymbol{\widetilde {I}}$ .

Equations (2.9a )–(2.9d ) define a beam-shaped axisymmetric acoustic field. As (i) the beam expands radially at a rate $S$ along $x$ , and (ii) the beam radius matches the transducer radius (0.5 in dimensionless coordinates) at $x=0$ , an approximate description of the beam radius $R_{\textit{beam}}$ is

(2.10) \begin{equation} R_{\textit{beam}} = 0.5 + S x , \end{equation}

from which we define the confinement ratio $C_R$ :

(2.11) \begin{equation} C_R = \frac { D/2 }{\left . R_{\textit{beam}} \right |_{x=L} } . \end{equation}

From now, we shall use $C_R$ as the main indicator of the expected radial confinement of the flow; small (respectively large) $C_R$ values correspond to cases of high (respectively weak) confinement.

To investigate the effect of confinement on the flow stability, we shall consider six different cavity sizes by varying either $C_R$ or $N$ (table 1). We take $C_R \in [ 1, 6 ]$ to be consistent with the experiments of Kamakura et al. (Reference Kamakura, Sudo, Matsuda and Kumamoto1996), Mitome (Reference Mitome1998) and Moudjed et al. (Reference Moudjed, Botton, Henry, Ben Hadid and Garandet2014a ) for which $C_R=2.9$ , $8.5$ and $4.1$ , respectively. As these experiments were made for $N \lt 1$ , we set $N=0.25$ for five cases. Nevertheless, from a process standpoint, it is preferable to have $N \geqslant 1$ to ensure that most of the acoustic kinetic energy radiated by the transducer is used to create flow kinetic energy (Vincent et al. Reference Vincent, Henry, Kumar, Botton, Pothérat and Miralles2025). We shall thus investigate a last case for which $ ( N, C_R ) = ( 1, 6 )$ (table 1). For all six cases, the parameters $S = 0.03$ and $N/L = 0.003$ are kept constant; these values are based on the 2 MHz water experiments of Moudjed et al. (Reference Moudjed, Botton, Henry, Ben Hadid and Garandet2014a , Reference Moudjed, Botton, Henry, Millet, Garandet and Ben Hadid2015). Therefore, the structure of $\boldsymbol{F_{\textit{ac}}}$ in our work is consistent with the forcing fields of former experiments. These experiments covered a wide range of ${\textit{Gr}}_{\textit{ac}}$ ; it typically ranged from $5.6 \times 10^2$ (Mitome Reference Mitome1998) to $1.2 \times 10^5$ (Moudjed et al. Reference Moudjed, Botton, Henry, Ben Hadid and Garandet2014a ). In the present work, we shall adjust ${\textit{Gr}}_{\textit{ac}}$ on a case-by-case basis to identify the instability onset.

Table 1. Values of the parameters defining each set-up. The parameter $C_R$ is defined as the ratio between the cavity radius $D/2$ and the approximate beam radius evaluated at $x=L$ using (2.10). The values of $S$ and of the acoustic pressure attenuation coefficient $N/L$ are chosen to match the 2 MHz water experiments of Moudjed et al. (Reference Moudjed, Botton, Henry, Ben Hadid and Garandet2014a ).

2.2. Linear stability analysis

Although the base flow is 2-D–2-C (two-dimensional, two-component), it may nevertheless be destabilised by non-axisymmetric 3-D–3-C perturbations under certain forcing and geometrical conditions. We shall therefore decompose the flow fields into a steady axisymmetric part, corresponding to the base flow, and an infinitesimal time-dependent perturbation:

(2.12) \begin{align} \boldsymbol{u} &= \boldsymbol{U}\left ( x, r \right ) + \boldsymbol{u'}\left ( x, r, \theta , t \right )\! , \end{align}

(2.13) \begin{align} p &= P\left ( x, r \right ) + p'\left ( x, r, \theta , t \right ) \! . \end{align}

The perturbations $ ( \boldsymbol{u'}, p' )$ are purely hydrodynamic i.e. they do not correspond to any acoustic signals as feedback effects of the flow on the acoustic fields are discarded. The linear stability equations are then obtained by inserting the decompositions (2.12) and (2.13) into (2.3) and (2.2) describing the mean flow, and by linearising the resulting equations around the base flow. Since $ ( \boldsymbol{U}, P )$ satisfy (2.3) and (2.2), the linear equations for $ ( \boldsymbol{u'}, p' )$ are

(2.14) \begin{align} \frac { \partial \boldsymbol{u'} }{ \partial t } + \left ( \boldsymbol{U} \boldsymbol{\cdot } \boldsymbol{\nabla } \right ) \boldsymbol{u'} + \left ( \boldsymbol{u'} \boldsymbol{\cdot }\boldsymbol{\nabla } \right ) \boldsymbol{U} &= - \boldsymbol{\nabla } p' + {\nabla }^2 \boldsymbol{u'} , \end{align}

(2.15) \begin{align} \boldsymbol{\nabla } \boldsymbol{\cdot } \boldsymbol{u'} &= 0 . \end{align}

As the base flow is invariant along $\boldsymbol{e_{\theta }}$ , a generic three-dimensional perturbation $\boldsymbol{q'}$ is expressed using the Fourier series:

(2.16) \begin{equation} \boldsymbol{q'} \left ( x, r, \theta , t \right ) = \sum _{m=-\infty }^{+\infty } \widehat { \boldsymbol{q'} }_m \left ( x, r, t \right )e^{\textrm {i} m \theta } , \end{equation}

where the integer $m$ is the azimuthal wavenumber. Since (2.14) and (2.15) are linear, the LSA is thus decoupled for each $m$ . The linearity of (2.14) and (2.15) decouples the Fourier modes, so that the LSA equations can be solved individualy for each $m$ . Finally, the perturbations satisfy the same boundaries conditions as the base flow, except at $r=0$ , where $m$ -dependent kinematic conditions are imposed to allow for three-dimensional perturbations to develop whilst avoiding the geometrical singularity (Blackburn & Sherwin Reference Blackburn and Sherwin2004).

The linear stability problem is addressed using a timestepper approach. From (2.14) and (2.15), the perturbation at a time $t + \tau$ results from the action of an evolution operator $\mathcal{A}(\tau )$ applied to the same perturbation at a previous time $t$ (Barkley, Blackburn & Sherwin Reference Barkley, Blackburn and Sherwin2008):

(2.17) \begin{equation} \boldsymbol{q'}\left ( t + \tau \right ) = \mathcal{A}\left ( \tau \right ) \boldsymbol{q'}\left ( t \right ) . \end{equation}

The complex eigenvalue $\zeta _k$ of $\mathcal{A} ( \tau )$ , associated with the eigenvector $\boldsymbol{q'}_k$ , is linked to the eigenvalue $ ( \sigma _k + \mathrm{i} \omega _k )$ of the linear problem (2.14)–(2.15) through an exponential mapping:

(2.18) \begin{equation} \zeta _k = e^{\left ( \sigma _k + \mathrm{i} \omega _k \right ) \tau } , \end{equation}

from which the growth rate $\sigma _k$ and frequency $\omega _k$ can be retrieved:

(2.19a,b) \begin{align} \sigma _k = \frac { \textrm {ln}\left ( \vert \zeta _k \vert \right ) }{\tau }, \quad \omega _k = \frac {\varphi _k}{\tau } , \end{align}

where $\varphi _k$ is the complex argument of $\zeta _k$ . A given base flow $ ( \boldsymbol{U}, P )$ is said to be unstable if, for any $m$ , there is at least one perturbation $ ( \boldsymbol{u'}, p' )$ for which $\sigma \gt 0$ (or, equivalently, $\vert \zeta \vert \gt 1$ ). Conversely, the base flow is stable if the growth rates associated with all the eigenvalues of the linear problem are strictly negative. The minimum value of ${\textit{Gr}}_{\textit{ac}}$ for which $\sigma = 0$ is reached is called the critical acoustic Grashof number ${\textit{Gr}}_{\textit{ac}}^c$ and indicates the onset of an instability.

3.1. Numerical methods

Our work involves four different numerical calculations. First, we calculate $\boldsymbol{F_{\textit{ac}}}$ at the grid points in the fluid domain. Second, we compute the steady axisymmetric 2-D–2-C streaming flow driven by $\boldsymbol{F_{\textit{ac}}}$ by time-stepping (2.2)–(2.3) until a steady-state is reached. Third, we perform the LSA of the axisymmetric 2-D–2-C flow for each case listed in table 1. The leading eigenmode and its eigenvalue are computed by solving the eigenvalue problem for $0 \leqslant m \leqslant 10$ . This allows us to determine ${\textit{Gr}}_{\textit{ac}}^c$ and identify the shape of the unstable perturbation for each case. Fourth, nonlinear 3-D–3-C unsteady simulations are run for ${\textit{Gr}}_{\textit{ac}} \gt {\textit{Gr}}_{\textit{ac}}^c$ for two purposes: (i) to validate the LSA and (ii) to determine the bifurcation nature (sub- or supercritical). The latter is inferred by modelling the growth and saturation of a perturbation near ${\textit{Gr}}_{\textit{ac}}^c$ with a Stuart–Landau equation (Landau & Lifschitz Reference Landau and Lifschitz1987), the parameters of which being evaluated from time series of the flow variables (Sheard, Thompson & Hourigan Reference Sheard, Thompson and Hourigan2004).

The forcing field is obtained by numerically evaluating the integrals in (2.9) at the grid points. The transducer is discretised by an homogeneous distribution of $N_s$ point sources in both the radial and the azimuthal directions. The total number of point sources used to estimate $\boldsymbol{ \widetilde {I} }$ (hence, $\boldsymbol{F_{\textit{ac}}}$ ) at any grid point in the fluid domain is thus $N_s^2$ .

Both the steady axisymmetric 2-D–2-C and the unsteady 3-D–3-C flows are computed using the open-source spectral-element code Semtex (Blackburn et al. Reference Blackburn, Lee, Albrecht and Singh2019). The $ (x, r )$ plane is split into non-overlapping elements, inside which the flow variables are expanded using a tensor product of Lagrange polynomials defined on the Gauss–Lobatto–Legendre points (Blackburn & Sherwin Reference Blackburn and Sherwin2004). Convergence of the solution is then achieved by increasing the polynomial degree $N_p$ of the expansion basis. For the 3-D–3-C unsteady simulations, the azimuthal direction is discretised using a Fourier spectral method (Blackburn & Sherwin Reference Blackburn and Sherwin2004).

The $ (x, r )$ plane is discretised using a grid of quadrilateral elements (figure 2). For instance, the mesh used for the $ ( N, C_R )= ( 0.25, 4 )$ comprises 24 and 9 elements in the axial and radial directions, respectively. For each case, the grid is refined near the upstream and downstream walls, so that the ratio of the lengths of two neighbouring elements does not exceed 1.35. In the radial direction, the element height is increased by a factor of 1.3 as the distance from the axis is increased. The grid density in the radial direction is set such that there is at least one element in $0 \leqslant r \leqslant 0.5$ . We also adapted the total number of elements to keep a similar spatial resolution of $\boldsymbol{F_{\textit{ac}}}$ in the very-near acoustic field ( $x \leqslant 1$ ) for each case (table 2).

Figure 2. Typical mesh used to discretise the fluid domain. The circular transducer is placed at $x=0$ and its axis is aligned with the $r=0$ line. The beam, of approximate radius $R_{\textit{beam}}$ (2.10), is shown in green. The displayed mesh corresponds to the $(N, C_R) = ( 0.25, 4 )$ case.

Table 2. Characteristics of the meshes used for each case. See table 1 for a complete list of the computational parameters defining each case.

Equations (2.2)–(2.3) are discretised in time using the operator splitting scheme of Karniadakis, Israeli & Orszag (Reference Karniadakis, Israeli and Orszag1991). The scheme integrates explicitly the inertia terms of (2.3) and implicitly its viscous terms. For the unsteady 3-D–3-C simulations, the temporal integration scheme is used in its third-order formulation. For the 2-D–2-C base flow calculations, as only steady solutions of (2.2)–(2.3) are sought, the scheme is used in its first-order version. A steady state is then considered to be reached when the time series of the flow variables, recorded at different locations within the fluid domain, do not display variations with time up to a precision of seven significant figures.

We carry out the LSA using the open-source code DOG (Direct Optimal Growth, Barkley et al. (Reference Barkley, Blackburn and Sherwin2008)), which is based on a time-stepper method and a spectral element spatial discretisation. The code marches forward in time the linear stability equations (2.14)–(2.15) using the second-order formulation of the operator splitting scheme of Karniadakis et al. (Reference Karniadakis, Israeli and Orszag1991). The eigenvalues and associated eigenvectors of the stability problem are then recovered by applying the iterative algorithm of Barkley et al. (Reference Barkley, Blackburn and Sherwin2008).

Finally, for a given set of values of $N$ , $C_R$ and ${\textit{Gr}}_{\textit{ac}}$ , the initial condition for the 2-D–2-C flow calculations is either set to $\boldsymbol{u}=0$ or to a previous solution obtained for the same $ ( N, C_R )$ but for a different ${\textit{Gr}}_{\textit{ac}}$ to accelerate the convergence towards a steady state. The unsteady 3-D–3-C simulations of the ${\textit{Gr}}_{\textit{ac}} \gt {\textit{Gr}}_{\textit{ac}}^c$ regimes are initialised with the steady axisymmetric base flow duplicated along $\boldsymbol{e}_{\theta }$ , plus a Gaussian white noise. The noise is evenly distributed over all the azimuthal wavenumbers, thus ensuring that the initial disturbance is not biased towards a preferred shape. The standard deviation of the white noise was set to 1 for the 3-D–3-C unsteady simulations run for $C_R = 6$ . As we reduced $C_R$ , such standard deviation yielded very large initial disturbances, making the early exponential growth of the perturbation impossible to observe. We thus reduced the noise standard deviation to 0.1 for $C_R \in \left \{ 3, 4 \right \}$ , and to 0.01 for $C_R \in \left \{ 1, 2 \right \}$ .

3.2. Grid sensitivity analysis

We assessed the suitability of the grids discretising the transducer and the fluid domain. For the transducer discretisation, we previously derived an analytical expression for the on-axis intensity $\widetilde {I}_x(x,r=0)$ , which we compared with the numerical evaluation of (2.9) using the methodology outlined in § 3.1 (Vincent et al. Reference Vincent, Henry, Kumar, Botton, Pothérat and Miralles2025). We showed that discretising the transducer with $N_s=300$ grid points in both the radial and azimuthal directions ensures a local error of less than 1 % between the analytical and numerical intensity profiles on the axis (represented by the black curve and red stars, respectively, in figure 3). All the results presented in this work were thus obtained with $N_s=300$ .

Figure 3. On-axis normalised acoustic intensity profile $\widetilde{I}_x (x, r=0)$ for $( N, C_R )=( 0.25, 4 )$ . The analytical profile, proposed by Vincent et al. (Reference Vincent, Henry, Kumar, Botton, Pothérat and Miralles2025), is shown in black. The transition between the near and far acoustic fields occurs at the last intensity peak located at approximately the Fresnel distance $L_F \approx 1.22 / (4 S )$ (purple dashed vertical line). The values obtained by numerically evaluating (2.9) at the collocation points of the mesh (with an expansion basis of degree $N_p=8$ in each element) are represented by red stars. The numerical points are obtained by discretising the transducer with $N_s = 300$ point sources in both the radial and the azimuthal directions.

As shown in figure 3, the forcing field features strong gradients in the near acoustic field. These sharp intensity variations occur over very short length scales and would require an extremely fine grid to be accurately captured. As these length scales are significantly smaller than those of the studied streaming flow (Vincent et al. Reference Vincent, Henry, Kumar, Botton, Pothérat and Miralles2025), the very near acoustic field (i.e. $x \leqslant 1$ ) is undersampled to keep the flow calculations tractable. We nevertheless performed $p$ -refinement to determine the minimum $N_p$ above which the variations of the leading eigenvalue remained below 0.5 %. An example of grid sensitivity analysis for the $ ( N, C_R ) = ( 0.25, 4 )$ case at ${\textit{Gr}}_{\textit{ac}} = 7000$ (the largest ${\textit{Gr}}_{\textit{ac}}$ for that set-up) is shown in table 3. The leading perturbation is an unstable non-oscillatory $m=2$ mode. We chose $N_p=8$ , since variations of $\sigma$ with respect to the most refined grid were at most approximately $ 0.1$ % for larger $N_p$ . We then repeated this $p$ -refinement study for all $ ( N, C_R )$ to ensure that our stability results are mesh-independent. We also checked that the acoustic field is sufficiently well resolved to inject the right amount of momentum into the fluid. Relying on momentum balances, the error made on the injected momentum is at most 3.9 % despite the very near acoustic field being undersampled (see the Appendix).

Table 3. Evolution of the leading mode growth rate $\sigma _{max}$ with the polynomial degree $N_p$ of the expansion basis. The growth rates are obtained for $( N, C_R )=( 0.25, 4 )$ with ${\textit{Gr}}_{\textit{ac}} = 7000$ and $m=2$ . The error is relative to the value of $\sigma _{max}$ obtained for $N_p = 14$ .

For the 3-D–3-C unsteady simulations, $\boldsymbol{u}$ and $p$ were expanded along $\boldsymbol{e_{\theta }}$ using $N_F$ Fourier modes. We assessed the effect of $N_F$ by computing the flow kinetic energy over the entire cavity after the perturbation had saturated. For $ ( N, C_R ) = ( 0.25, 1 )$ with ${\textit{Gr}}_{\textit{ac}} = 15\,500$ , for instance, the relative difference between the kinetic energies computed with either $N_F = 32$ or 64 was approximately $ 0.03$ %. Even smaller differences were obtained for the ${\textit{Gr}}_{\textit{ac}} \gt {\textit{Gr}}_{\textit{ac}}^c$ simulations of the $C_R \in \left \{ 4, 6 \right \}$ cases. The 3-D–3-C unsteady simulations were thus performed with $N_F = 32$ for all cases.

Finally, we adapted the time step $\Delta t$ to each case so that the Courant–Friedrichs–Lewy (CFL) number remained below unity to ensure the numerical stability. For instance, the $ ( N, C_R ) = ( 0.25, 4 )$ case at ${\textit{Gr}}_{\textit{ac}} = 7000$ and with $N_p = 8$ was time-stepped with $\Delta t = 10^{-4}$ , yielding a maximum CFL number of 0.476.

4.1. Impinging jet structure of the streaming flow in a weakly confined setting

Before studying its stability, we shall first discuss the steady two-dimensional and axisymmetric base flow. Figure 4 shows a typical base flow considered for LSA, with the $ ( N, C_R )= ( 0.25, 6 )$ case taken as an illustrative example. The main flow features are identified by the patterns of streamlines together with critical point theory (Chong, Perry & Cantwell Reference Chong, Perry and Cantwell1990). The critical point analysis is commonly used to identify structures within complex 3-D flows (Hunt et al. Reference Hunt, Abell, Peterka and Woo1978). It was first introduced to analyse the wall-stress fields of flows around bluff bodies in wind tunnels, but we showed that its application to the complex flow field produced by 3-D DNS (direct numerical simulation) of separated flows around bluff bodies (Dousset & Pothérat Reference Dousset and Pothérat2010, Reference Dousset and Pothérat2012) or complex conduits (Pothérat & Zhang Reference Pothérat and Zhang2018) offered a powerful way to reliably extract the main structures governing the flow’s dynamics. The main idea is to extract 2-D streamlines or stresslines along symmetry planes or solid walls, and to classify the critical points, where both components of velocity or stress cancel, into nodes and saddles. First, the theoretical constraint on the number of saddles and nodes offers a convenient way to validate the numerical simulations (Hunt et al. Reference Hunt, Abell, Peterka and Woo1978). Second, these points offer a simple parametrisation of the entire flow field. When a flow transitions between regimes, topological changes are straightforwardly captured by changes in the number of critical points, and the nature of these points often points to the physical mechanism underpinning the transition. We use it here for the first time in the context of stability analysis, and apply it separately to the base flow and to the perturbations.

Figure 4. Map of the steady velocity magnitude $\Vert \boldsymbol{U} \Vert$ computed for $( N, C_R ) = ( 0.25, 6 )$ at ${\textit{Gr}}_{\textit{ac}} = 6.4\times 10^3$ . The streamlines are displayed together with the base flow critical points: half-saddles (zero wall skin friction, red squares) are labelled as $S'_i$ and $N_i$ are nodes ( $\boldsymbol{U} = \boldsymbol{0}$ , purple discs). The green line depicts the approximate beam radius given by (2.10). The base flow is axisymmetric and is reflected about the $r=0$ axis for the sake of visualisation.

The critical points where $\boldsymbol{U} = \boldsymbol{0}$ are shown in purple (labelled $S_i$ for saddles and $N_i$ for nodes) and points of zero skin friction on the solid boundaries, displayed in red, are half-saddles labelled as $S'_i$ . The base flow field typically comprises a jet aligned with the cavity axis. Being driven by the acoustic forcing, the jet originates at $S'_1$ and impinges the wall facing the transducer. The jet impingement at $S'_8$ then gives rise to a large roll centred on $N_4$ . This roll features weak velocity amplitudes, and spreads over the entire space available between the jet and the lateral cavity wall. The presence of this roll causes the return flow to separate from and reattach to the wall at $S'_5$ and $S'_4$ , respectively. A recirculation bubble centred at $N_3$ lies between these two points.

The distribution and type of critical points within a flow field obey a strict topological rule. More precisely, the number of nodes $\varSigma _{N}$ , of half-nodes $\varSigma _{N'}$ , of saddles $\varSigma _{S}$ and half-saddles $\varSigma _{S'}$ in a slice of $\boldsymbol{U}$ are related to each other through (Hunt et al. Reference Hunt, Abell, Peterka and Woo1978; Foss Reference Foss2004)

(4.1) \begin{equation} 2\varSigma _{N} + \varSigma _{N'} - 2\varSigma _{S} - \varSigma _{S'} = \chi , \end{equation}

where $\chi$ is a parameter defining the topology of the surface in which the streamline patterns are studied ( $\chi = 2$ for the plane shown in figure 4 (Foss Reference Foss2004)). Inspection of the streamlines in figure 4 gives

(4.2) \begin{equation} \varSigma _{N} = 8, \quad \varSigma _{N'} = 0 , \quad \varSigma _{S} = 0 , \quad \varSigma _{S'} = 14 , \end{equation}

so that the left-hand side of (4.1) is indeed equal to $\chi = 2$ . The topological rule (4.1) is thus satisfied, meaning either that the entire flow topology is captured or at least that unresolved parts of the flow are consistently cut off. Indeed, only one out of the infinite series of corner vortices theorised by Moffatt (Reference Moffatt1964) is captured in each corner. Increasing the resolution would reveal a greater number of them. By ‘consistently’, we mean that with each additional vortex, a node and two half-saddles would add to the tally of critical points, without invalidating (4.1).

4.2. Effect of confinement on the base flow topology and on the jet velocity

Reducing the cavity size by either varying $N$ or $C_R$ may affect the 2-D–2-C base flow structure. However, any topological modification of the base flow shall comply with the topological rule (4.1). Figure 5 shows the steady axisymmetric field obtained at the same ${\textit{Gr}}_{\textit{ac}}$ as in figure 4, but for the smaller value of $C_R=4$ . The increased radial confinement yields two additional critical points: a node $N_5$ at $x \approx 66$ and a saddle $S_1$ at $x \approx 72$ . This new critical point combination keeps the left-hand side of the topological rule (4.1) identical to the $C_R=6$ case (figure 4); the base flow for $C_R = 4$ shown in figure 5 is thus again valid from a topological perspective. In addition to these examples, a similar change of the base velocity field topology has only been observed when reducing $C_R$ from 2 to 1 and for forcing magnitudes close to their respective ${\textit{Gr}}_{\textit{ac}}^c$ .

Figure 5. Map of the steady velocity magnitude $\Vert \boldsymbol{U} \Vert$ computed for ( $N = 0.25$ , $C_R=4$ ) for ${\textit{Gr}}_{\textit{ac}} = 6.4\times 10^3$ . The streamlines are displayed together with the base flow critical points: half-saddles (zero wall skin friction, red squares) are labelled as $S'_i$ , $N_i$ are nodes ( $\boldsymbol{U} = \boldsymbol{0}$ , purple discs) and $S_i$ are saddles ( $\boldsymbol{U} = \boldsymbol{0}$ , purple squares). The green line depicts the approximate beam radius given by (2.10). The base flow is axisymmetric and is reflected about the $r=0$ axis for the sake of visualisation.

In addition to potentially modifying the base flow topology, reducing $C_R$ also affects the jet velocity. Figure 6 shows several profiles of the centreline base velocity $U_x (x,r=0 )$ computed for $1 \leqslant C_R \leqslant 6$ and ${\textit{Gr}}_{\textit{ac}} = 6400$ . In a weakly confined setting, the jet strongly accelerates over $0 \leqslant x \lesssim 20$ for all $C_R$ as a result of a balance between $\boldsymbol{F_{\textit{ac}}}$ and the inertial forces (Moudjed et al. Reference Moudjed, Botton, Henry, Ben Hadid and Garandet2014a ; Vincent et al. Reference Vincent, Henry, Kumar, Botton, Pothérat and Miralles2025). Decreasing $C_R$ enhances the jet deceleration at larger $x$ . For $2 \leqslant C_R \leqslant 6$ , this confinement-induced jet deceleration is marginal: $U_x (x,r=0 )$ differs by at most $5.5$ %, and its decrease along the jet results from a balance between $\boldsymbol{F_{\textit{ac}}}$ and viscous forces (Vincent et al. Reference Vincent, Henry, Kumar, Botton, Pothérat and Miralles2025). However, the jet deceleration is significantly stronger for $C_R = 1$ : near the jet impingement, the centreline velocity is 50 % less than for the $C_R \geqslant 2$ cases. This enhanced jet deceleration is caused by the build-up of a streamwise pressure gradient as $C_R$ is decreased (inset of figure 6). This pressure gradient acts against the acoustic forcing and drives the return flow so that the net mass flux through a cross-section is zero (Rudenko & Sukhorukov Reference Rudenko and Sukhorukov1998).

Figure 6. Profiles of the on-axis velocity $U_x(x,r=0)$ along the jet illustrating the effect of flow confinement on the jet velocity for $N=0.25$ and ${\textit{Gr}}_{\textit{ac}} = 6400$ . The profiles are obtained for $1 \leqslant C_R \leqslant 6$ . The dashed vertical line corresponds to the Fresnel distance $L_F$ marking the transition between the near ( $x \lt L_F$ ) and the far ( $x \gt L_F$ ) acoustic fields. Inset shows variations with $C_R$ of the streamwise pressure gradient $\langle \mathrm{d}P / \mathrm{d} x \rangle$ of the base flow at the mid cavity length $L/2$ and averaged over the cross-sectional area.

5.1. Non-oscillatory instabilities in weakly confined settings

We now turn to the linear stability analysis of the steady 2-D–2-C flows similar to those shown in figures 4 and 5. The variations of the growth rate of the leading eigenmode $\sigma$ with ${\textit{Gr}}_{\textit{ac}}$ and $m$ for $N \in \left \{ 0.25, 1 \right \}$ and $4 \leqslant C_R \leqslant 6$ are shown in figure 7. For $4 \leqslant C_R \leqslant 6$ , the same non-oscillatory $m=2$ mode is found to destabilise the flow; we shall denote the associated mode branch as ‘Branch I’. The unstable mode branch remains the same for $4 \leqslant C_R \leqslant 6$ , but ${\textit{Gr}}_{\textit{ac}}^c$ increases from $2191$ for $ ( N, C_R ) = ( 1, 6 )$ to $6237$ for $ ( N, C_R ) = ( 0.25, 4 )$ . Confining the streaming jet flow thus has a stabilising effect, similarly to confined flows past bluff bodies (Mondal & Mahbub Alam Reference Mondal and Alam2023).

Figure 7. Growth rate $\sigma$ of the leading eigenmode as a function of the azimuthal wavenumber $m$ for $3 \leqslant C_R \leqslant 6$ : (a) ( $N, C_{R}$ ) = ( $1,6$ ), (b) ( $N, C_{R}$ ) = ( $0.25,6$ ), (c) ( $N, C_{R}$ ) = ( $0.25,4$ ) (d) ( $N, C_{R}$ ) = ( $0.25,3$ ). Oscillatory and non-oscillatory modes are represented by filled circles and triangles, respectively. The insets show the evolution of $\sigma$ with ${\textit{Gr}}_{\textit{ac}}$ for the indicated value of $m$ . The critical Grashof number ${\textit{Gr}}_{\textit{ac}}^c$ , obtained through spline interpolation of $\sigma$ , is reported in red. The filled black symbols in the inset of panel (b) are additional eigenvalue computations made to improve the interpolation accuracy near $\sigma = 0$ .

Figure 8. Leading mode for $( N, C_R )=( 0.25, 4 )$ at ${\textit{Gr}}_{\textit{ac}} = 6400$ . The mode is non-oscillatory ( $\omega = 0$ ) and unstable ( $\sigma = 1.59 \times 10^{-2} \gt 0$ ). (a) Axial velocity perturbation $u'_x$ in the $( x, r )$ plane, along with the streamlines and critical points of the base velocity $\boldsymbol{U}$ . (b) Velocity perturbation $\boldsymbol{u'}$ , with the background colour corresponding to its azimuthal component $u'_{\theta }$ . (c)–(h) Slices of $\boldsymbol{u'}$ in constant- $x$ planes located by the blue vertical lines in panels (a) and (b), with either $u_x'$ or the streamwise vorticity perturbation $\omega _x'$ as background colour (negative in blue, positive in red). The purple circles locate the points where $U_x = 0$ . In all figures, the green solid line represents the approximate beam radius (2.10) and the blue solid line locates the radius where $U_x$ is 50 % of its on-axis value. All the figures on a given row share the same colour levels.

The topology of leading eigenmodes from Branch I is illustrated for $ ( N, C_R ) = ( 0.25, 4 )$ in figure 8, with the critical points of the base flow overlaid. From the topology in the $ (x,r )$ plane (figures 8 a and 8 b), the perturbation intensity, as measured by $u^\prime _x$ , is strongest in a region stretching from the impingement region of the jet to the shear layer between recirculations centred on $N_4$ and $N_5$ . Radially, the perturbation is located within the acoustic beam (the approximate acoustic beam radius $R_{\textit{beam}}$ is marked by a green line) but remains outside the jet (the points where $U_x$ is 50 % of the on-axis velocity are marked by the blue line). The perturbation then crosses the shear layer at the boundary of the jet to reach the region between the two recirculations around $N_4$ and $N_5$ . Hence, the mechanism is not associated with a shear layer instability, but rather coincides precisely with the strong velocity gradients in the impingement region centred on $S^\prime _8$ , where the flow turns from its incoming streamwise direction to a radial one along the end wall. The azimuthal structure is best illustrated by contours of streamwise vorticity in constant $x$ planes containing either $N_5$ , $S_1$ or $N_4$ (figure 8 f–h) showing a breakdown into an $m=2$ modulation. These figures also clearly show that the streamwise vortices associated with these structures do not follow the shear layer (marked by the purple line corresponding to $U_x=0$ ) and so exclude a shear layer instability mechanism.

Figure 9. Lines of the skin friction perturbation stresses on the downstream wall, with the background colour corresponding to the azimuthal component of the stress vector (positive in red, negative in blue). These lines are shown along with the critical points (blue symbols) at $N = 0.25$ and for different $C_R$ : (a) $C_R = 6$ ; (b) $4$ and (c) $3$ . The approximate beam radius (2.10) is represented by the green circle. The base flow half-saddle $S'_7$ on the downstream wall is located by the solid red circle. The solid blue circle approximates the separation line going through the saddles found at a same distance from the origin. For each case, the mode shapes are computed at forcing magnitudes ${\textit{Gr}}_{\textit{ac}}$ slightly above the instability onset.

To best capture the differences between the modes associated with different branches of instabilities, we use the lines of the skin friction perturbation on the downstream wall and classify their critical points. For case $(N, C_R)= (0.25, 4)$ (figure 9 b), the friction lines include the following numbers of critical points:

(5.1) \begin{align} \varSigma _N = 8 , \quad \varSigma _{N'} = 4 , \quad \varSigma _S = 9 , \quad \varSigma _{S'} = 0. \end{align}

Since $\chi = 2$ for the downstream wall (Foss Reference Foss2004), the topological rule (4.1) is thus again satisfied, validating the mode structure shown in figure 9(b). The topology features attachment and detachment nodes where regions of zero azimuthal stress cause the flow to respectively attach to and separate from the wall (Délery Reference Délery2001). The half-saddle of the base flow $S'_7$ in the meridional plane translates into a circular line separating a ring of Moffat corners in the outmost part of the cavity and an adjacent inner ring of four structures with closed friction lines (red line in the plots of figure 9). The friction lines patterns thus feature four centres in this region (figure 9 a). These points are located between this outer separator and an inner separator, represented by the solid blue circle. In this case, the outer separator is impermeable for both the base flow and the perturbation, i.e. the perturbation field involves no mass flux between the outer region of the Moffatt vortices of the base flow and the main four-vortices structures of the perturbation. Decreasing $C_R$ narrows the gap between the two separators, as seen by comparing the cases with $C_R=6$ and $C_R=4$ , respectively shown in figures 9(a) and 9(b). Decreasing $C_R$ reduces the space available for the friction lines to form closed loops and forms streamwise vortices that are ever more stretched in the azimutal direction.

For $C_R=3$ , another local maximum of $\sigma (m)$ is the first to become unstable at ${\textit{Gr}}_{\textit{ac}}^c=7752$ (figure 7 d). The leading modes of this new branch, which we denote as ‘Branch II’, correspond to non-oscillatory $m=4$ perturbations. This change of unstable branch is the first significant effect of confinement on the flow stability and highlights a minimum $D$ below which the leading perturbation is altered.

To understand the change of leading mode branch, we shall turn again to the patterns of the friction perturbation on the downstream wall. Here, the confinement is increased to the point where the streamwise vortices between the separators of Branch I are stretched to breaking point. This happens as the centres collapse on the inner separator (blue circle), making it impossible to form closed friction loops without crossing the separator. In other words, further reducing $C_R$ necessarily requires the skin friction lines patterns to be rearranged for basic topological rules to not be violated. This topological rearrangement is made possible by a different value of $m$ , corresponding to a new unstable branch. The resulting Branch II mode features a new set of saddles and nodes (figure 9 c). The critical point at $r=0$ is a ‘high-order’ critical point: both the base flow axisymmetry and the $m=4$ symmetry of the mode cause the gradient of the skin friction vector to vanish, meaning that higher-order derivatives are required to classify this point. We refer to the latter as a ‘star saddle’. Being different in nature to classical four-branch saddles of topological weight $w_{4s}=-2$ (i.e. the constant by which $\varSigma _S$ is multiplied in the topological rule (4.1)), this type of critical point must have its own topological weight. Assuming that the topological rule (4.1) is satisfied, the topological weight of this eight-branch star saddle must be $w_{8s}=-6$ .

5.2. Rise of oscillatory instabilities in confined settings

As $C_R$ is further decreased, the modes destabilising the flow become oscillatory. Indeed, the $C_R=2$ flow is destabilised at ${\textit{Gr}}_{\textit{ac}}^c = 10\,513$ by the growth of an oscillatory $m=2$ perturbation (figure 10 a). This mode is associated with a new branch that we label as ‘Branch III’. Interestingly, Branch III and Branch II coexist for this $C_R$ , and the leading Branch II mode becomes unstable at ${\textit{Gr}}_{\textit{ac}} \approx 10\,630 \approx 1.01 {\textit{Gr}}_{\textit{ac}}^c$ . The $\sigma$ of the Branch II mode increases more rapidly with ${\textit{Gr}}_{\textit{ac}}$ than the $\sigma$ of the Branch III mode, so that Branch II becomes dominant for ${\textit{Gr}}_{\textit{ac}} \geqslant 11\,000 \approx 1.04 {\textit{Gr}}_{\textit{ac}}^c$ .

Figure 10. Growth rate $\sigma$ of the leading eigenmode as a function of the azimuthal wavenumber $m$ for (a) $(N, C_R) = (0.25, 2)$ and (b) $( 0.25, 1 )$ . Oscillatory and non-oscillatory modes are represented by filled circles and triangles, respectively. The insets show the evolution of $\sigma$ with ${\textit{Gr}}_{\textit{ac}}$ for the indicated value of $m$ . The critical Grashof number ${\textit{Gr}}_{\textit{ac}}^c$ , obtained through spline interpolation of $\sigma$ , is reported in red. The filled black symbols in the insets of panel (a) are additional eigenvalue computations made to improve the interpolation accuracy in the vicinity of $\sigma = 0$ .

Figure 11. Leading mode for $( N, C_R )=( 0.25, 2 )$ at ${\textit{Gr}}_{\textit{ac}} =10\,600$ . The mode is oscillatory ( $\omega \neq 0$ ) and unstable ( $\sigma = 2.54 \times 10^{-2} \gt 0$ ). (a) Azimuthal velocity perturbation $u'_{\theta }$ in the $( x, r )$ plane, along with the critical points of the base velocity $\boldsymbol{U}$ . (b) Detailed view near the impingement (region framed in red in panel a) of the axial velocity perturbation $u'_x$ together with the base flow streamlines. (c) Details of $\boldsymbol{u'}$ near the impingement, with $u'_{\theta }$ taken as background colour. In all plots, the vertical dashed blue lines locate the slices shown in figure 12. The green solid line represents the approximate beam radius (2.10) and the blue solid line locates the radius where $U_x$ is 50 % of its on-axis value.

Figure 12. Leading velocity perturbation $\boldsymbol{u'}$ for $( N, C_R )=( 0.25, 2 )$ at ${\textit{Gr}}_{\textit{ac}} =10\,600$ . The mode is displayed in constant- $x$ planes located by the vertical dashed blue lines shown in figure 11 and going through the base flow node $N_5$ (slice 1, left) and saddle $S_1$ (slice 2, right). Either the streamwise velocity perturbation $u_x'$ (panels a–b) or the streamwise vorticity perturbation $\omega _x'$ (panels c–d) is used as background colour (negative in blue, positive in red). In all figures, the purple circles locate the points where $U_x = 0$ , the green circle represents the approximate beam radius (2.10) and the blue circle locates the radius where $U_x$ is 50 % of its on-axis value. Figures on a given row share the same colour levels.

For $C_R=2$ , the increased flow confinement and the greater ${\textit{Gr}}_{\textit{ac}}^c$ cause the flow structures resulting from the jet impingement to shrink down and to be further concentrated near the downstream wall (figure 11). In particular, $N_5$ and $S_1$ now lie inside the main forcing region; that was not the case for larger $C_R$ . Therefore, part of the return flow near $N_5$ and $S_1$ now faces the acoustic forcing.

The velocity perturbation field displays distinct features compared with the previous $C_R \geqslant 3$ cases. First, $\vert u'_{\theta } \vert$ is no longer maximum at the impingement, but near $S_1$ instead (figure 11 c). The flow is thus no longer destabilised by the jet impingement, but by the shear layer $N_5$ and $S_1$ instead. This is confirmed by the shape of $\boldsymbol{u'}$ , and the contours of $\omega _x$ at $N_5$ and $S_1$ (figure 12): as for $C_R \geqslant 3$ , the perturbation is still clustered near the shear layer around the jet where $U_x \simeq 0$ . Unlike at lower confinement, the instability extends along a surface closely following this layer. Second, $u'_{\theta }$ displays a wave-like structure in the $x$ -direction, the amplitude of which being significantly reduced as the distance from the downstream wall increases. This is consistent with the findings of Henry et al. (Reference Henry, Vincent, Miralles, Botton and Hadid2022), who reported similar unstable waves for a plane 2-D streaming flow forced by a non-diffracting and non-attenuated beam of square cross-section. These features, reminiscent of a Kelvin–Helmoltz instability, suggest that a high confinement instability originates in the shear layer, rather than at the impingement. Finally, for $C_R = 2$ , the saddles of the skin friction perturbation on the downstream wall no longer coincide with the base flow half-saddle $S'_7$ (figure 13 a). In other words, reducing the cavity radius creates a perturbation flux between the Moffatt roll delimited by $S'_7$ and the remaining of the flow near $x = L$ .

Figure 13. Lines of skin friction stresses exerted by the leading perturbation on the downstream wall for (a) $( N, C_R ) = ( 0.25, 2 )$ and (b) $( N, C_R ) = ( 0.25, 1 )$ . In each plot, the background colour represents the azimuthal component of the stress vector (positive in red, negative in blue). The friction lines are shown along with the zero mode friction critical points (blue symbols) and the approximate beam size (green circle, as defined by (2.10)). The base flow half-saddle $S'_7$ is also reported in red.

Figure 14. Leading mode for $( N, C_R )=( 0.25, 1 )$ at ${\textit{Gr}}_{\textit{ac}} = 15\,500$ . The mode is oscillatory ( $\omega \neq 0$ ) and unstable ( $\sigma = 1.701 \times 10^{-1} \gt 0$ ). (a) Azimuthal velocity perturbation $u'_{\theta }$ in the $( x, r )$ plane, along with the critical points of the base velocity $\boldsymbol{U}$ . (b) Detailed view near the impingement (region framed in red in panel a) of the axial velocity perturbation $u'_x$ together with the base flow streamlines. (c) Details of $\boldsymbol{u'}$ near the impingement, with $u'_{\theta }$ taken as background colour.

Figure 15. Leading velocity pertubation $\boldsymbol{u'}$ for $( N, C_R )=( 0.25, 1 )$ at ${\textit{Gr}}_{\textit{ac}} = 15\,500$ . The mode is shown in constant- $x$ planes located by the blue vertical lines in figure 14, with either $u_x'$ (panels a–b) or the streamwise vorticity perturbation $\omega _x'$ (panels c–d) as background colour (negative in blue, positive in red). The purple circles locate the points where $U_x = 0$ . In all figures, the green solid line represents the approximate beam radius (2.10) and the blue solid line locates the radius where $U_x$ is 50 % of its on-axis value. Figures on a given row share the same colour levels.

For the most confined case ( $C_R=1$ ), a new mode branch arises (figure 10 b). This branch, labelled as ‘Branch IV’, is defined by a leading $m=3$ oscillatory mode. The onset occurs at ${\textit{Gr}}_{\textit{ac}} = 15\,310$ . At ${\textit{Gr}}_{\textit{ac}} \approx {\textit{Gr}}_{\textit{ac}}^c$ , both the enhanced confinement and the greater strength of the return flow resulting from the larger ${\textit{Gr}}_{\textit{ac}}^c$ significantly affect the base flow topology near the downstream wall (figure 14). First, the size of the recirculating flow structures is greatly reduced, e.g. the secondary recirculation bubble on the lateral wall nearly vanishes. Second, the reduced $C_R$ yields two new critical points $N_6$ and $S_2$ in the spatial structure of $\boldsymbol{U}$ . These new critical points appear at approximately the same radial location as $N_5$ and $S_1$ .

Maps of $\boldsymbol{u'}$ in the $ ( x, r )$ plane show that $\vert u'_x \vert$ and $\vert u'_{\theta } \vert$ are maximum near $N_5$ and $S_1$ , i.e. in regions of high shear in the bulk of the base flow (figures 14 b and 14 c). This confirms the trend observed for $C_R = 2$ : confining the flow shifts the locus of the instability from the impingement to the shear layer near $N_5$ and $S_1$ . This shift is also favoured by the build-up of an adverse pressure gradient slowing down the jet (inset of figure 6). In addition, $\boldsymbol{u'}$ is again marginal in the jet (blue line and circles in figure 14) and creates strong $\omega _x$ at the shear layer (figure 15). Finally, the leading Branch IV mode features a ‘star saddle’ critical point at $r=0$ (figure 13 b), as does the leading mode for $C_R=3$ (Branch II).

To conclude, for $N \in \{ 0.25, 1 \}$ and $3 \leqslant C_R \leqslant 6$ , the primary instabilities are caused by non-oscillatory perturbations originating from the jet impingement. Further confining the flow ( $1 \leqslant C_R \leqslant 2$ ) not only relocates the origin of the instability to the shear layer between the jet and the impingement-induced recirculation structures, but also causes the leading perturbation to become oscillatory.

6.1. Model for the nonlinear growth and saturation of the perturbations

To complete the stability study, we shall now determine the criticality of the bifurcations identified with LSA. Our approach relies on the model of Landau & Lifschitz (Reference Landau and Lifschitz1987) describing the growth and saturation of a perturbation near the onset of instability. This model has been extensively used to classify the bifurcations arising in numerous situations, such as the flow past rings (Sheard et al. Reference Sheard, Thompson and Hourigan2004), the wake behind a sphere (Thompson, Leweke & Provansal Reference Thompson, Leweke and Provansal2001) and behind one or multiple cylinders (Henderson & Barkley Reference Henderson and Barkley1996; Henderson Reference Henderson1997; Carmo et al. Reference Carmo, Sherwin, Bearman and Willden2008). The Stuart–Landau model reads

(6.1) \begin{equation} \frac { \mathrm{d} A}{ \mathrm{d} t} = \left ( \sigma + \mathrm{i} \omega \right ) A - l \left ( 1 + \mathrm{i} c \right ) \vert A \vert ^2 A + O(A^5) , \end{equation}

where $A$ is a complex perturbation, and $l$ and $c$ are two real constants. The mechanism responsible for the saturation of $A$ can be easily understood from the sign of $l$ in (6.1). For $l \gt 0$ , the cubic term in the right-hand side of (6.1) acts against the growth of the perturbation: this situation corresponds to a supercritical bifurcation. In contrast, for $l \lt 0$ , higher-order terms are required for $A$ to saturate: the bifurcation is thus subcritical. Therefore, only the sign of $l$ is needed to determine the nature (or criticality) of the bifurcation.

As explained by Sheard et al. (Reference Sheard, Thompson and Hourigan2004) and Kumar & Pothérat (Reference Kumar and Pothérat2020), $l$ may be inferred from a more convenient form of (6.1). Plugging $A(t) = \vert A(t) \vert e^{ \mathrm{i} \phi (t) }$ into the Stuart–Landau model (6.1) gives, after separating the real and imaginary part,

(6.2) \begin{align} \frac { \mathrm{d} \log \vert A \vert }{ \mathrm{d} t} &= \sigma - l \vert A \vert ^2 \, + O ( \vert A \vert ^4 ), \end{align}

(6.3) \begin{align} \frac { \mathrm{d} \phi }{ \mathrm{d} t} &= \omega - l c \vert A \vert ^2 + O ( \vert A \vert ^4 ) . \end{align}

The nature of the bifurcation is determined by plotting $(\mathrm{d} \log \vert A \vert / \mathrm{d} t)$ against $\vert A \vert ^2$ ; the slope of the resulting curve is $-l$ for sufficiently small $\vert A \vert ^2$ and $\sigma$ is obtained from the $y$ -intercept. From this, we compare $\sigma$ with its value obtained from the LSA results and so assess the accuracy of the estimate for $l$ . Note that for a non-oscillatory perturbation, $\phi (t) = 0$ , so only (6.2) is needed.

The saturated state is characterised using the steady-state forms of (6.2) and (6.3). For a supercritical ( $l \gt 0$ ) bifurcation for instance, the first two terms in the right-hand side of the amplitude (6.2) are dominant, so that the saturated perturbation amplitude is

(6.4) \begin{equation} \vert A_{\textit{sat}} \vert = \sqrt { \frac {\sigma }{l} } . \end{equation}

Furthermore, if the saturated state is a time-periodic signal of constant amplitude $\vert A_{\textit{sat}} \vert$ , then $\mathrm{d} \phi / \mathrm{d} t$ reduces to the saturated frequency $\omega _{\textit{sat}}$ . Using (6.4), the steady-state version of (6.3) gives, for the Landau constant $c$ ,

(6.5) \begin{equation} c = \frac { \omega - \omega _{\textit{sat}} }{\sigma } , \end{equation}

thus allowing for all the quantities in the truncated Stuart–Landau model (6.1) to be fully characterised.

We shall now define $A$ in terms of flow quantities. We followed Kumar & Pothérat (Reference Kumar and Pothérat2020), and based $\vert A \vert$ on a local value of a flow variable and, more particularly, on $u'_{\theta }(t)$ . Other choices are possible, including global measures of the flow unsteadiness (Thompson et al. Reference Thompson, Leweke and Provansal2001; Sheard et al. Reference Sheard, Thompson and Hourigan2004; Sherwin & Blackburn Reference Sherwin and Blackburn2005; Sapardi et al. Reference Sapardi, Hussam, Pothérat and Sheard2017). For oscillatory bifurcations, we considered the envelope of $u'_{\theta }(t)$ . Choosing where to record $u'_{\theta }(t)$ is based on the need to obtain clean signals. To maximise the signal-to-noise ratio, we thus chose points where LSA predicted large $\vert u'_{\theta } \vert$ . Depending on the values of $C_R$ , these points are located near the critical points $N_5$ and $S_1$ , and near the impingement. The precise locations of these points are given in table 4.

Table 4. Coordinates $( x, r, \theta )$ of the points where the time series of the azimuthal velocity perturbation $u'_{\theta }$ are recorded in the nonlinear 3-D–3-C simulations.

From a numerical point of view, the nonlinear 3-D–3-C unsteady simulations were run for ${\textit{Gr}}_{\textit{ac}}$ values slightly larger than ${\textit{Gr}}_{\textit{ac}}^c$ . The relative gap to ${\textit{Gr}}_{\textit{ac}}^c$ is defined by the criticality parameter $r_c$ :

(6.6) \begin{equation} r_c = \frac {{\textit{Gr}}_{\textit{ac}}}{{\textit{Gr}}_{\textit{ac}}^c} - 1 , \end{equation}

so that $r_c \gt 0$ corresponds to a ${\textit{Gr}}_{\textit{ac}} \gt {\textit{Gr}}_{\textit{ac}}^c$ flow regime and $r_c \lt 0$ to a ${\textit{Gr}}_{\textit{ac}} \lt {\textit{Gr}}_{\textit{ac}}$ regime. Finally, $(\mathrm{d} \log \vert A \vert / \mathrm{d} t)$ was computed from the recorded time series using centred finite differences.

6.2. Topology of the perturbed states from the nonlinear 3-D–3-C unsteady simulations

We first assessed the relevance of the LSA to the full nonlinear solution. For that purpose, we plotted the contours of the azimuthal velocity perturbation $u_\theta ^\prime$ for four representative cases, one per instability branch. These were obtained slightly above onset, by susbstracting the base axisymmetric solution from the full nonlinear solution. The results are illustrated in figure 16. Comparing with the topology of the leading eigenmodes in figures 8, 9, 11–15, confirms that the topologies of the linear and nonlinear solutions are practically identical. See also the comparison between the leading oscillatory mode from LSA at $(N,C_R)=(0.25,1)$ and the travelling wave that results from this mode in the 3-D–3-C nonlinear simulations in figure 21. In all cases, the leading eigenmode comes very close to the full nonlinear solution. In particular, the nonlinear 3-D–3-C unsteady simulations recover precisely the different branches of instability found with the LSA and their corresponding azimuthal wavenumbers. We shall see in the next section that oscillating frequencies are also closely recovered.

Figure 16. Nonlinear 3-D–3-C unsteady simulations of the velocity perturbation represented by contours of azimuthal velocity perturbation $u^\prime _\theta$ set at $\pm 0.2\textrm {max}{|u^\prime _\theta |}$ for the leading eigenmode corresponding to unstable branches (a) IV ( $N=0.25,\,C_R=1$ ), (b) III ( $N=0.25$ , $C_R=2$ ), (c) II ( $N=0.25$ , $C_R=3$ ) and (d) I ( $N=1$ , $C_R=6$ ). The black surface represents the position of the approximate beam radius given by (2.10).

6.3. Nature of the bifurcations

The Stuart–Landau analysis has been carried out for all the set-ups listed in table 1; the cases defined by $N=0.25$ and $C_R \in \left \{ 1, 2, 4 \right \}$ are shown in figure 17 as examples. For each case, the initial instants of $\vert A \vert$ are dominated by noise (figures 17 a, 17 c and 17 e). The Stuart–Landau equation (6.2) for $\vert A \vert$ is then fitted to the time series over a time interval free of the initial noise; that range is highlighted in red in all plots of figure 17.

Figure 17. Stuart–Landau analysis for (a)–(b) $(N = 0.25, C_R = 4)$ at ${\textit{Gr}}_{\textit{ac}} = 6400$ ( $r_c = 0.0262$ ), (c)–(d) $(N = 0.25, C_R = 2)$ at ${\textit{Gr}}_{\textit{ac}} = 10\,600$ ( $r_c = 0.0083$ ) and (e)–(f) $(N = 0.25, C_R = 1)$ at ${\textit{Gr}}_{\textit{ac}} = 15\,500$ ( $r_c = 0.0124$ ). (a,c,e) Growth and saturation of a white-noise-triggered perturbation injected at $t=0$ , monitored by the time series of the azimuthal velocity perturbation $u_{\theta }'$ . (b,d,f) Evolution of $\mathrm{d} ( \log \vert A \vert ) / \mathrm{d}t$ with $\vert A \vert ^2$ , where $\vert A \vert$ is the amplitude of the perturbation plotted in panels (a,c,e) (for the $C_R = 2$ and $1$ cases, $\vert A \vert$ is based on the envelope of $u_{\theta }'$ ). The growth rate $\sigma$ and the Landau coefficient $l$ are obtained by fitting (6.2) to the red curve. For each case, the fitting range is highlighted in red.

For $4 \leqslant C_R \leqslant 6$ , $l \gt 0$ near the instability onset. The leading Branch I mode thus destabilises the flow through a supercritical circular pitchfork bifurcation (Touihri, Ben Hadid & Henry Reference Touihri, Ben Hadid and Henry1999). Further confining the flow yields $l \lt 0$ for $C_R=3$ , indicating that an unstable Branch II mode triggers a subcritical circular pitchfork bifurcation. Although the leading mode becomes oscillatory for $C_R=2$ , $l$ is negative (figure 17 d): the associated Hopf bifurcation is thus subcritical. Finally, $l \gt 0$ for $C_R=1$ (figure 17 f); hence, the unstable Branch IV mode is associated with a supercritical Hopf bifurcation. For all cases, the difference between the $\sigma$ obtained with the nonlinear 3-D–3-C simulations and LSA is at most 6 % (table 5). This excellent agreement (i) validates the LSA and (ii) adds further confidence on the validity of the Stuart–Landau analyses.

Finding subcritical bifurcations for $C_R \in \left \{ 2, 3 \right \}$ raises the question of whether the corresponding instabilities can be triggered for $r_c \lt 0$ . For both cases, introducing white noise for $- 0.05 \leqslant r_c \leqslant -0.01$ did not cause the flow to deviate from its axisymmetric base state. Triggering these bifurcations at ${\textit{Gr}}_{\textit{ac}} \lt {\textit{Gr}}_{\textit{ac}}^c$ would require different tools, e.g. the analysis of the transient growth of non-modal perturbations (Schmid Reference Schmid2007). Finding the minimum seeds to trigger these subcritical bifurcations is nevertheless out of the scope of the present work and shall be the subject of a dedicated study.

6.4. Saturated states and secondary instabilities

We shall finally characterise the saturated states following the perturbation growth at $r_c \gt 0$ .

Table 5. Comparison between the growth rates $\sigma$ and frequencies $\omega$ returned by linear stability analysis (LSA) with those obtained from nonlinear 3-D–3-C unsteady simulations of ${\textit{Gr}}_{\textit{ac}} \gt {\textit{Gr}}_{\textit{ac}}^c$ regimes ( $r_c \gt 0$ , see 6.6). The relative errors between the values of $\sigma$ and $\omega$ obtained by LSA and the nonlinear 3-D–3-C simulations are given by $\varepsilon _{\sigma }$ and $\varepsilon _{\omega }$ , respectively.

Figure 18. Frequencies obtained from the time-series of the azimuthal velocity perturbation $u'_{\theta }(t)$ for (a) $(N=0.25, C_R=1)$ , $\textit{Gr}_{\textit{ac}}=15\,500\,(r_{c}=0.0124)$ and (b) $(N=0.25, C_R=2)$ , $\textit{Gr}_{\textit{ac}}=10\,600{}(r_{c}=0.0083)$ . The corresponding time-series are those displayed in figures 17(e) and 17(c), respectively. The frequencies are computed from FFTs made in the linear regime ( $\omega$ , blue) and in the saturated regime ( $\omega _{\textit{sat}}$ , red).

For Hopf bifurcations, the frequency $\omega _{\textit{sat}}$ of the saturated perturbation may differ from $\omega$ in the linear regime (regime of exponential perturbation growth), yielding a non-zero $c$ (6.5). A change of frequency is captured by taking the fast Fourier transform (FFT) of $u'_{\theta }(t)$ in both the linear and the saturated regimes (figure 18).

For $ ( N, C_R ) = (0.25, 2 )$ , $\omega _{\textit{sat}}$ differs by 41.8 % from $\omega$ (figure 18 b). In that case, $c = 44.83 \pm 14.63$ , and the discrepancy between the values of $\omega$ predicted by LSA and the nonlinear 3-D–3-C simulations is 0.30 % (table 5). For $ (N, C_R ) = (0.25, 1 )$ , the $\omega$ and $\omega _{\textit{sat}}$ obtained with the nonlinear 3-D–3-C simulations in the linear and saturated regimes, respectively, are almost identical (figure 18 a). They also barely differ from the frequency predicted by LSA: the observed discrepancies are 0.29 % (table 5) and 0.87 % in the linear and saturated regimes, respectively. Given these small differences, $\omega$ and $\omega _{\textit{sat}}$ are considered to be equal; hence, $c = 0$ .

For $( N, C_R) = (0.25, 2)$ , the large uncertainty of $c$ is due to the rather low resolution of the saturated regime spectrum (figure 18 a). That low resolution results from the short duration of the saturated regime: the latter indeed lasts approximately 10 oscillation periods, after which $u'_{\theta }$ changes abruptly (figure 19). The saturated state is thus itself unstable and causes the flow to transition towards a highly unsteady and chaotic flow regime. Unfortunately, that secondary instability develops too abruptly to be characterised by means of a Stuart–Landau analysis.

A secondary instability was also observed for $ (N, C_R ) = (0.25, 3 )$ at $r_c = 0.0320$ (figure 20 a). For that case, the saturated state is unstable to a non-oscillatory perturbation. In contrast to the $ (N, C_R ) = (0.25, 2 )$ case, the slow deviation from the saturated state allows for a Stuart–Landau analysis to be made (figure 20 b). As the primary bifurcation, the secondary instability is found to be subcritical. However, it remains unclear whether the temporal variations of $u_{\theta }'$ observed for $t \geqslant 110$ correspond to a transient regime before reaching a stable steady state or, instead, are the beginning of a chaotic regime. Significantly longer simulations would be necessary to answer this question.

Figure 19. Time-series of the azimuthal velocity perturbation $u_{\theta }'$ for $( N, C_R )=( 0.25, 2 )$ and ${\textit{Gr}}_{\textit{ac}} = 10\,600$ ( $r_c = 0.0083$ ). The signal is recorded at $(x, r, \theta ) = ( 81, 1.3, 0 )$ .

Figure 20. Stuart–Landau analysis for $( N, C_R )=( 0.25, 3 )$ at ${\textit{Gr}}_{\textit{ac}} = 8000$ ( $r_c=0.0320$ ). (a) Temporal evolution of the azimuthal velocity perturbation $u_{\theta }'$ obtained from nonlinear 3-D–3-C unsteady simulations and recorded at $( x, r, \theta ) = ( 80, 2, 0 )$ . The perturbation originates from white noise injected at $t=0$ . (b) Evolution of $\mathrm{d} ( \log \vert A \vert ) / \mathrm{d}t$ with $\vert A \vert ^2$ (red curve), where $\vert A \vert = \vert u'_{\theta } \vert$ and is extracted from panel (a). The black dashed line is obtained by fitting the Stuart–Landau equation (6.2) to the red curve over the range highlighted in red. The positive slope indicates a subcritical transition.

Figure 21. Comparison between the leading velocity perturbation $\boldsymbol{u'}$ obtained with LSA (left) and the perturbation obtained with nonlinear 3-D–3-C unsteady simulations (right) for (a) $( N, C_R ) = ( 0.25, 1 )$ and (b) $( N, C_R ) = ( 0.25, 2 )$ . The results from LSA and nonlinear 3-D–3-C unsteady simulations are compared at the location of the axisymmetric base flow saddle $S_1$ ( $x=80.15$ in panel a and $x=77.52$ in panel b). The indicated $\sigma$ and $\omega$ are those of LSA and the nonlinear results are extracted from the early exponential growth of the perturbation. Features of the steady base velocity $\boldsymbol{U}$ are reported in both plots: the points where the sign of the streamwise velocity $U_x$ changes (purple circles) and the points where $U_x$ is equal to 50 % of its on-axis value (blue circle). The green circle represents the approximate edge of the beam defined by (2.10). See the supplementary movies (movies 1 and 2) for the animations of the velocity perturbations obtained from the nonlinear 3-D–3-C simulations.

Finally, the leading oscillatory modes for $( N, C_R ) = ( 0.25, 1 )$ and $( 0.25, 2 )$ are defined by a pair of complex conjugate eigenvalues corresponding to counter-propagating waves in the $\boldsymbol{e_{\theta }}$ -direction. Depending on the complex amplitudes of these unstable modes, their combination results in either a standing wave or a travelling wave (Clune & Knobloch Reference Clune and Knobloch1993). Although these complex amplitudes cannot be obtained from LSA, they do appear in the nonlinear 3-D–3-C unsteady simulations. The leading LSA mode shapes are compared with snapshots of the 3-D–3-C unsteady simulations in figure 21 for $ ( N, C_R ) = ( 0.25, 1 )$ and $ ( 0.25, 2 )$ . The snapshots are taken from the early stage of the exponential growth of the perturbation to avoid nonlinear effects on the perturbation shape. For both cases, a strong agreement is found between the Branches III and IV mode shapes obtained by LSA and by the nonlinear 3-D–3-C simulations, and both perturbations correspond to waves travelling along $\boldsymbol{e}_{\theta }$ (see the supplementary movies 1 and 2 available at https://doi.org/10.1017/jfm.2025.10810).

To conclude, the thorough characterisation of the Stuart–Landau model (6.1) for all cases (table 6) showed that the sub- or supercritical nature of the primary bifurcation is closely related to the flow confinement. For the subcritical primary bifurcations, the saturated states were found to be unstable, causing the flow to transition to an unsteady state.

Table 6. Summary of the Stuart–Landau analysis based on local measurements of $u'_{\theta }$ ; these measurements are obtained from nonlinear 3-D–3-C unsteady simulations of ${\textit{Gr}}_{\textit{ac}} \gt {\textit{Gr}}_{\textit{ac}}^c$ regimes ( $r_c \gt 0$ , see (6.6)). The $l$ estimates are bounded by their 95 % confidence interval. The confidence intervals for $c$ further involve the resolution of the frequency spectra from which $\omega$ and $\omega _{\textit{sat}}$ are determined. See table 4 for the coordinates of the points where $u'_{\theta }$ is recorded.