# Stability of acoustic streaming jets confined in cylindrical cavities

Bjarne Vincent, Abhishek Kumar, Daniel Henry, Sophie Miralles, Valéry

Botton, Alban Pothérat

► To cite this version:

Bjarne Vincent, Abhishek Kumar, Daniel Henry, Sophie Miralles, Valéry Botton, et al. Stability of acoustic streaming jets confined in cylindrical cavities. 2025. hal-05322822

HAL Id: hal-05322822

https://hal.science/hal-05322822v1

Preprint submitted on 20 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Banner appropriate to article type will appear here in typeset article

1

2

# Stability of acoustic streaming jets confined in cylindrical cavities

3

4

Bjarne Vincent¹,², Abhishek Kumar², Daniel Henry¹, Sophie Miralles¹†, Valéry Botton¹, and Alban Pothérat²‡

5

¹INSA Lyon, CNRS, Ecole Centrale de Lyon, Universite Claude Bernard Lyon 1, Laboratoire de

6

Mécanique des Fluides et d'Acoustique, UMR5509, 69621, Villeurbanne France

7

²Fluid and Complex Systems Research Centre, Coventry University, Coventry CV1 5FB, UK

8

(Received xx; revised xx; accepted xx)

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

We study the stability of a steady Eckart streaming jet flowing in a closed cylindrical cavity. This configuration is a generic representation of industrial processes where driving flows in cavity by means of acoustic forcing offers a contactless way of stirring or controling flows. Successfully doing so however requires sufficient insight into the topology induced by the acoustic beam. This, in turn, raises the more fundamental question of whether the basic jet topology is stable and, when it is not, of the alternative states that end up being acoustically forced. To answer these questions we consider a flow forced by an axisymmetric diffracting beam of attenuated sound waves emitted by a plane circular transducer at one cavity end. At the opposite end, the jet impingement drives recirculating structures spanning nearly the entire cavity radius. We rely on Linear Stability Analysis (LSA) together with three-dimensional nonlinear simulations to identify the flow destabilisation mechanisms and to determine the bifurcation criticalities. We show that flow destabilisation is closely related to the impingement-driven recirculating structures, and that the ratio $C_R$ between the cavity and the maximum beam radii plays a key role on the flow stability. In total, we identified four mode types destabilising the flow. For $4 \le C_R \le 6$, a non-oscillatory perturbation rooted in the jet impingement triggers a supercritical bifurcation. For $C_R = 3$, the flow destabilises through a subcritical non-oscillatory bifurcation, and we explain the topological change of the unstable perturbation by analysing its critical points. Further reducing $C_R$ increases the shear within the flow, and gradually moves the instability origin to the shear layer between impingement-induced vortices: for $C_R = 2$, an unstable travelling wave grows out of a subcritical bifurcation, which becomes supercritical for $C_R = 1$. For each geometry, the nonlinear 3D simulations confirm both the topology and the growth rate of the unstable perturbation returned by LSA. This study offers fundamental insight into the stability of acoustically-driven flows in general but also opens possible pathways to either induce turbulence acoustically, or to avoid it in realistic configurations.

Key words: —

<footnote>† Email address for correspondence: sophie.miralles@insa-lyon.fr</footnote>

<footnote>‡ Email address for correspondence: alban.potherat@coventry.ac.uk</footnote>
35

36

# 1. Introduction

37

38

39

40

This work is concerned with the stability of a jet flow driven by an axisymmetric beam of attenuated travelling sound waves in a closed cavity. We focus here on the effect of confinement on the stability of the jet, and shed light on the very unstable nature of such flows.

41

42

43

44

45

46

47

Sound-driven flows are referred to as acoustic streaming, and result from a nonlinear process producing momentum out of the attenuation of sound waves. Acoustic streaming can be divided into two categories, depending on whether the flow is forced by either standing (Rayleigh-Schlichting streaming (Rayleigh 1884; Westervelt 1953)) or travelling (Eckart streaming (Eckart 1948)) sound waves. Since this work focuses exclusively on flows driven by travelling sound waves, any mention of acoustic streaming from now shall refer to Eckart streaming, unless specified otherwise.

48

49

50

51

52

53

54

As they allow for remote flow creation, both Rayleigh and Eckart streaming recently received a renewed interest in fields for which geometrical constraints prevent fluid actuation by mechanical means, e.g., using a propeller. This is typically the case in microfluidics, where Rayleigh and Eckart streaming are commonly used to promote flow in microcavities and microchannels (Hagsäter et al. 2007; Frommelt et al. 2008; Friend & Yeo 2011; Muller & Bruus 2015), in droplets (Alghane et al. 2012; Riaud et al. 2017), and around oscillating bubbles (Cleve et al. 2019; Doinikov et al. 2022; Fauconnier et al. 2022).

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

At a larger scale, acoustic streaming offers an advantageous mean of actuating fluids whilst avoiding the contamination inherent to mechanical contact. This situation commonly arises in metals and semiconductors solidification processes, during which flows in the melt may affect the quality of the solid ingot. By stabilising convective flows within the melt (Dridi et al. 2008b), acoustic streaming is known to enhance the homogenity of the solid ingots (Kozhemyakin 2003). Solidification processes may also benefit from a local control of the grain growth obtained by directing a streaming jet towards the solidification front (Lebon et al. 2019). More recently, Absar et al. (2017) showed that nanocomposites manufacturing can also benefit from acoustic streaming to homogenise the distribution of fibers and particles within the metal matrix. All these applications of acoustic streaming to materials manufacturing rely on low-frequency and high-power ultrasounds (typically a few kilohertz and hundreds of watts) which, close to the source, create a cloud of cavitation bubbles acting as the main sound attenuation mechanism driving the flow. Although this setting leads to strong jets, it suffers from the loss of beam coherence at short distances from the source. This prevents any remote control of the flow. We are instead interested in using beams of low-power and megahertz sound waves to extend the range over which controlled flow actuation is possible (Kamakura et al. 1996; Mitome 1998). By acting deep into the bulk of the flow, such a choice of parameters may open new possibilities of using acoustic streaming in solidification processes. Examples include the local enhancement of mass transfer at the solidification front with an elongated jet (El Ghani et al. 2021), and the creation of three-dimensional flows using multiple reflections of a single beam to stir the melt (Vincent et al. 2024b). Still, whether the purpose is to drive a complex chaotic flow (Cambonie et al. 2017; Launay et al. 2019; Qu et al. 2022) or to stabilise the melt (Dridi et al. 2008b), these applications require a detailed understanding of the flow patterns for a given configuration and beam intensity. Hydrodynamic instabilities may lead to the formation of flow patterns with a topology that significantly differs from that of the forcing beam and even to turbulence. Whether such states are desired or not, understanding how a streaming jet may destabilise is key to designing optimal industrial setups involving acoustic streaming.
83

84

85

86

87

88

89

90

91

92

93

The destabilisation of Eckart streaming jets has been observed in several experiments, such as in those of Green et al. (2016). The authors studied a streaming flow in a vertical cylindrical cavity filled with a mixture of diethyl phtalate and ethanol. The circular transducer, placed at the top endwall of the cavity, fired bursts of ultrasonic waves travelling towards the opposite end of the fluid volume, where a layer of sound-absorbing material prevented the backward reflection of the incident sound waves. This setting created a steady streaming flow, which however transitioned towards another state after some time due to a thermal instability arising from the sound-absorbing plate progressively heating up. Moudjed et al. (2014b) also reported strong distortion of the streaming jet when increasing the forcing magnitude. These severe and intermittent jet oscillations systematically occurred near the wall facing the transducer, suggesting that the impingement plays a major role on the flow stability.

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

These experimental observations of acoustic streaming destabilisation motivated theoretical stability analyses. The first acoustic streaming stability work traces back to Dridi et al. (2010), who considered an analytical 1D-1C (one-dimensional, one-component) jet velocity profile in a 2D infinite layer. For various ratios of the source height to the layer thickness, the authors showed that the flow is destabilised by a wave travelling in the direction facing the acoustic forcing. The frequency of that wave decreases as the layer height is increased with respect to the source height. Later, multidimensional simulations in a closed cavity revealed that the destabilisation is indeed rooted in the jet impingement (Ben Hadid et al. 2012), as observed by (Moudjed et al. 2014b). Since then, numerous studies focused on how acoustic streaming may stabilise flows driven by other phenomena such as lateral (Dridi et al. 2008a) and vertical (i.e., Rayleigh-Bénard) heating (Henry et al. 2022), gradient of chemical concentration (Lyubimova & Skuridin 2019) and thermodiffusion (Charrier-Mojtabi et al. 2012, 2019). All these studies, however, suffer from a rather crude modelling of the acoustic field driving the flow. First, none of these works considered beam diffraction. The radial enlargement of the beam not only prescribes a radial length scale to the jet (Moudjed et al. 2014a; Vincent et al. 2024a), but also locally reduces the magnitude of the acoustic forcing due to the conservation of acoustic power along the beam. As such, these two effects create longitudinal gradients of the streamwise jet velocity (Kamakura et al. 1996; Mitome 1998; Dentry et al. 2014) that cannot be recovered in a non-diffracting approximation. Second, none of the former stability studies accounted for the decay of the forcing along the beam caused by sound attenuation. This approximation violates two fundamental properties of streaming jets: First, attenuation causes the streamwise jet momentum flux to saturate at infinite distances from the source. Second, attenuation defines the length scale over which the jet builds up momentum (Lighthill 1978). As a consequence, this non-attenuated approximation actually severely overestimates the jet velocity, and fails to capture the jet's structure. Therefore, the causes of destabilisation of realistic jets forced by a diffracting beam of attenuated sound waves still remains an open question.

121

122

123

124

125

126

127

128

129

130

131

132

The aim of this work is precisely to determine how a laminar streaming jet driven by an attenuated and diffracting ultrasound beam may become unstable. We shall consider the simplest geometry for this purpose, i.e. a closed cylindrical cavity with impermeable and sound-absorbing boundaries. The acoustic source, placed at one cavity end, radiates an attenuated beam driving a swirl-free axisymmetric jet. The key novelty is to use a realistic model for the jet accounting for both attenuation and diffraction based on our recent work (Vincent et al. 2024a). With these models for the cavity and acoustic beam, the problem loses streamwise invariance. The base flow must be determined by numerical simulations and its bi-global stability is carried out numerically too. To assess the influence of confinement, we consider different cavity sizes based on the attenuation length of the sound waves and on the maximum beam diameter. We address the problem by means of numerical methods based on high-order spectral elements to answer the following questions:
<figure><img src="image_2.png" /></figure>

Figure 1: Sketch of the investigated setup 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 decay 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 u. All the boundaries are no-slip walls and completely absorb the acoustic waves.

(i) What are the mechanisms underpinning the flow's instability and what are the conditions of their onset?

(ii) What is the type (oscillatory or not) of the unstable perturbation, and how does its topology change when varying the cavity size with respect to the size of the forcing field?

(iii) What is the nature of the bifurcations leading to these modes, and in particular, are they subcritical or supercritical?

(iv) What are the nonlinear states arising out of the instability?

This work is organised as follows: we shall first define the problem along with the governing equations in § 2, before presenting in § 3 the computational methods together with the grid convergence analyses. We then analyse the effect of confinement on the flow structure in § 4, before turning to the LSA in § 5. Then, we shall determine the bifurcation criticalities for all cavity sizes in § 6. Concluding remarks are given in § 7.

## 2. Formulating the stability problem for jets driven by attenuated and diffracting progressive waves

### 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 $e_x$, $e_r$ and $e_θ$ refer to the axial, radial and azimuthal directions, respectively. The spatial coordinates in the corresponding directions are x, r and θ; 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.* 2024a). In particular, all lengths are normalised by the dimensional source diameter $D_s$, and we shall use $ρD_s^3$ as mass scale (where ρ is the dimensional fluid density) and $D_s^2/v$ as time scale (where v stands for the dimensional kinematic viscosity). From now on, we shall only rely on nondimensional parameters based on these scales.

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 $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 **u**. The injected momentum is scaled by M, the maximum dimensionless momentum flux to be acquired by the jet if L → ∞. 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 5 nondimensional parameters that are defined in terms of dimensional parameters as follows (Vincent et al. 2024a):

$$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_{ac}}{\rho c v^2},$$

where L_c and D_c are the dimensional cavity length and diameter, respectively. The dimensional acoustic pressure attenuation coefficient is denoted by α and the dimensional speed of sound by c. The dimensional acoustic power radiated by the transducer is denoted by P\_ac and λ is the dimensional wavelength.

While the acoustic field is characterised by small length scales and fast time scales, the streaming flow associated to them is significantly slower than the speed of sound and involves length scales that greatly exceed the acoustic wavelength (Lighthill 1978; Orosco & Friend 2022). These scale differences allow to model the mean motion as that of an incompressible flow due to an external momentum source (Moudjed et al. 2014a):

$$\nabla \cdot \mathbf{u} = 0, \qquad (2.1)$$

$$\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla) \mathbf{u} = -\nabla p + \nabla^2 \mathbf{u} + F_{ac}, \qquad (2.2)$$

$$F_{ac} = Gr_{ac} \tilde{I}. \qquad (2.3)$$

Here, p and u respectively refers 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, F\_ac, couples the acoustics and the hydrodynamics problems. The spatial structure of F\_ac is set by the normalised acoustic intensity $\tilde{I}$, which is bounded by unity for an unattenuated sound field. The magnitude of the forcing is defined by the Grashof number

$$Gr_{ac} = \frac{32\alpha P_{ac} D_s}{\pi \rho c v^2} = \frac{32 N M}{\pi L}$$

comparing the magnitude of the acoustic forcing F\_ac to the magnitude of the viscous forces at the scale of the transducer (El Ghani et al. 2021).

The cavity walls impose the no-slip boundary condition:

$$u = 0$$

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

$$u_{\theta} = \frac{\partial u_{x}}{\partial \theta} = \frac{\partial u_{r}}{\partial \theta} = \frac{\partial p}{\partial \theta} = 0,$$

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

$$u_r = \frac{\partial u_x}{\partial r} = \frac{\partial p}{\partial r} = 0$$
203 on the axis ($r = 0$).

204 To model the forcing field $F_{ac}$, we assume that:

205 (i) The length and time scales associated with the acoustic problem are significantly smaller than those of the streaming flow, as verified in the experiments of (Kamakura et al. 1996; Mitome 1998; Moudjed et al. 2014a). As a result, both the streaming flow velocity levels and gradients do not alter the propagation of sound (Pierce 1990). Therefore, $F_{ac}$ is time-independent.

209 (ii) We consider the acoustic waves to be irrotational. This common assumption in acoustics is known to hold in the bulk of the fluid (Coulouvrat 1992; Kinsler et al. 2000; Blackstock 2000); however, sound waves are no longer irrotational in a thin viscous layer called 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. (2014a), on which we shall base the structure of $F_{ac}$, the ABL thickness is $\delta \approx 1.3 \times 10^{-5}$. As $\delta$ is several orders of magnitudes 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. 2010).

210 (iii) We assume that 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. (2014a), and the acoustic nonlinearity parameter (Hamilton & Blackstock 2024). 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) < L_s$.

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

220 (v) Sound attenuation is due to the viscosity of the medium only. Coulouvrat (1992) 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. 2017) (this approximation however does not hold for liquid metals). We also discard any sound attenuation caused by relaxation.

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

230 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. 2024a). We showed that $\tilde{I}$ can be obtained from a set of integrals defining the acoustic pressure $\tilde{p}_{ac}$ and velocity $\tilde{u}_{ac}$ at any point $(x, r)$ in the fluid domain:

$$ \tilde{I} = \mathfrak{R} \{ \tilde{p}_{ac} (\tilde{u}_{ac}^*) \}, \quad (2.4a) $$

$$ \tilde{p}_{ac} = \frac{k^2 S}{8 \times 1.22\pi^2} \int_0^{1/2} \int_0^{2\pi} \frac{e^{-ik\|d\|}}{\|d\|} r' dr' d\theta', \quad (2.4b) $$

$$ \tilde{u}_{ac} = \frac{1}{4\pi} \int_{0}^{1/2} \int_{0}^{2\pi} (1 + ik\|d\|) \frac{e^{-ik\|d\|}}{\|d\|^2} \frac{d}{\|d\|} r' dr' d\theta', \quad (2.4c) $$

$$ k = \frac{2.44\pi}{S} - i\frac{N}{L}, \quad (2.4d) $$

243 where * refers to the complex conjugate and $\mathfrak{R}$ to the real part. In equations (2.4), $d = x e_x + [r - r' \cos(\theta')] e_r - r' \sin(\theta') 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, $\tilde{I}$ is mainly shaped by the complex wave number $k$: the real part
<table><caption>Table 1: Values of the parameters defining each setup. 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 equation (2.5). 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. (2014a).</caption><thead><tr><th>C<sub>R</sub></th><th>N</th><th>S</th><th>L</th><th>D</th></tr></thead><tbody><tr><td>1</td><td>0.25</td><td>0.03</td><td>83.3</td><td>6.02</td></tr><tr><td>2</td><td>0.25</td><td>0.03</td><td>83.3</td><td>12.04</td></tr><tr><td>3</td><td>0.25</td><td>0.03</td><td>83.3</td><td>18.06</td></tr><tr><td>4</td><td>0.25</td><td>0.03</td><td>83.3</td><td>24.08</td></tr><tr><td>6</td><td>0.25</td><td>0.03</td><td>83.3</td><td>36.12</td></tr><tr><td>6</td><td>1</td><td>0.03</td><td>333.3</td><td>126.51</td></tr></tbody></table>

of $k$ sets the beam diffraction, whilst the imaginary part of $k$ involves the sound attenuation coefficient $N/L$ and thus adds attenuation to $\tilde{I}$.

Equations (2.4) 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_{beam}$ is:

$$R_{beam} = 0.5 + Sx, \qquad (2.5)$$

from which we define the confinement ratio $C_R$:

$$C_R = \frac{D/2}{R_{beam}|_{x=L}}. \qquad (2.6)$$

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 (respective 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. (1996), Mitome (1998) and Moudjed et al. (2014a) for which $C_R = 2.9, 8.5$ and $4.1$, respectively. As these experiments were made for $N < 1$, we set $N = 0.25$ for five cases. Nevertheless, from a process standpoint, it is preferable to have $N \ge 1$ to ensure that most of the acoustic kinetic energy radiated by the transducer is used to create flow kinetic energy (Vincent et al. 2024a). 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. (2014a) and Moudjed et al. (2015). Therefore, the structure of $F_{ac}$ in our work is consistent with the forcing fields of former experiments. These experiments covered a wide range of $Gr_{ac}$; it typically ranged from $5.6 \times 10^2$ (Mitome 1998) to $1.2 \times 10^5$ (Moudjed et al. 2014a). In the present work, we shall adjust $Gr_{ac}$ on a case-by-case basis to identify the instability onset.

## 2.2. Linear stability analysis

Although the base flow is 2D-2C (two-dimensional, two-component), it may nevertheless be destabilised by non-axisymmetric 3D-3C 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:

$$u = U(x, r) + u'(x, r, \theta, t), \qquad (2.7)$$
279

$$p = P(x, r) + p'(x, r, \theta, t). \qquad (2.8)$$

280

281

282

283

284

285

The perturbations ($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.7) and (2.8) into equations (2.2) and (2.1) describing the mean flow, and by linearising the resulting equations around the base flow. Since ($U, P$) satisfy equations (2.2) and (2.1), the linear equations for ($u', p'$) are thus:

286

$$\frac{\partial u'}{\partial t} + (U \cdot \nabla) u' + (u' \cdot \nabla) U = -\nabla p' + \nabla^2 u', \qquad (2.9)$$

287

$$\nabla \cdot u' = 0. \qquad (2.10)$$

288

289

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

290

$$q'(x, r, \theta, t) = \sum_{m=-\infty}^{+\infty} \hat{q}'_m(x, r, t) e^{im\theta}, \qquad (2.11)$$

291

292

293

294

295

296

297

where the integer $m$ is the azimuthal wavenumber. Since equations (2.9) and (2.10) are linear, the LSA is thus decoupled for each $m$. The linearity of equations (2.9) and (2.10) decouples the Fourier modes, so that the LSA (Linear Stability Analysis) equations can be solved individually 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 2004).

298

299

300

The linear stability problem is addressed using a timesteper approach. From equations (2.9) and (2.10), 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 et al. 2008):

301

$$q'(t + \tau) = \mathcal{A}(\tau) q'(t).$$

302

303

The complex eigenvalue $\zeta_k$ of $\mathcal{A}(\tau)$, associated with the eigenvector $q'_k$, is linked to the eigenvalue $(\sigma_k + i\omega_k)$ of the linear problem (2.9)-(2.10) through an exponential mapping:

304

$$\zeta_k = e^{(\sigma_k + i\omega_k)\tau},$$

from which the growth rate $\sigma_k$ and frequency $\omega_k$ can be retrieved:

305

306

307

308

309

310

$$\sigma_k = \frac{\ln(|\zeta_k|)}{\tau}, \quad \omega_k = \frac{\varphi_k}{\tau}, \qquad (2.12a, b)$$

where $\varphi_k$ is the complex argument of $\zeta_k$. A given base flow ($U, P$) is said to be unstable if, for any $m$, there is at least one perturbation ($u', p'$) for which $\sigma > 0$ (or, equivalently, $|\zeta| > 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 $Gr_{ac}$ for which $\sigma = 0$ is reached is called the critical acoustic Grashof number $Gr_{ac}^c$ and indicates the onset of an instability.
<figure><img src="image_3.png" /><figcaption>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_{beam}$ (equation (2.5)), is shown in green. The displayed mesh corresponds to the $(N, C_R) = (0.25, 4)$ case.</figcaption></figure>

# 3. Computational methodology

## 3.1. *Numerical methods*

Our work involves four different numerical calculations. First, we calculate $F_{ac}$ at the grid points in the fluid domain. Second, we compute the steady axisymmetric 2D-2C streaming flow driven by $F_{ac}$ by time-stepping equations (2.2)-(2.1) until a steady-state is reached. Third, we perform the LSA of the axisymmetric 2D-2C flow for each case listed in table 1. The leading eigenmode and its eigenvalue are computed by solving the eigenvalue problem for $0 \le m \le 10$. This allows to determine $Gr_{ac}^c$ and identify the shape of the unstable perturbation for each case. Fourth, nonlinear 3D-3C unsteady simulations are run for $Gr_{ac} > Gr_{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 $Gr_{ac}^c$ with a Stuart-Landau equation (Landau & Lifschitz 1987), the parameters of which being evaluated from time series of the flow variables (Sheard et al. 2004).

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

Both the steady axisymmetric 2D-2C and the unsteady 3D-3C flows are computed using the open-source spectral-element code Semtex (Blackburn et al. 2019). 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 2004). Convergence of the solution is then achieved by increasing the polynomial degree $N_p$ of the expansion basis. For the 3D-3C unsteady simulations, the azimuthal direction is discretised using a Fourier spectral method (Blackburn & Sherwin 2004).

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 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 \le r \le 0.5$. We also adapted the total number of elements to keep a similar spatial resolution of $F_{ac}$ in the very-near acoustic field ($x \le 1$) for each case (table 2).

Equations (2.1)-(2.2) are discretised in time using the operator splitting scheme of Karniadakis et al. (1991). The scheme integrates explicitly the inertia terms of equation (2.2), and implicitly its viscous terms. For the unsteady 3D-3C simulations, the temporal integration scheme is used in its third-order formulation. For the 2D-2C base flow calculations, as only steady solutions of equations (2.1)-(2.2) are sought, the scheme is used in its first-order
<table><caption>Table 2: Characteristics of the meshes used for each case. See table 1 for a complete list of the computational parameters defining each case.</caption><thead><tr><th>N</th><th>C<sub>R</sub></th><th>Elements in axial direction</th><th>Elements in radial direction</th><th>Total number of elements</th></tr></thead><tbody><tr><td>0.25</td><td>1</td><td>24</td><td>5</td><td>120</td></tr><tr><td>0.25</td><td>2</td><td>24</td><td>7</td><td>168</td></tr><tr><td>0.25</td><td>3</td><td>24</td><td>9</td><td>216</td></tr><tr><td>0.25</td><td>4</td><td>24</td><td>9</td><td>216</td></tr><tr><td>0.25</td><td>6</td><td>24</td><td>10</td><td>240</td></tr><tr><td>1</td><td>6</td><td>34</td><td>14</td><td>476</td></tr></tbody></table>

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. (2008)), 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.9)-(2.10) using the second-order formulation of the operator splitting scheme of Karniadakis et al. (1991). The eigenvalues and associated eigenvectors of the stability problem are then recovered by applying the iterative algorithm of Barkley et al. (2008).

Finally, for a given set of values of $N$, $C_R$ and $Gr_{ac}$, the initial condition for the 2D-2C flow calculations is either set to $u = 0$, or to a previous solution obtained for the same ($N, C_R$) but for a different $Gr_{ac}$ to accelerate the convergence towards a steady state. The unsteady 3D-3C simulations of the $Gr_{ac} > Gr_{ac}^c$ regimes are initialised with the steady axisymmetric base flow duplicated along $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 3D-3C 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 \{3, 4\}$, and to 0.01 for $C_R \in \{1, 2\}$.

## 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 $\tilde{I}_x(x, r = 0)$, which we compared to the numerical evaluation of equation (2.4) using the methodology outlined in § 3.1 (Vincent et al. 2024a). 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$.

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. 2024a), the very near acoustic field (i.e., $x \ll 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
<figure><img src="image_4.png" /><figcaption>Figure 3: On-axis normalised acoustic intensity profile $\tilde{I}(x, r = 0)$ for $(N, C_R) = (0.25, 4)$. The analytical profile, proposed by Vincent et al. (2024a), 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/(4S)$ (purple dashed vertical line). The values obtained by numerically evaluating equations (2.4) 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.</figcaption></figure>

<table><caption>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 $Gr_{ac} = 7000$ and $m = 2$. The error is relative to the value of $\sigma_{max}$ obtained for $N_p = 14$.</caption><thead><tr><th>N<sub>p</sub></th><th>σ<sub>max</sub></th><th>Relative difference (%)</th></tr></thead><tbody><tr><td>5</td><td>0.10249</td><td>33.11</td></tr><tr><td>6</td><td>0.07379</td><td>4.17</td></tr><tr><td>7</td><td>0.07792</td><td>1.20</td></tr><tr><td>8</td><td>0.07695</td><td>0.05</td></tr><tr><td>9</td><td>0.07710</td><td>0.13</td></tr><tr><td>10</td><td>0.07709</td><td>0.12</td></tr><tr><td>11</td><td>0.07696</td><td>0.05</td></tr><tr><td>12</td><td>0.07702</td><td>0.03</td></tr><tr><td>13</td><td>0.07701</td><td>0.02</td></tr><tr><td>14</td><td>0.07700</td><td>—</td></tr></tbody></table>

leading eigenvalue remained below 0.5 %. An example of grid sensitivity analysis for the $(N, C_R) = (0.25, 4)$ case at $Gr_{ac} = 7000$ (the largest $Gr_{ac}$ for that setup) 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 $\approx 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 appendix A).

For the 3D-3C unsteady simulations, $u$ and $p$ were expanded along $e_\theta$ using $N_F$ Fourier
<figure><img src="image_5.png" /><figcaption>Figure 4: Map of the steady velocity magnitude $||U||$ computed for ($N, C_R$) = (0.25, 6) at $Gr_{ac}$ = 6.4 × 10³. 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 ($U = 0$, purple discs). The green line depicts the approximate beam radius given by equation (2.5). The base flow is axisymmetric and is reflected about the $r = 0$ axis for the sake of visualisation.</figcaption></figure>

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 $Gr_{ac}$ = 15500 for instance, the relative difference between the kinetic energies computed with either $N_F$ = 32 or 64 was ≈ 0.03 %. Even smaller differences were obtained for the $Gr_{ac} > Gr_{ac}^c$ simulations of the $C_R \in \{4, 6\}$ cases. The 3D-3C 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 $Gr_{ac}$ = 7000 and with $N_p$ = 8 was time-stepped with $\Delta t$ = 10⁻⁴, yielding a maximum CFL number of 0.476.

## 4. Two-dimensional axisymmetric base flow

### 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 illustrative example. The main flow features are identified by the patterns of streamlines together with critical point theory (Chong et al. 1990). The critical point analysis is commonly used to identify structures within complex 3D flows (Hunt et al. 1978). 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 3D DNS (Direct Numerical Simulations) of separated flows around bluff bodies (Dousset & Pothérat 2010, 2012) or complex conduits (Pothérat, A. and Zhang, L. 2018) offered a powerful way to reliably extract the main structures governing the flow's dynamics.
The main idea is to extract 2D 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. 1978). Second, these points offers 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.

The critical points where $U = 0$ are shown in purple (labelled $S_i$ for saddles and $N_i$ for nodes, respectively) 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 $\Sigma_N$, of half-nodes $\Sigma_{N'}$, of saddles $\Sigma_S$ and half-saddles $\Sigma_{S'}$ in a slice of $U$ are related to each other through (Hunt et al. 1978; Foss 2004):

$$2\Sigma_N + \Sigma_{N'} - 2\Sigma_S - \Sigma_{S'} = \chi, \qquad (4.1)$$

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 2004)). Inspection of the streamlines in figure 4 gives:

$$\Sigma_N = 8, \quad \Sigma_{N'} = 0, \quad \Sigma_S = 0, \quad \Sigma_{S'} = 14,$$

so that the left-hand side of equation (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 (1964) 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 equation (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 2D-2C 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 $Gr_{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. Asides of 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 $Gr_{ac}^c$.

Besides 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 \le C_R \le 6$ and $Gr_{ac} = 6400$. In a weakly-confined setting, the jet strongly
<figure><img src="image_6.png" /><figcaption>Figure 5: Map of the steady velocity magnitude $||U||$ computed for ($N = 0.25, C_R = 4$) for $Gr_{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 ($U = 0$, purple discs) and $S_i$ are saddles ($U = 0$, purple squares). The green line depicts the approximate beam radius given by equation (2.5). The base flow is axisymmetric and is reflected about the $r = 0$ axis for the sake of visualisation.</figcaption></figure>

accelerates over $0 \le x \le 20$ for all $C_R$ as a result of a balance between $F_{ac}$ and the inertia forces (Moudjed et al. 2014a; Vincent et al. 2024a). Decreasing $C_R$ enhances the jet deceleration at larger $x$. For $2 \le C_R \le 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 $F_{ac}$ and viscous forces (Vincent et al. 2024a). 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 \ge 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 1998).

## 5. Stability analysis

### 5.1. Non-oscillatory instabilities in weakly-confined settings

We now turn to the linear stability analysis of the steady 2D-2C flows similar to those shown in figures 4 and 5. The variations of the growth rate of the leading eigenmode $\sigma$ with $Gr_{ac}$ and $m$ for $N \in \{0.25, 1\}$ and $4 \le C_R \le 6$ are shown in figure 7. For $4 \le C_R \le 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 \le C_R \le 6$, but $Gr_{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 2023).

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 (figure 8 (a) and (b)), the perturbation intensity as measured by $u'_x$ is strongest in a region stretching from the impingement region of the jet to the shear layer between recirculations
<figure><img src="image_7.png" /><figcaption>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 $Gr_{ac} = 6400$. The profiles are obtained for $N = 0.25$ and $1 \le C_R \le 6$. The dashed vertical line corresponds to the Fresnel distance $L_F$ marking the transition between the near ($x < L_F$) and the far ($x > L_F$) acoustic fields. Inset: variations with $C_R$ of the streamwise pressure gradient $\langle dP/dx \rangle$ of the base flow at the mid cavity length $L/2$ and averaged over the cross-sectional area.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: Growth rate $\sigma$ of the leading eigenmode as a function of the azimuthal wavenumber $m$ for $3 \le C_R \le 6$. Oscillatory and non-oscillatory modes are represented by filled circles and triangles, respectively. The insets show the evolution of $\sigma$ with $Gr_{ac}$ for the indicated value of $m$. The critical Grashof number $Gr_{ac}^c$, obtained through spline interpolation of $\sigma$, is reported in red. The filled black symbols in the inset of (b) are additional eigenvalue computations made to improve the interpolation accuracy near $\sigma = 0$.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8: Leading mode for $(N, C_R) = (0.25, 4)$ at $Gr_{ac} = 6400$. The mode is non-oscillatory ($\omega = 0$) and unstable ($\sigma = 1.59 \times 10^{-2} > 0$). (a) Axial velocity perturbation $u'_x$ in the $(x, r)$ plane, along with the streamlines and critical points of the base velocity $U$. (b) Velocity perturbation $u'$, with the background color corresponding to its azimuthal component $u'_\theta$. (c) to (h): Slices of $u'$ in constant-$x$ planes located by the blue vertical lines in (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.5), 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.</figcaption></figure>

<footnote n="486">486 centred on $N_4$ and $N_5$. Radially, the perturbation is located within the acoustic beam (the approximate acoustic beam radius $R_{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).</footnote>

<footnote n="489">489 The perturbation then crosses the shear layer at the boundary of the jet to reach the region</footnote>
between the two recirculations around $N_4$ and $N_5$. Hence, the mechanism is not associated to a shear layer instability but rather coincides precisely with the strong velocity gradients in the impingement region centred on $S'_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 to 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.

To best capture the differences between the modes associated to 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:

$$ \Sigma_N = 8, \quad \Sigma_{N'} = 4, \quad \Sigma_S = 9, \quad \Sigma_{S'} = 0. $$

Since $\chi = 2$ for the downstream wall (Foss 2004), 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 2001). 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 outer region of the Moffatt vortices of the base flow and the main 4-vortices structures of the perturbation. Decreasing $C_R$ narrows down the gap between the two separators, as seen by comparing cases with $C_R = 6$ and $C_R = 4$, respectively shown in figures 9 (a) and (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 azimuthal direction.

For $C_R = 3$, another local maximum of $\sigma(m)$ is first to become unstable at $Gr_{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 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 $\Sigma_S$ is multiplied in the topological rule (4.1)), this type of critical point
$(N, C_R) = (0.25, 6): \text{Gr}_{ac} = 4750, m = 2,$
$\sigma = 2.27 \times 10^{-3}, \omega = 0 \text{ (Branch I)}$

<figure><img src="image_10.png" /></figure>

$(N, C_R) = (0.25, 4): \text{Gr}_{ac} = 6400, m = 2,$
$\sigma = 1.59 \times 10^{-2}, \omega = 0 \text{ (Branch I)}$

<figure><img src="image_11.png" /></figure>

$(N, C_R) = (0.25, 3): \text{Gr}_{ac} = 8000, m = 4,$
$\sigma = 7.25 \times 10^{-2}, \omega = 0 \text{ (Branch II)}$

<figure><img src="image_12.png" /><figcaption>Figure 9: Lines of the skin friction perturbation stresses on the downstream wall, with the background color 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 for $N = 0.25$ and for different $C_R$: (a) $C_R = 6$, (b) 4 and (c) 3. The approximate beam radius (equation (2.5)) 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 $Gr_{ac}$ slightly above the instability onset.</figcaption></figure>

538 must have its own topological weight. Assuming that the topological rule (4.1) is satisfied,

539 the topological weight of this 8-branch star saddle must be $w_{8s} = -6$.

540

5.2. *Rise of oscillatory instabilities in confined settings*

541 As $C_R$ is further decreased, the modes destabilising the flow become oscillatory. Indeed,

542 the $C_R = 2$ flow is destabilised at $Gr_{ac}^c = 10513$ by the growth of an oscillatory $m = 2$

543 perturbation (figure 10 (a)). This mode is associated with a new branch that we label as
<figure><img src="image_13.png" /><figcaption>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 $Gr_{ac}$ for the indicated value of $m$. The critical Grashof number $Gr_{ac}^c$, obtained through spline interpolation of $\sigma$, is reported in red. The filled black symbols in the insets of (a) are additional eigenvalue computations made to improve the interpolation accuracy in the vicinity of $\sigma = 0$.</figcaption></figure>

“Branch III”. Interestingly, Branch III and Branch II coexist for this $C_R$, and the leading Branch II mode becomes unstable at $Gr_{ac} \approx 10630 \approx 1.01Gr_{ac}^c$. The $\sigma$ of the Branch II mode increases more rapidly with $Gr_{ac}$ than the $\sigma$ of the Branch III mode, so that Branch II becomes dominant for $Gr_{ac} \ge 11000 \approx 1.04Gr_{ac}^c$.

For $C_R = 2$, the increased flow confinement and the greater $Gr_{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 to the previous $C_R \ge 3$ cases. First, $|u'_\theta|$ 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 $u'$ and the contours of $\omega_x$ at $N_5$ and $S_1$ (figure 12): as for $C_R \ge 3$, the perturbation is still clustered near the shear layer around the jet where $U_x \approx 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. (2022), who reported similar unstable waves for a plane 2D 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 the 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$.

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 $Gr_{ac} = 15310$. At $Gr_{ac} \approx Gr_{ac}^c$, both the enhanced confinement and the greater strength of the return flow resulting from the larger $Gr_{ac}^c$ significantly affect the base flow topology near the downstream wall (figure 14). First, the size of the recirculating flow structures is greatly
$(N, C_R) = (0.25, 2): \text{Gr}_{ac} = 10600, m = 2, \sigma = 2.54 \times 10^{-2}, \omega = 3.9462 \text{ (Branch III)}$

<figure><img src="image_14.png" /><figcaption>Figure 11: Leading mode for $(N, C_R) = (0.25, 2)$ at $Gr_{ac} = 10600$. The mode is oscillatory ($\omega \neq 0$) and unstable ($\sigma = 2.54 \times 10^{-2} > 0$). (a) Azimuthal velocity perturbation $u'_\theta$ in the $(x, r)$ plane, along with the critical points of the base velocity $U$. (b) Detailed view near the impingement (region framed in red in (a)) of the axial velocity perturbation $u'_x$ together with the base flow streamlines. (c) Details $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.5), and the blue solid line locates the radius where $U_x$ is 50% of its on-axis value.</figcaption></figure>

<footnote n="574">574 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 $U$. These new critical points appear at approximately the same radial location as $N_5$ and $S_1$.</footnote>

<footnote n="577">577 Maps of $u'$ in the $(x, r)$ plane show that $|u'_x|$ and $|u'_\theta|$ are maximum near $N_5$ and $S_1$, i.e., in regions of high shear in the bulk of the base flow (figure 14 (b) and (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</footnote>
<figure><img src="image_15.png" /></figure>

Figure 12: Leading velocity perturbation $\mathbf{u}'$ for ($N, C_R$) = (0.25, 2) at $Gr_{ac}$ = 10600. 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'(\text{(a)-(b)})$ or the streamwise vorticity perturbation $\omega_x'(\text{(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.5), 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.

581

582

583

584

of an adverse pressure gradient slowing down the jet (inset of figure 6). Besides, $\mathbf{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).

585

586

587

588

589

To conclude, for $N \in \{0.25, 1\}$ and $3 \le C_R \le 6$, the primary instabilities are caused by non-oscillatory perturbations originating from the jet impingement. Further confining the flow ($1 \le C_R \le 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.

590

# 6. Bifurcation characterisation

591

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

592

593

594

595

596

597

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 (1987) 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. 2004), the wake behind a sphere (Thompson et al. 2001) and behind one or multiple cylinders (Henderson & Barkley 1996; Henderson 1997;
$(N, C_R) = (0.25, 2): Gr_{ac} = 10600, m = 2,$
$\sigma = 2.54 \times 10^{-2}, \omega = 3.9462$ (Branch III)

$(N, C_R) = (0.25, 1): Gr_{ac} = 15500, m = 3,$
$\sigma = 1.7 \times 10^{-1}, \omega = 21.472$ (Branch IV)

<figure><img src="image_16.png" /><figcaption>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 equation (2.5)). The base flow half saddle $S'_7$ is also reported in red.</figcaption></figure>

598 Carmo et al. 2008). The Stuart-Landau model reads:

$$ \frac{dA}{dt} = (\sigma + i\omega) A - l(1 + ic) |A|^2 A + O(A^5). \qquad (6.1) $$

600 where A is a complex perturbation, and $l$ and $c$ are two real constants. The mechanism 601 responsible for the saturation of A can be easily understood from the sign of $l$ in equation (6.1).

602 For $l > 0$, the cubic term in the right-hand side of equation (6.1) acts against the growth of 603 the perturbation: this situation corresponds to a supercritical bifurcation. On the contrary, 604 for $l < 0$, higher-order terms are required for A to saturate: the bifurcation is thus subcritical. 605

606 Therefore, only the sign of $l$ is needed to determine the nature (or criticality) of the bifurcation.

607 As explained in Sheard et al. (2004) and Kumar & Pothérat (2020), $l$ may be inferred from a 608 more convenient form of equation (6.1). Plugging $A(t) = |A(t)|e^{i\phi(t)}$ into the Stuart-Landau 609 model (6.1) gives, after separating the real and imaginary part:

$$ \frac{d \log |A|}{dt} = \sigma - l|A|^2 + O(|A|^4), \qquad (6.2) $$

$$ \frac{d\phi}{dt} = \omega - lc|A|^2 + O(|A|^4). \qquad (6.3) $$

611 The nature of the bifurcation is determined by plotting $(d \log |A|/dt)$ against $|A|^2$; the slope 612 of the resulting curve is $-l$ for sufficiently small $|A|^2$, and $\sigma$ is obtained from the y-intercept.

613 From this we compare $\sigma$ to its value obtained from the LSA results and so assess the 614 accuracy of the estimate for $l$. Note that for an non-oscillatory perturbation, $\phi(t) = 0$, so 615 only equation (6.2) is needed.

616 The saturated state is characterised using the steady-state forms of equations (6.2) and 617 (6.3). For a supercritical ($l > 0$) bifurcation for instance, the first two terms in the right- 618 hand side of the amplitude equation (6.2) are dominant, so that the saturated perturbation
$(N, C_R) = (0.25, 1): Gr_{ac} = 15500, m = 3, \sigma = 1.701 \times 10^{-1}, \omega = 21.472$ (Branch IV)

<figure><img src="image_17.png" /><figcaption>Figure 14: Leading mode for $(N, C_R) = (0.25, 1)$ at $Gr_{ac} = 15500$. The mode is oscillatory ($\omega \neq 0$) and unstable ($\sigma = 1.701 \times 10^{-1} > 0$). (a) Azimuthal velocity perturbation $u'_\theta$ in the $(x, r)$ plane, along with the critical points of the base velocity $U$. (b) Detailed view near the impingement (region framed in red in (a)) of the axial velocity perturbation $u'_\chi$ together with the base flow streamlines. (c) Details of $u'$ near the impingement, with $u'_\theta$ taken as background colour.</figcaption></figure>

619 amplitude is:

$$ |A_{sat}| = \sqrt{\frac{\sigma}{l}}. \qquad (6.4) $$

620 Furthermore, if the saturated state is a time-periodic signal of constant amplitude $|A_{sat}|$, then $d\phi/dt$ reduces to the saturated frequency $\omega_{sat}$. Using equation (6.4), the steady-state version of equation (6.3) gives, for the Landau constant $c$:

$$ c = \frac{\omega - \omega_{sat}}{\sigma}, \qquad (6.5) $$

625 thus allowing for all the quantities in the truncated Stuart-Landau model (6.1) to be fully characterised.
<figure><img src="image_18.png" /><figcaption>Figure 15: Leading velocity perturbation $u'$ for ($N, C_R$) = (0.25, 1) at $Gr_{ac}$ = 15500. The mode is shown in constant-$x$ planes located by the blue vertical lines in figure 14, with either $u'_x$ ((a)-(b)) or the streamwise vorticity perturbation $\omega'_x$ ((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.5), 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.</figcaption></figure>

We shall now define $A$ in terms of flow quantities. We followed Kumar & Pothérat (2020) and based $|A|$ 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. 2001; Sheard et al. 2004; Sherwin & Blackburn 2005; Sapardi et al. 2017). 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 $|u'_\theta|$. 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.

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

$$r_c = \frac{Gr_{ac}}{Gr_{ac}^c} - 1, \qquad (6.6)$$

so that $r_c > 0$ corresponds to a $Gr_{ac} > Gr_{ac}^c$ flow regime, and $r_c < 0$ to a $Gr_{ac} < Gr_{ac}$ regime. Finally, $(d \log |A|/dt)$ was computed from the recorded time series using centred finite differences.
<table><caption>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 3D-3C simulations.</caption><thead><tr><th>N</th><th>C<sub>R</sub></th><th>(x, r, θ)</th></tr></thead><tbody><tr><td>0.25</td><td>1</td><td>(82, 1.8, 0)</td></tr><tr><td>0.25</td><td>2</td><td>(81, 1.3, 0)</td></tr><tr><td>0.25</td><td>3</td><td>(80, 2, 0)</td></tr><tr><td>0.25</td><td>4</td><td>(78, 6, 0)</td></tr><tr><td>0.25</td><td>6</td><td>(75, 2, 0)</td></tr><tr><td>1</td><td>6</td><td>(308, 30, 0)</td></tr></tbody></table>

## 6.2. Topology of the perturbed states from the nonlinear 3D-3C 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}$ for four representative cases, one per instability branch. These were obtained slightly above onset, by subtracting the base axisymmetric solution from the full nonlinear solution. The results are illustrated in figure 16. Comparing to the topology of the leading eigenmodes in figures 8, 9, 11, 12, 13, 14, 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 3D-3C nonlinear simulations on figure 21. In all cases, the leading eigenmode comes very close to the full nonlinear solution. In particular, the nonlinear 3D-3C 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.

## 6.3. Nature of the bifurcations

The Stuart-Landau analysis has been carried for all the setups listed in table 1; the cases defined by $N = 0.25$ and $C_R \in \{1, 2, 4\}$ are shown in figure 17 as examples. For each case, the initial instants of $|A|$ are dominated by noise (figure 17 (a), (c) and (e)). The Stuart-Landau equation (6.2) for $|A|$ 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.

For $4 \le C_R \le 6$, $l > 0$ near the instability onset. The leading Branch I mode thus destabilises the flow through a supercritical circular pitchfork bifurcation (Touihri et al. 1999). Further confining the flow yields $l < 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 > 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 3D-3C 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 \{2, 3\}$ raises the question of whether the corresponding instabilities can be triggered for $r_c < 0$. For both cases, introducing white noise for $-0.05 \le r_c \le -0.01$ did not cause the flow to deviate from its axisymmetric base state. Triggering these bifurcations at $Gr_{ac} < Gr_{ac}^c$ would require different tools, e.g. the analysis of the transient growth of non-modal perturbations (Schmid 2007). 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.
<figure><img src="image_19.png" /><figcaption>Branch IV</figcaption></figure>

$Gr_{ac} = 15500 (r_c = 0.012)$

$t = 130$

<figure><img src="image_20.png" /><figcaption>Branch III</figcaption></figure>

$Gr_{ac} = 10600 (r_c = 0.008)$

$t = 189$

<figure><img src="image_21.png" /><figcaption>Branch II</figcaption></figure>

$Gr_{ac} = 8000 (r_c = 0.032)$

$t = 90$

<figure><img src="image_22.png" /><figcaption>Branch I</figcaption></figure>

$Gr_{ac} = 2500 (r_c = 0.141)$

$t = 2000$

<table><caption>Figure 16: Nonlinear 3D-3C unsteady simulations of the velocity perturbation represented by contours of azimuthal velocity perturbation $u'_\theta$ set at $\pm 0.2 \max|u'_\theta|$ for the leading eigenmode corresponding to unstable branches IV (a), III (b), II (b) and I (d). The black surface represents the position of the approximate beam radius given by equation (2.5).</caption><thead><tr><th>N</th><th>C<sub>R</sub></th><th>Gr<sub>ac</sub></th><th>r<sub>c</sub></th><th>σ<br>(LSA)</th><th>σ<br>(nonlinear<br>3D-3C)</th><th>ε<sub>σ</sub><br>(%)</th><th>ω<br>(LSA)</th><th>ω<br>(nonlinear<br>3D-3C)</th><th>ε<sub>ω</sub><br>(%)</th><th>ω - ω<sub>sat</sub></th></tr></thead><tbody><tr><td>0.25</td><td>1</td><td>15500</td><td>0.0124</td><td>1.701 × 10<sup>-1</sup></td><td>1.703 × 10<sup>-1</sup></td><td>0.12</td><td>21.472</td><td>21.409</td><td>0.29</td><td>0</td></tr><tr><td>0.25</td><td>2</td><td>10600</td><td>0.0083</td><td>2.540 × 10<sup>-2</sup></td><td>2.601 × 10<sup>-2</sup></td><td>2.35</td><td>3.946</td><td>3.958</td><td>0.30</td><td>1.166</td></tr><tr><td>0.25</td><td>3</td><td>8000</td><td>0.0320</td><td>7.245 × 10<sup>-2</sup></td><td>7.690 × 10<sup>-2</sup></td><td>5.79</td><td>0</td><td>0</td><td>—</td><td>—</td></tr><tr><td>0.25</td><td>4</td><td>6400</td><td>0.0262</td><td>1.591 × 10<sup>-2</sup></td><td>1.594 × 10<sup>-2</sup></td><td>0.19</td><td>0</td><td>0</td><td>—</td><td>—</td></tr><tr><td>0.25</td><td>6</td><td>5000</td><td>0.0635</td><td>1.384 × 10<sup>-2</sup></td><td>1.383 × 10<sup>-2</sup></td><td>0.07</td><td>0</td><td>0</td><td>—</td><td>—</td></tr><tr><td>1</td><td>6</td><td>2500</td><td>0.1411</td><td>2.755 × 10<sup>-3</sup></td><td>2.896 × 10<sup>-3</sup></td><td>4.87</td><td>0</td><td>0</td><td>—</td><td>—</td></tr></tbody></table>

<caption>Table 5: Comparison between the growth rates $\sigma$ and frequencies $\omega$ returned by linear stability analysis (LSA) with those obtained from nonlinear 3D-3C unsteady simulations of $Gr_{ac} > Gr_{ac}^c$ regimes ($r_c > 0$, see equation (6.6)). The relative errors between the values of $\sigma$ and $\omega$ obtained by LSA and the nonlinear 3D-3C simulations are given by $\epsilon_\sigma$ and $\epsilon_\omega$, respectively.</caption>
<figure><img src="image_23.png" /></figure>

Figure 17: Stuart-Landau analysis for (a)-(b) ($N = 0.25, C_R = 4$) at $Gr_{ac} = 6400$ ($r_c = 0.0262$), (c)-(d) ($N = 0.25, C_R = 2$) at $Gr_{ac} = 10600$ ($r_c = 0.0083$) and (e)-(f) ($N = 0.25, C_R = 1$) at $Gr_{ac} = 15500$ ($r_c = 0.0124$). Left: 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$. Right: evolution of $d(\log|A|)/dt$ with $|A|^2$, where $|A|$ is the amplitude of the perturbation plotted on the left figures (for the $C_R = 2$ and 1 cases, $|A|$ is based on the envelope of $u'_\theta$). The growth rate $\sigma$ and the Landau coefficient $l$ are obtained by fitting equation (6.2) to the red curve. For each case, the fitting range is highlighted in red.

## 6.4. Saturated states and secondary instabilities

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

For Hopf bifurcations, the frequency $\omega_{sat}$ of the saturated perturbation may differ from $\omega$ in the linear regime (regime of exponential perturbation growth), yielding a non-zero $c$ (equation (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).
<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /><figcaption>Figure 18: Frequencies obtained from the time-series of the azimuthal velocity perturbation $u'_\theta(t)$ for (a) $(N, C_R) = (0.25, 1)$ and (b) $(0.25, 2)$. The corresponding time-series are those displayed in figure 17 (e) and (c), respectively. The frequencies are computed from FFTs made in the linear regime ($\omega$, blue) and in the saturated regime ($\omega_{sat}$, red).</figcaption></figure>

For $(N, C_R) = (0.25, 2)$, $\omega_{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 3D-3C simulations is 0.30 % (table 5). For $(N, C_R) = (0.25, 1)$, the $\omega$ and $\omega_{sat}$ obtained with the nonlinear 3D-3C 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 of 0.29 % (table 5) and 0.87 % in the linear and saturated regimes, respectively. Given these small differences, $\omega$ and $\omega_{sat}$ are considered to be equal, hence $c = 0$.

For $(N, C_R) = (0.25, 2)$, the large uncertainty on $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. Contrarily to the $(N, C_R) = (0.25, 2)$ case, the slow deviation from the saturated state allows
<table><caption>Table 6: Summary of the Stuart-Landau analysis based on local measurements of $u'_{\theta}$; these measurements are obtained from nonlinear 3D-3C unsteady simulations of $Gr_{ac} > Gr_{ac}^c$ regimes ($r_c > 0$, see equation (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_{sat}$ are determined. See table 4 for the coordinates of the points where $u'_{\theta}$ is recorded.</caption><thead><tr><th>N</th><th>C<sub>R</sub></th><th>Gr<sub>ac</sub></th><th>r<sub>c</sub></th><th>l</th><th>c</th><th>Branch</th><th>Bifurcation type</th></tr></thead><tbody><tr><td>0.25</td><td>1</td><td>15500</td><td>0.0124</td><td>(2.856 ± 0.012) × 10<sup>-2</sup></td><td>0</td><td>IV</td><td>Supercritical Hopf</td></tr><tr><td>0.25</td><td>2</td><td>10600</td><td>0.0083</td><td>-1.547 ± 0.001</td><td>44.83 ± 14.63</td><td>III</td><td>Subcritical Hopf</td></tr><tr><td>0.25</td><td>3</td><td>8000</td><td>0.0320</td><td>-0.794 ± 0.003</td><td>0</td><td>II</td><td>Subcritical circular pitchfork</td></tr><tr><td>0.25</td><td>4</td><td>6400</td><td>0.0262</td><td>0.391 ± 0.004</td><td>0</td><td>I</td><td>Supercritical circular pitchfork</td></tr><tr><td>0.25</td><td>6</td><td>5000</td><td>0.0635</td><td>0.142 ± 0.001</td><td>0</td><td>I</td><td>Supercritical circular pitchfork</td></tr><tr><td>1</td><td>6</td><td>2500</td><td>0.1411</td><td>(1.971 ± 0.006) × 10<sup>-2</sup></td><td>0</td><td>I</td><td>Supercritical circular pitchfork</td></tr></tbody></table>

<figure><img src="image_26.png" /><figcaption>Figure 19: Time-series of the azimuthal velocity perturbation $u'_{\theta}$ for ($N, C_R$) = (0.25, 2) and $Gr_{ac} = 10600$ ($r_c = 0.0083$). The signal is recorded at ($x, r, \theta$) = (81, 1.3, 0).</figcaption></figure>

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 \ge 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.

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 $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 1993). Although these complex amplitudes cannot be obtained from LSA, they do appear in the nonlinear 3D-3C unsteady simulations. The leading LSA mode shapes are compared to snapshots of the 3D-3C 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
<figure><img src="image_27.png" /></figure>

Figure 20: Stuart-Landau analysis for (N, C<sub>R</sub>) = (0.25, 3) at Gr<sub>ac</sub> = 8000 (r<sub>c</sub> = 0.0320). (a) Temporal evolution of the azimuthal velocity perturbation u'<sub>θ</sub> obtained from nonlinear 3D-3C unsteady simulations and recorded at (x, r, θ) = (80, 2, 0). The perturbation originates from white noise injected at t = 0. (b) Evolution of d (log |A|) / dt with |A|^2 (red curve), where |A| = |u'<sub>θ</sub>| and is extracted from (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.

716 perturbation to avoid nonlinear effects on the perturbation shape. For both cases, a strong
717 agreement is found between the Branches III and IV mode shapes obtained by LSA and
718 by the nonlinear 3D-3C simulations, and both perturbations correspond to waves travelling
719 along *e*<sub>θ</sub> (see the supplementary movies 1 and 2).

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

724 **7. Conclusion**

725 We studied the stability of a laminar Eckart streaming jet flowing in a closed cylindrical
726 cavity. A plane circular transducer placed at one end of the cavity radiates an axisymmetric
727 diffracting sound beam, whose attenuation drives a steady axisymmetric jet impinging the
728 cavity wall facing the transducer. To our knowledge, this work is the first stability analysis of
729 a streaming jet forced by a realistic sound field that includes the effects of both diffraction
730 and attenuation on its spatial structure. As such, this study provides a generic framework to
731 understand the conditions in which realistic acoustic streaming jets are stable, and their path
732 to unsteadiness when they are not.

733 We assessed the effect of the cavity size on the streaming flow stability by either varying
734 the cavity length with respect to the attenuation length (N = 1 and 0.25), and by changing
735 the ratio C<sub>R</sub> between the cavity radius and the beam radius at the downstream cavity wall
736 (C<sub>R</sub> ∈ [1, 6]). For each case, we (a) identified the destabilisation mechanism by analysis of
737 the flow topology and its critical points, (b) computed the critical forcing magnitude Gr<sub>ac</sub><sup>c</sup>
738 above which the flow is linearly unstable (Linear Stability Analysis), and (c) determined the
739 nature of the bifurcation by means of Stuart-Landau analyses of data obtained from nonlinear
740 3D-3C (three dimensions, three components) unsteady simulations. The main outcomes of
741 this work are summarised below:

742 (i) For all cavity sizes, recirculating structures originating from the jet impingement are
743 formed near the downstream wall. The primary bifurcations result from the destabilisation
744 of these structures, rather than from the destabilisation of the jet itself.
$$ (a) \qquad (N, C_R) = (0.25, 1), \quad Gr_{ac} = 15500, \\ \qquad \sigma = 1.701 \times 10^{-1}, \quad \omega = 21.472 \text{ (Branch IV)} $$

<figure><img src="image_28.png" /></figure>

$$ (b) \qquad (N, C_R) = (0.25, 2), \quad Gr_{ac} = 10600, \\ \qquad \sigma = 2.54 \times 10^{-2}, \quad \omega = 3.9462 \text{ (Branch III)} $$

<figure><img src="image_29.png" /><figcaption>Figure 21: Comparison between the leading velocity perturbation $u'$ obtained with LSA (left) and the perturbation obtained with nonlinear 3D-3C 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 3D-3C unsteady simulations are compared at the location of the axisymmetric base flow saddle $S_1$ ($x = 80.15$ in (a) and $x = 77.52$ in (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 $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 equation (2.5). See the supplementary movies (movies 1 and 2) for the animations of the velocity perturbations obtained from the nonlinear 3D-3C simulations.</figcaption></figure>
<figure><img src="image_30.png" /><figcaption>Figure 22: Variations of the critical Grashof number $Gr_{ac}^c$ with the radial confinement ratio $C_R$ defined by equation (2.6). The $N = 0.25$ cases are represented by filled discs and are interpolated by cubic splines. The $Gr_{ac}^c$ obtained for $(N, C_R) = (1, 6)$ is also reported (black square) and illustrates the destabilising effect of increasing $N$. The dashed vertical lines are arbitrarily placed in between cases having different primary bifurcation types.</figcaption></figure>

(ii) Reducing the cavity size has a stabilising effect (figure 22): the critical $Gr_{ac}$ for the onset of instability ranges from $Gr_{ac}^c = 2191$ for $(N, C_R) = (1, 6)$ (less-confined case) to $Gr_{ac}^c = 15310$ for $(N, C_R) = (0.25, 1)$ (most confined case). For $4 \le C_R \le 6$, the leading perturbation is non-oscillatory (Branch I). This unstable mode displays large perturbation velocities near the jet impingement. The first significant effect of confinement occurs for $C_R = 3$: the leading mode remains non-oscillatory, but its topology changes (Branch II). The reason for the change of topology is that the mode corresponding to Branch I with an $m = 2$ azimuthal wavenumber becomes radially so confined that it becomes azimuthally stretched to the point of breaking up. At this point, an $m = 4$ mode from Branch II emerges. Flows in further confined settings are destabilised by oscillatory perturbations travelling in the azimuthal direction (Branch IV for $C_R = 1$). For $C_R = 2$, these perturbations also propagate backwards, i.e., towards the source (Branch III). Overall, reducing $C_R$ (a) increases the shear in a layer separating counter-rotating structures in the bulk of the flow, and (b) enhances an adverse pressure gradient slowing down the jet. As a consequence, reducing $C_R$ moves the locus of instability from the jet impingement to the shear layer between the jet and the recirculating structures.

(iii) Finally, changing the flow confinement significantly affects the nature of the primary bifurcations. Whilst the unstable modes of Branches I and IV respectively trigger supercritical pitchfork and Hopf bifurcations, the bifurcations associated with Branch II and III modes are subcritical.

The subcritical nature of bifurcations for the $C_R \in \{2, 3\}$ cases raises the question of whether these bifurcations could be triggered $Gr_{ac}^c < Gr_{ac}^c$. The nonlinear 3D-3C unsteady simulations carried out between 0.95 $Gr_{ac}^c$ and 0.99 $Gr_{ac}^c$ revealed that white-noise-based perturbations could not ignite subcritical transitions towards another state. Whether this scenario could occur remains uncertain at this stage and shall require further analysis, for instance by finding the initial disturbances maximising the energy growth over a finite time range (Schmid 2007). Other examples of convective of shear flows exhibit similar subcritical bifurcations where transition away from the steady state is not easily triggered in the subcritical regime (Kumar & Pothérat 2020; Camobreco et al. 2020, 2021, 2023; Huang et al. 2024). The stability and transition to turbulence in these flows are very different to the classical shear flows such as Couette or Poiseuille flows, where the subcritical nature of the bifurcation favours transition to turbulence via mechanisms that radically differ from
those of the LSA. Camobreco et al. (2023) and Huang et al. (2024) showed that transition to turbulence in quasi-2D and 2D shear flows could be triggered only in mildly subcritical regimes and that the mode mediating the transition was in fact the leading LSA eigenmode. In acoustically-driven jets, the instability mainly originates from shear, but the base flow is significantly more complex than in these examples. Hence an interesting question for further work is whether when the transition is weakly subcritical, the leading eigenmodes underpin any transition to unsteadiness and possibly to turbulence, as they do in 2D shear flows.

For $C_R \in \{2, 3\}$, nonlinear 3D-3C unsteady simulations run for $Gr_{ac} > Gr_{ac}^c$ showed that the growth and saturation of the initial white-noise perturbation were followed by the development of a secondary instability. For $C_R = 2$ in particular, that secondary instability triggered a transition towards an unsteady state. This raises the question of whether this transition is the trace of an edge state, i.e., if there exists a threshold for the initial perturbation energy below which the saturated state transitions back to the initial steady axisymmetric state. Answering these questions would require further analysis.

Finally, the values of $Gr_{ac}^c$ we obtained correspond to relatively low forcing magnitudes. If we use the length, mass and time scales of the 2 MHz water experiments of Moudjed et al. (2014a) made with a 28.5 mm-in-diameter source, the instability onsets correspond to critical acoustic powers $P_{ac}^c$ ranging from 0.11 to 0.78 W. Such $P_{ac}^c$ values lie at the low end of the ranges commonly explored in water experiments. For instance, Moudjed et al. (2014a) did their measurements for $0.7 \ W \le P_{ac} \le 5.6 \ W$, Frenkel et al. (2001) for $2 \ W \le P_{ac} \le 10 \ W$, and Kamakura et al. (1996) for $P_{ac} = 0.68 \ W$. Besides shedding light on the rather unstable nature of acoustic streaming jets, our study also demonstrates the relevance of numerical simulations and LSA to capture the onset of instability, which is difficult to catch in water experiments. This work is thus an important milestone towards the thorough understanding of the rich dynamics of streaming flows (Cambonie et al. 2017; Launay et al. 2019) and their transition to turbulence. Conversely, the strong variation of $Gr_{ac}^c$ with confinement offers a convenient way to manipulate the nature of the flow forced in specific configurations: where turbulence is desired for mixing or other purposes, immersing the acoustic transducer with a plate cutting the jet in an open region of the flow would favour the instabilities conducive to turbulence. If, by contrast, a well-controlled laminar flow is sought, the transducer could be housed at one end of a thin cylinder open towards the forcing region at the other end: such a configuration would suppress the shear instability at the edge of the jet whilst still injecting the required amount of momentum into the flow. The analysis carried out in this paper offers predictive tools that may be easily adapted to design such a device. The implementation of such a technique to control the flow state raises further fundamental questions: The transducer is usually switched on to its nominal power from 0, rather than gently ramped up, which may drive a complex transient. The nature of this transient and its stability is also likely to depend on the initial flow conditions, at the point the transducer is switched on. Specific stability analyses based on adjoint methods would be required to understand how. In the end, if the desired effect is to maintain the flow in either a laminar or turbulent state, a form of control may be required. The topology of the acoustic force may indeed be optimised to that end, also by means of adjoint methods (Luchini & Bottaro 2014), as recently proposed in mixed baroclinic flows relevant to the continuous casting of Aluminium (Kumar & Pothérat 2025).

**Supplementary data.** Supplementary movies are available at ?????

**Acknowledgements.** The support from the PMCS2I (Ecole Centrale de Lyon) and from the High Performance Computing Centre of the Faculty of Engineering, Environment and Computing (Coventry University) for the numerical calculations is gratefully acknowledged. The authors would like to thank Laurent Pouilloux (Ecole Centrale de Lyon) and Alex Pedcenko (Coventry University) for their availability and for providing help at any stage of the project. The authors also wish to thank Hugh M. Blackburn for the support on the use of the spectral-element code Semtex.
827 For the purpose of Open Access, a CC-BY public copyright licence has been applied by the authors to the
828 present document and will be applied to all subsequent versions up to the Author Accepted Manuscript
829 arising from this submission.

830 **Funding.** This work was carried out as part of the BRASSOA project supported by the Institut Carnot
831 Ingénierie@Lyon and by a Royal Society International Exchange grant (Ref. IES\R2\202212). AP
832 acknowledges support from EPSRC through grant No. EP/X010937/1.

833 **Declaration of interests.** The authors report no conflict of interest.

834 **Author ORCIDs.** B. Vincent, https://orcid.org/0000-0002-9937-7494; A. Kumar, https://orcid.org/0000-0002-6026-5727; D. Henry, https://orcid.org/0000-0002-7231-7918; S. Miralles, https://orcid.org/0000-0002-4701-8609; V. Botton, https://orcid.org/0000-0002-0180-3089; A. Pothérat, https://orcid.org/0000-0001-8691-5241

838 **Appendix A. Effect of under-resolving the near-acoustic field on the numerical**
839 **precision of the acoustic force**

840 The numerical meshes used in this study are too coarse to precisely resolve the near-acoustic
841 field in the region $x \le \frac{1.22}{4S}$, as illustrated in figure 3. To determine the spatial resolution of $\tilde{I}$
842 for simulating streaming jets, we start with an important theoretical result originally derived
843 by Lighthill (1978) as the simplest reference for comparison. Lighthill derived an analytical
844 expression for the axial momentum flux $M_{\text{th}}$ of a steady axisymmetric streaming jet flowing
845 in a semi-infinite $x \ge 0$ domain, without any pressure forces. Using the transducer diameter
846 $D_s$ as length scale, $\rho D_s^3$ as mass scale ($\rho$ is the fluid density) and $D_s^2/\nu$ as time scale ($\nu$ is
847 the kinematic viscosity), Lighthill's result reads:

$$
848 M_{\text{th}} = \int_{\theta=0}^{2\pi} \int_{r=0}^{+\infty} u_x^2 r \, dr \, d\theta = M \left(1 - e^{-2 \frac{N}{L} x}\right), \quad (\text{A 1})
$$

849 where

$$
850 \qquad M = \frac{\pi L G r_{ac}}{32N} \tag{A 2}
$$

851 is the maximum axial momentum flux, achieved in the limit $x \to +\infty$. To determine whether
852 undersampling the very near acoustic field was acceptable, we computed the undersampled
853 momentum flux $M_{\text{num}}$ at each $x$ and for each case:

$$
854 \qquad M_{\text{num}} = \int_{\theta=0}^{2\pi} \int_{r=0}^{D/2} u_x^2 r \, dr \, d\theta, \qquad (\text{A 3})
$$

855 and compared its value to Lighthill's theoretical prediction (A 1). The profiles of $M_{\text{num}}/M$
856 are shown in figure 23 for all $(N, C_R)$ and for $1000 \le Gr_{ac} \le 15500$. The $M_{\text{th}}$ profile is
857 also reported in black. For the least confined case (i.e., $(N, C_R) = (1, 6)$, plotted with stars),
858 there is an excellent agreement between $M_{\text{num}}$ and $M_{\text{th}}$: the two profiles differ by only 2.2 %
859 on $x \in [0, 0.75L]$. This excellent agreement means that the injected momentum is correct,
860 despite the near acoustic field being undersampled. The larger discrepancy observed for
861 $x \in [0.75L, L]$ (blue-shaded area in figure 23) is caused by the recirculation imposed by the
862 downstream wall; such wall is absent in the free-jet setting for which Lighthill's theoretical
863 prediction (A 1) is valid.

864 However, further confining the flow increases the discrepancy between $M_{\text{num}}$ and $M_{\text{th}}$ (see
865 inset of figure 23). For $(N, C_R) = (0.25, 6)$ and $Gr_{ac} = 4750$, that discrepancy is 3.5 % on
866 $x \in [0, 75L]$. In the most confined setting $(N, C_R) = (0.25, 1)$ with $Gr_{ac} = 6400$, $M_{\text{num}}$
867 and $M_{\text{th}}$ differ by 53.0 % on $x \in [0, 0.75L]$. The increase of that discrepancy with the flow
868 confinement has several causes:
<figure><img src="image_31.png" /></figure>

Figure 23: Profiles of the numerical axial momentum flux $M_{\text{num}}$ for all cases and different $Gr_{ac}$. The profiles are normalised by the theoretical maximum axial momentum flux $M$ in an infinitely-long cavity (A 2)). The profile of the theoretical flux $M_{\text{th}}$, as originally derived by Lighthill (J. Sound Vib., 1987), is shown in black. The blue-shaded regions correspond to the last 25 % of the cavity length for the $N = 0.25$ (left region) and $N = 1$ (right region). These regions are affected by the recirculation imposed by the downstream wall. Inset: detailed view of the region $x \in [0, 50]$.

(i) Laterally confining the flow (i.e., reducing $C_R$) increases the adverse pressure gradient, as shown in figure 6 of our manuscript. For the most confined cases, the magnitude of the pressure forces becomes comparable to the magnitude of the acoustic force. The pressure forces can no longer be ignored, thus invalidating Lighthill's model (A 1).

(ii) The adverse pressure gradient is responsible for driving the return flow. As we reduce $C_R$, the return flow is accelerated and increases shear on the lateral wall. These shear forces are not accounted for in Lighthill's model (A 1), which assumes the fluid to be at rest as $r \to +\infty$.

(iii) All cavity walls are sound-absorbing walls, meaning that the acoustic field is truncated by the walls. Reducing $C_R$ increases the loss of acoustic power caused by truncation of the sound field. As a consequence, the RHS of Lighthill's model overestimates the net acoustic force in a confined setup.

For the most confined setups, the theoretical model of Lighthill (A 1) can thus no longer be used to determine whether we inject the right amount of momentum within the fluid.

To overcome the limitations of Lighthill's model, we shall now account for the pressure, acoustic and viscous forces in the momentum balance. The control volume is a cylinder of length $x$, diameter $D$ (i.e., the same diameter as the cavity) and whose base lies at $x = 0$. Owing to the no-slip boundary condition on the walls, the momentum balance along the
<figure><img src="image_32.png" /><figcaption>Figure 24: Profiles of $M_{\text{RHS}}$, defined as the right-hand side of the momentum balance (A 4). All cases are represented for several values of $Gr_{ac}$. The black curve represents the exact expression of $M_{\text{RHS}}$ in a semi-infinite fluid domain, whilst the points represented by filled symbols are obtained from the simulation results. The blue-shaded areas represent the last 25 % of the cavity length for the $N = 0.25$ (left area) and $N = 1$ (right area) cases. Inset: detailed view at $x \in [0, 50]$.</figcaption></figure>

axial direction reads:

$$ \underbrace{M \left(1 - e^{-2Nx/L}\right)}_{\text{Net acoustic force in a semi-infinite domain}} \approx \underbrace{\int_{0}^{2\pi} \int_{0}^{D/2} u_{x}^{2}(x, r) r \, dr \, d\theta}_{\text{Inertia}} + \underbrace{\int_{0}^{2\pi} \int_{0}^{D/2} [p(x, r) - p(0, r)] r \, dr \, d\theta}_{\text{Pressure forces}} \\ - \underbrace{\int_{0}^{x} \int_{0}^{2\pi} \frac{D}{2} \frac{\partial u_{x}}{\partial r}(x', \frac{D}{2}) \, dx' \, d\theta}_{\text{Viscous forces on the lateral wall}} \qquad (\text{A 4}) \\ - \underbrace{\int_{0}^{2\pi} \int_{0}^{D/2} 2 \left[ \frac{\partial u_{x}}{\partial x}(x, r) - \frac{\partial u_{x}}{\partial x}(0, r) \right] r \, dr \, d\theta}_{\text{Normal viscous forces}} $$

where we replaced the volume integral of the acoustic force by its exact expression (A 1) for a continuous force field in a semi-infinite domain. Denoting the right-hand side of the momentum balance (A 4) by $M_{\text{RHS}}$, we check whether the momentum balance (A 4) is satisfied. If so, this confirms that the acoustic field is sufficiently well discretised to capture the correct injection of momentum into the fluid volume.

Figure 24 compares the profiles of $M_{\text{RHS}}$ to the approximate expression $M(1 - e^{-2Nx/L})$ of the integrated force term. All data points collapse on the theoretical net acoustic force (i.e., the LHS of equation (A 4)) curve shown in black. The maximum discrepancy between
898 the numerics and the theory falls down to 3.9 % for the ($N, C_R$) = (0.25, 1) case with $Gr_{ac}$ = 15500. Note that this is a rather pessimistic estimate, as the analytical expression of the net acoustic force does not account for the truncation of the acoustic field by the lateral walls. Therefore, this excellent agreement confirms that the solution with undersampled near-field captures the correct momentum injected into the fluid.

903 To conclude, we used momentum balances to verify that the spatial resolution of the forcing field is sufficient to capture the correct momentum source due to the acoustic force within the fluid domain. The forcing thus does not need to be spatially-averaged in the near acoustic field.

REFERENCES

907 ABSAR, S., PASUMARTHI, P. & CHOI, H. 2017 Numerical and experimental studies about the effect of acoustic streaming on ultrasonic processing of metal matrix nanocomposites (MMNCs). *J. Manuf. Process.* **28**, 515–522.

908 ALGHANE, M., FU, Y. Q., CHEN, B. X., LI, Y., DESMULLIEZ, M. P. Y. & WALTON, A. J. 2012 Frequency effect on streaming phenomenon induced by Rayleigh surface acoustic wave in microdroplets. *J. Appl. Phys.* **112**.

910 BARKLEY, D., BLACKBURN, H. M. & SHERWIN, S. J. 2008 Direct optimal growth analysis for timesteppers. *Intl J. Numer. Meth. Fluids* **57**, 1435–1458.

915 BEN HADID, H., DRIDI, W., BOTTON, V., MOUDJED, B. & HENRY, D. 2012 Instabilities in the Rayleigh-Bénard-Eckart problem. *Phys. Rev. E* **86**, 016312.

917 BLACKBURN, H.M., LEE, D., ALBRECHT, T. & SINGH, J. 2019 Semtex: a spectral element-Fourier solver for the incompressible Navier-Stokes equations in cylindrical or cartesian coordinates. *Comput. Phys. Commun.* **245**, 106804.

920 BLACKBURN, H.M. & SHERWIN, S.J. 2004 Formulation of a Galerkin spectral element-Fourier method for three-dimensional incompressible flows in cylindrical geometries. *J. Comput. Phys.* **197**, 759–778.

922 BLACKSTOCK, D.T. 2000 *Fundamentals of Physical Acoustics*. John Wiley & Sons, Inc.

923 CAMBONIE, T., MOUDJED, B., BOTTON, V., HENRY, D. & BEN HADID, H. 2017 From flying wheel to square flow: Dynamics of a flow driven by acoustic forcing. *Phys. Rev. Fluids* **2**, 123901.

925 CAMOBRECO, C.J., POTHÉRAT, A. & SHEARD, G.J 2020 Subcritical route to turbulence via the Orr mechanism in a quasi-two-dimensional boundary layer. *Phys. Rev. Fluids* **5** (11), 113902.

927 CAMOBRECO, C.J., POTHÉRAT, A. & SHEARD, G.J 2021 Transition to turbulence in quasi-two-dimensional MHD flow driven by lateral walls. *Phys. Rev. Fluids* **6** (1), 013901.

929 CAMOBRECO, C. J., POTHÉRAT, A. & SHEARD, G. J. 2023 Subcritical transition to turbulence in quasi-two-dimensional shear flows. *J. Fluid Mech.* **963**, R3.

931 CARMO, B. S., SHERWIN, S. J., BEARMAN, P. W. & WILLDEN, R. H. J. 2008 Wake transition in the flow around two circular cylinders in staggered arrangements. *J. Fluid Mech.* **597**, 1–29.

933 CHARRIER-MOJTABI, M-C., FONTAINE, A. & MOJTABI, A. 2012 Influence of acoustic streaming on thermo-diffusion in a binary mixture under microgravity. *Intl. J. Heat Mass Transfer* **55**, 5992–5999.

935 CHARRIER-MOJTABI, M-C., JACOB, X., DOCHY, T. & MOJTABI, A. 2019 Species separation of a binary mixture under acoustic streaming. *Eur. Phys. J. E* **42**, 64.

937 CHONG, M. S., PERRY, A. E. & CANTWELL, B. J. 1990 A general classification of three-dimensional flow fields. *Phys. Fluids A* **2**, 765–777.

939 CLEVE, S., GUÉDRA, M., MAUGER, C., INSERRA, C. & BLANC-BENON, P. 2019 Microstreaming induced by acoustically trapped, non-spherically oscillating microbubbles. *J. Fluid Mech.* **875**, 597–621.

941 CLUNE, T. & KNOBLOCH, E. 1993 Pattern selection in rotating convection with experimental boundary conditions. *Phys. Rev. E* **47**, 2536–2550.

943 COULOUVRAT, FRANÇOIS 1992 On the equations of nonlinear acoustics. *Journal Acoustique* **5**, 321–359.

944 DÉLERY, J. M. 2001 Robert Legendre and Henri Werlé: Toward the elucidation of three-dimensional separation. *Ann. Rev. Fluid Mech.* **33**, 129–154.

946 DENTRY, M. B., YEO, L. Y. & FRIEND, J. R. 2014 Frequency effects on the scale and behavior of acoustic streaming. *Phys. Rev. E* **89**, 013203.

948 DOINIKOV, A. A., REGNAULT, G., MAUGER, C., BLANC-BENON, P. & INSERRA, C. 2022 Acoustic
949 microstreaming produced by two interacting gas bubbles undergoing axisymmetric shape oscillations.

950 *J. Fluid Mech.* **931**, A19.

951 DOUSSET, V. & POTHÉRAT, A. 2010 Formation mechanism of hairpin vortices in the wake of a truncated
cylinder in a duct. *J. Fluid. Mech.* **653**, 519–536.

952
953 DOUSSET, V. & POTHÉRAT, A. 2012 Characterisation of the wake of a truncated cylinder in a duct under a
strong axial magnetic field. *J. Fluid Mech.* **691**, 341–367.

955 DRIDI, W., HENRY, D. & BEN HADID, H. 2008a Influence of acoustic streaming on the stability of a laterally
heated three-dimensional cavity. *Phys. Rev. E* **77**, 046311.

957 DRIDI, W., HENRY, D. & BEN HADID, H. 2008b Influence of acoustic streaming on the stability of melt flows
in horizontal Bridgman configurations. *J. Cryst. Growth* **310**, 1546–1551.

959 DRIDI, W., HENRY, D. & HADID, H. BEN 2010 Stability of buoyant convection in a layer submitted to acoustic
streaming. *Phys. Rev. E* **81**, 056309.

961 ECKART, C. 1948 Vortices and streams caused by sound waves. *Phys. Rev.* **73**, 68–76.

962 EL GHANI, N., MIRALLES, S., BOTTON, V., HENRY, D., BEN HADID, H., TER-OVANESSIAN, B. & MARCELIN,
S. 2021 Acoustic streaming enhanced mass transfer at a wall. *Intl. J. Heat Mass Transfer* **172**.

964 FAUCONNIER, M., MAUGER, C., BÉRA, J-C. & INSERRA, C. 2022 Nonspherical dynamics and microstreaming
of a wall-attached microbubble. *J. Fluid Mech.* **935**, A22.

966 FOSS, J. F. 2004 Surface selections and topological constraint evaluations for flow field analyses. *Exp. Fluids* **37**, 883–898.

968 FRENKEL, V., GURKA, R., LIBERZON, A., SHAVIT, U. & KIMMEL, E. 2001 Preliminary investigations of
ultrasound induced acoustic streaming using particle image velocimetry. *Ultrasonics* **39**, 153–156.

970 FRIEND, J. & YEO, L. Y. 2011 Microscale acoustofluidics: microfluidics driven via acoustics and ultrasonics.
*Rev. Mod. Phys.* **83**, 647–704.

972 FROMMELT, T., KOSTUR, M., WENZEL-SCHÄFER, M., TALKNER, P., HÄNGGI, P. & WIXFORTH, A. 2008
Microfluidic mixing via acoustically driven chaotic advection. *Phys. Rev. Lett.* **100**, 034502.

974 GREEN, A., MARSHALL, J. S., MA, D. & WU, J. 2016 Acoustic streaming and thermal instability of flow
generated by ultrasound in a cylindrical container. *Phys. Fluids* **28**, 104105.

976 HAGSÄTER, S. M., JENSEN, T. G., BRUUS, H. & KUTTER, J. P. 2007 Acoustic resonances in microfluidic chips:
full-image micro-PIV experiments and numerical simulations. *Lab Chip* **7**.

978 HAMILTON, M. F. & BLACKSTOCK, D. T. 2024 *Nonlinear Acoustics*. Springer Nature Switzerland.

979 HENDERSON, R. D. 1997 Nonlinear dynamics and pattern formation in turbulent wake transition. *J. Fluid Mech.* **352**, 65–112.

981 HENDERSON, R. D. & BARKLEY, D. 1996 Secondary instability in the wake of a circular cylinder. *Phys. Fluids* **8**, 1683–1685.

983 HENRY, D., VINCENT, B., MIRALLES, S., BOTTON, V. & HADID, H. BEN 2022 Steady and oscillatory flows
generated by Eckart streaming in the two-dimensional Rayleigh–Bénard configuration. *J. Fluid Mech.* **952**, A28.

986 HUANG, Z., GAO, R., GAO, YY & XI, G 2024 Subcritical transitional flow in two-dimensional plane poiseuille
flow. *J. Fluid Mech.* **994**, A6.

988 HUNT, J. C. R., ABELL, C. J., PETERKA, J. A. & WOO, H. 1978 Kinematical studies of the flows around free
or surface-mounted obstacles; applying topology to flow visualization. *J. Fluid Mech.* **86**, 179–200.

990 KAMAKURA, T., SUDO, T., MATSUDA, K. & KUMAMOTO, Y. 1996 Time evolution of acoustic streaming from
a planar ultrasound source. *J. Acoust. Soc. Am.* **100**, 132–138.

992 KARNIADAKIS, G. E., ISRAELI, M. & ORSZAG, S. A. 1991 High-order splitting methods for the incompressible
Navier-Stokes equations. *J. Comput. Phys.* **97**, 414–443.

994 KINSLER, L. E., FREY, A. R., COPPENS, A. B. & SANDERS, J. V. 2000 *Fundamentals of Acoustics*, 4th edn.
John Wiley & Sons, Inc.

996 KOZHEMYAKIN, G.N. 2003 Imaging of convection in a Czochralski crucible under ultrasound waves. *J. Cryst.*
*Growth* **257** (3), 237–244.

998 KUMAR, A. & POTHÉRAT, A. 2020 Mixed baroclinic convection in a cavity. *J. Fluid Mech.* **885**, A40.

999 KUMAR, A. & POTHÉRAT, A. 2025 Suppressing instabilities in mixed baroclinic flow using an actuation
based on receptivity. *J. Fluid Mech.* **1013**, A20.

1001 LANDAU, L. D. & LIFSCHITZ, E. M. 1987 *Fluid Mechanics*, 2nd edn. Pergamon Books Ltd.

1002 LAUNAY, G., CAMBONIE, T., HENRY, D., POTHÉRAT, A. & BOTTON, V. 2019 Transition to chaos in an
acoustically driven cavity flow. *Phys. Rev. Fluids* **4**, 044401.

1004 LEBON, G., SALLOUM-ABOU-JAOUDE, G., ESKIN, D., TZANAKIS, I., PERICLEOUS, K. & JARRY, P. 2019
1005 Numerical modelling of acoustic streaming during the ultrasonic melt treatment of direct-chill (DC) casting. *Ultrason. Sonochem.* **54**, 171–182.

1007 LIGHTHILL, J. 1978 Acoustic streaming. *J. Sound Vib.* **61**, 391–418.

1008 LUCHINI, P. & BOTTARO, A. 2014 Adjoint Equations in Stability Analysis. *Ann. Rev. Fluid Mech.* **46**, 493–517.

1009 LYUBIMOVA, T. P. & SKURIDIN, R. V. 2019 The effect of acoustic wave on the stability of stationary convective flow of binary fluid in a horizontal layer subjected to horizontal temperature and concentration gradients. *Intl. J. Heat Mass Transfer* **132**, 789–801.

1012 MITOME, H. 1998 The mechanism of generation of acoustic streaming. *Electron. Commun. Japan* **81**.

1013 MOFFATT, H. K. 1964 Viscous and resistive eddies near a sharp corner. *J. Fluid Mech.* **18**, 1–18.

1014 MONDAL, R. & MAHBUB ALAM, MD. 2023 Blockage effect on wakes of various bluff bodies: A review of confined flow. *Ocean Eng.* **286**, 115592.

1016 MOUDJED, B., BOTTON, V., HENRY, D., BEN HADID, H. & GARANDET, J.-P. 2014a Scaling and dimensional analysis of acoustic streaming jets. *Phys. Fluids* **26**, 093602.

1018 MOUDJED, B., BOTTON, V., HENRY, D., MILLET, S., GARANDET, J.-P. & BEN HADID, H. 2014b Oscillating acoustic streaming jet. *Appl. Phys. Lett.* **105**.

1020 MOUDJED, B., BOTTON, V., HENRY, D., MILLET, S., GARANDET, J.-P. & BEN HADID, H. 2015 Near-field acoustic streaming jet. *Phys. Rev. E* **91**, 033011.

1022 MULLER, P. B. & BRUUS, H. 2015 Theoretical study of time-dependent, ultrasound-induced acoustic streaming in microchannels. *Phys. Rev. E* **92**, 063018.

1024 ROSCO, J. & FRIEND, J. 2022 Modeling fast acoustic streaming: steady-state and transient flow solutions. *Phys. Rev. E* **106**, 045101.

1026 PIERCE, A. D. 1990 Wave equation for sound in fluids with unsteady inhomogeneous flow. *J. Acoust. Soc. Am.* **87**, 2292–2299.

1028 POTHÉRAT, A. AND ZHANG, L. 2018 Dean flow and vortex shedding in a three-dimensional 180 degrees sharp bend. *arXiv preprint arXiv:1807.10950*.

1030 QU, J., HENRY, D., MIRALLES, S., BOTTON, V. & RAYNAL, F. 2022 Chaotic mixing in an acoustically driven cavity flow. *Phys. Rev. Fluids* **7**, 064501.

1032 RAYLEIGH, L. 1884 On the circulation of air observed in Kundt's tubes, and on some allied acoustical problems. *Phil. Trans. R. Soc. London* **175**, 1–21.

1034 RIAUD, A., BAUDOIN, M., BOU MATAR, O., THOMAS, J-L. & BRUNET, P. 2017 On the influence of viscosity and caustics on acoustic streaming in sessile droplets: an experimental and a numerical study with a cost-effective method. *J. Fluid Mech.* **821**, 384–420.

1037 RUDENKO, O.V & SUKHORUKOV, A.A. 1998 Nonstationary Eckart Streaming and Pumping of Liquid in an Ultrasonic Field. *Acoust. Phys.* **44**, 565–570.

1039 SAPARDI, A. M., HUSSAM, W. K., POTHÉRAT, A. & SHEARD, G. J. 2017 Linear stability of confined flow around a 180-degree sharp bend. *J. Fluid Mech.* **822**, 813–847.

1041 SCHMID, P. J. 2007 Nonmodal stability theory. *Ann. Rev. Fluid Mech.* **39**, 129–162.

1042 SHEARD, G. J., THOMPSON, M. C. & HOURIGAN, K. 2004 From spheres to circular cylinders: non-axisymmetric transitions in the flow past rings. *J. Fluid Mech.* **506**, 45–78.

1044 SHERWIN, S. J. & BLACKBURN, H. M. 2005 Three-dimensional instabilities and transition of steady and pulsatile axisymmetric stenotic flows. *J. Fluid Mech.* **533**, 297–327.

1046 THOMPSON, M.C., LEWEKE, T. & PROVANSAL, M. 2001 Kinematics and dynamics of sphere wake transition. *J. Fluids Struct.* **15**, 575–585.

1048 TOUIHRI, R., BEN HADID, H. & HENRY, D. 1999 On the onset of convective instabilities in cylindrical cavities heated from below. II. Effect of a magnetic field. *Phys. Fluids* **11**, 2089–2100.

1050 VINCENT, B., HENRY, D., KUMAR, A., BOTTON, V., POTHÉRAT, A. & MIRALLES, S. 2024a Phenomenology of laminar acoustic streaming jets, arXiv: 2409.15083.

1052 VINCENT, B., MIRALLES, S., HENRY, D., BOTTON, V. & POTHÉRAT, A. 2024b Experimental study of a helical acoustic streaming flow. *Phys. Rev. Fluids* **9**, 024101.

1054 WESTERVELT, P. 1953 The theory of steady rotational flow generated by a sound field. *J. Acoust. Soc. Am.* **25**, 60–67.