1. Introduction
Trailing-edge noise is an overriding concern for design of quiet air vehicles. At low and moderate Reynolds numbers, tonal noise becomes an important component of the acoustic spectrum. Several studies on trailing-edge aeroacoustics were conducted starting in the 1970s to examine the tonal noise generation by airfoils (Paterson et al. Reference Paterson, Vogt, Fink and Munch1973; Tam Reference Tam1974; Fink Reference Fink1975; Longhouse Reference Longhouse1977; Arbey & Bataille Reference Arbey and Bataille1983). Noise measurements were performed by Paterson et al. (Reference Paterson, Vogt, Fink and Munch1973) for symmetric NACA airfoils at various angles of attack over a Reynolds number range between 
$10^{5}$ and 
$10^{6}$. Their results showed the existence of multiple tones in a ladder-like structure in terms of frequency and freestream velocity. These authors also found strong two-point correlation of surface pressure along the spanwise direction, which indicated that the flow phenomenon associated with tonal noise generation could be modelled as two-dimensional (2-D).
Tam (Reference Tam1974) suggested that the ladder-like structure is due to a self-excited feedback loop between disturbances in the boundary layer and the airfoil wake. Fink (Reference Fink1975) assumed that the discrete tonal frequencies are related to the laminar boundary layer on the pressure side. In order to elucidate aspects of airfoil noise, Arbey & Bataille (Reference Arbey and Bataille1983) performed experiments in an open wind tunnel for different NACA airfoils at 
$\alpha = 0^{\circ }$ for 
$1 \times 10^{5} \le Re \le 7 \times 10^{5}$. The aforementioned studies showed that the noise spectrum has a broadband component with a main tonal peak plus a set of equidistant secondary tones due to a feedback mechanism closing at the point of maximum flow velocity along the airfoil. The broadband component was assumed to appear due to scattering of Tollmien–Schlichting (TS) instabilities. For airfoils at incidence, Lowson, Fiddes & Nash (Reference Lowson, Fiddes and Nash1994) found that the presence of secondary tones was related to a separation bubble developed on the airfoil pressure side. In this case, TS instabilities developing along the laminar boundary layer would lead to acoustic scattering on the trailing edge and acoustic waves would then propagate upstream closing the feedback loop, with the separation bubble acting as an amplifier of acoustic disturbances.
Nash, Lowson & McAlpine (Reference Nash, Lowson and McAlpine1999) performed experimental studies of airfoil noise for a NACA0012 profile up to a Reynolds number of 
$1.45\times 10^{6}$ and several angles of attack. A closed-section wind tunnel, with and without acoustic-absorbing lining on its walls, was used in the experiments and results from the hard-wall tunnel revealed multiple frequency peaks. However, the authors argued that these tonal peaks were correlated to resonant frequencies of the wind tunnel. Thus, they carried out measurements with lined walls simulating anechoic conditions and a single dominant tone was observed instead of several peaks. Furthermore, no ladder-like structure of tonal frequency was observed, in disagreement with previous studies of Paterson et al. (Reference Paterson, Vogt, Fink and Munch1973), Fink (Reference Fink1975) and Arbey & Bataille (Reference Arbey and Bataille1983). It is important to mention that secondary tones were often observed in experiments conducted in open-jet facilities.
More recently, Plogmann, Herrig & Würz (Reference Plogmann, Herrig and Würz2013) also found multiple tones in their experiments (NACA0012, 
$3.1 \times 10^{5} \le Re \le 1.5 \times 10^{6}$, 
$0^{\circ } \le \alpha \le 9^{\circ }$) and demonstrated that tripping the pressure side boundary layer leads to turbulent flow, eliminating the separation bubble and the secondary tones. These authors emphasise that the feedback loop is extremely sensitive to small variations in the flow conditions that, in turn, lead to changes in the tonal components. This occurs particularly at higher Reynolds numbers, where the flow is more prone to transition. The dependency of angle of attack and Reynolds number on tonal noise emission is highlighted by Pröbsting, Scarano & Morris (Reference Pröbsting, Scarano and Morris2015). These authors performed experiments for a NACA0012 airfoil for 
$0.3\times 10^{5} \le Re \le 2.3\times 10^{5}$ and effective angles of attack 
$0^{\circ } \le \alpha \le 6.3^{\circ }$. Tripping devices were applied separately on each side of the airfoil to identify their respective role in the noise generation and it was found that suction-side (pressure-side) events dominate at lower (higher) Reynolds numbers. Moreover, it was observed that, at low angles of attack, interactions between the two sides of the airfoil become increasingly important.
As discussed previously, early experiments report the presence of a separation bubble on the airfoil pressure side. Most of these investigations were conducted at higher Reynolds numbers, where the flow was likely turbulent on the suction side. These observations are in agreement with the large eddy simulations (LESs) of a NACA0012 at 
$\alpha = 5^{\circ }$ and 
$Re = 4 \times 10^{5}$ performed by Wolf, Azevedo & Lele (Reference Wolf, Azevedo and Lele2012a) and Ricciardi, Ribeiro & Wolf (Reference Ricciardi, Ribeiro and Wolf2019). Flow visualisation and proper orthogonal decomposition were used in these references to identify coherent structures shed from the pressure side near the trailing edge. Such flow structures were found responsible for the intense tonal noise generation despite the fact that a turbulent boundary layer developed on the airfoil suction side. In this case, the pressure side boundary layer was laminar due to the favourable pressure gradient. In agreement with the experiments of Plogmann et al. (Reference Plogmann, Herrig and Würz2013), Wolf et al. (Reference Wolf, Azevedo and Lele2012a) showed that the tonal component vanishes when both boundary layers are tripped.
At lower Reynolds numbers, a laminar separation bubble (LSB) exists on the airfoil suction side and is responsible for the overall flow dynamics and noise generation. In the context of airfoil flows, this flow feature has been studied for different purposes. For instance, direct numerical simulations were performed by Jones, Sandberg & Sandham (Reference Jones, Sandberg and Sandham2008) for a NACA0012 airfoil at Reynolds number 
$Re = 5 \times 10^{4}$, 
$M_\infty = 0.4$ and 
$\alpha = 5^{\circ }$. It was shown that, despite being absolutely stable by means of linear stability analysis, turbulence is self-sustained even with the absence of forcing. Later, experimental investigations by Pröbsting & Yarusevych (Reference Pröbsting and Yarusevych2015) with a NACA0012 airfoil at moderate Reynolds numbers 
$0.65 \times 10^{5} \le Re \le 4.5 \times 10^{5}$ and 
$\alpha = 2^{\circ }$ present intermittent laminar–turbulent transition that affects the advection of coherent structures from the LSB towards the trailing edge. Further analyses were also presented by Kurelek, Lambert & Yarusevych (Reference Kurelek, Lambert and Yarusevych2016), Kurelek, Kotsonis & Yarusevych (Reference Kurelek, Kotsonis and Yarusevych2018), Kurelek, Yarusevych & Kotsonis (Reference Kurelek, Yarusevych and Kotsonis2019), Michelis, Yarusevych & Kotsonis (Reference Michelis, Yarusevych and Kotsonis2018) and Pröbsting & Yarusevych (Reference Pröbsting and Yarusevych2021) to study the impact of acoustic excitation, three-dimensional (3-D) effects and shedding/merging of vortices from the suction side LSB.
The first numerical simulations on airfoil secondary tones were conducted by Desquesnes, Terracol & Sagaut (Reference Desquesnes, Terracol and Sagaut2007) considering 2-D flows over a NACA0012 airfoil for a Reynolds number 
$Re = 1\times 10^{5}$ at an angle of attack 
$\alpha = 5^{\circ }$, and for 
$Re = 2\times 10^{5}$ at 
$\alpha = 2^{\circ }$. In agreement with most experimental observations, multiple tonal peaks were observed. The previous authors also performed a local stability analysis assuming parallel flow and showed that the main tone frequency radiated to the far-field was close to that most amplified along the pressure side boundary layer. An assessment of the linear dynamics of wavepackets driving the feedback loop mechanism and the flow receptivity to acoustic forcing was investigated by time-marching the linearised Navier–Stokes (LNS) equations by Jones, Sandberg & Sandham (Reference Jones, Sandberg and Sandham2010) for a NACA0012 airfoil at 
$Re = 5\times 10^{4}$, 
$M_\infty = 0.4$ and 
$\alpha = 5^{\circ }$. Fosas de Pando, Schmid & Sipp (Reference Fosas de Pando, Schmid and Sipp2014b) applied bi-global stability analysis to investigate the dynamics coupling the boundary layers on both airfoil sides and the wake. The authors studied the 2-D flow over a NACA0012 airfoil for 
$Re = 2\times 10^{5}$, 
$M_\infty = 0.4$ and 
$\alpha = 2^{\circ }$ and found multiple frequencies in the eigenspectrum related not only to the main tonal peak but also to the secondary tones. Prior to this work, the stability analyses were limited a parallel flow assumption, and only a discussion of the dominant tonal frequency was presented. By means of adjoint and resolvent analyses, Fosas de Pando, Schmid & Lele (Reference Fosas de Pando, Schmid and Lele2014a) and Fosas de Pando, Schmid & Sipp (Reference Fosas de Pando, Schmid and Sipp2017) also identified sensitive regions of the flow which are prone to close the feedback loop mechanism. Hence, linear stability theory has proven to be an important methodology to investigate the generation of airfoil secondary tones and the feedback loop mechanism.
Another important observation from Desquesnes et al. (Reference Desquesnes, Terracol and Sagaut2007) regards the amplitude modulation of velocity fluctuations computed near the trailing edge. These authors discussed that such modulation is caused by interference of vortical structures from both sides of the airfoil which, combined with the feedback loop mechanism, would lead to the presence of multiple tones. Following the discussion on modulation of flow structures, Pröbsting, Serpieri & Scarano (Reference Pröbsting, Serpieri and Scarano2014) employed particle image velocimetry to study airfoils at Reynolds numbers 
$1\times 10^{5} \le Re \le 2.7\times 10^{5}$ and 
$2^{\circ } \le \alpha \le 4^{\circ }$. They investigated the mechanisms associated with tonal noise generation and the interference effects between suction and pressure sides of the airfoil. For 
$Re \approx 1.5\times 10^{5}$, at the lowest angle of attack, the authors showed that the amplitude modulation discussed by Desquesnes et al. (Reference Desquesnes, Terracol and Sagaut2007) was related to destructive interference of boundary layer instabilities occurring on both sides of the trailing edge. More recently, Ricciardi, Arias-Ramirez & Wolf (Reference Ricciardi, Arias-Ramirez and Wolf2020) performed 2-D simulations of a NACA0012 at 
$Re = 1.0\times 10^{5}$ at 
$\alpha = 3^{\circ }$ and showed that the multiple tones are related to modulation of the vortical structures developing on the suction side. In this case, the instantaneous main frequency alternates in time due to phase modulation of the flow structures shed by the LSB. As a consequence, multiple equidistant frequencies must appear in the Fourier transform to reconstruct the modulated signals. However, as pointed by Pröbsting et al. (Reference Pröbsting, Serpieri and Scarano2014), flow transition on the suction side occurs in experiments and their conclusions may be different from those observed in 2-D numerical simulations. In this regard, the simulation of Sanjose et al. (Reference Sanjose, Towne, Jaiswal, Moreau, Lele and Mann2019) for a thin-cambered airfoil at 
$Re = 1.5 \times 10^{5}$ exhibits flow transition where intermittency played a key role in the flow dynamics and noise emission for the configuration investigated. Thus, experimental (Pröbsting et al. Reference Pröbsting, Serpieri and Scarano2014) and numerical (Sanjose et al. Reference Sanjose, Towne, Jaiswal, Moreau, Lele and Mann2019; Nguyen et al. Reference Nguyen, Golubev, Mankbadi, Yakhina and Roger2021) findings contradict previous assumptions that 2-D simulations are sufficient to explain all the mechanisms of secondary tones.
In this study, a LES is performed to investigate a NACA0012 airfoil at 
$\alpha = 3^{\circ }$ immersed in a freestream flow with Mach number 
$M_{\infty } = 0.3$ and Reynolds number 
$Re = 5 \times 10^{4}$. This condition is chosen based on experimental results from Pröbsting et al. (Reference Pröbsting, Scarano and Morris2015) and Pröbsting & Yarusevych (Reference Pröbsting and Yarusevych2015). These authors showed that, for this flow configuration, suction-side events are responsible for the noise generation without a contribution from the pressure side. Moreover, the previous references show that the flow behaviour at this particular angle of attack is somehow independent of the Reynolds number up to a point where the relevant vortex dynamics switches to the pressure side. For the current flow setup, some key aspects are investigated including the dynamics of the suction side separation bubble and its vortex shedding, in addition to flow intermittency and its effect on amplitude modulation. Linear stability theory is applied to investigate the presence of sensitive regions and amplification mechanisms within the flow.
The present work is organised as follows: § 3.1 presents the LSBs which are observed over the suction side and near the trailing edge, on the pressure side. Then, flow visualisation is employed in § 3.2 to investigate vortex shedding aft of the suction side separation bubble. It is shown that vortex interaction results in merging or bursting such that coherent structures or turbulent packets are advected towards the trailing edge, leading to different levels of noise emission. Despite the intermittent occurrence of laminar–turbulent transition, the noise spectrum presents multiple equidistant secondary tones as discussed in § 3.3. To understand the role of flow instabilities on tonal noise generation, in § 3.4, the LNS equations are examined in its modal form using bi-global stability and resolvent analyses, and by the time evolution of disturbances which results in periodic wavepackets. Aspects of the acoustic feedback loop mechanism in both the linear and nonlinear solutions are examined in § 3.5 and we highlight the role of leading-edge disturbances from secondary diffraction in the closure of the acoustic feedback loop mechanism. The intermittent transition which leads to amplitude modulation of noise generation is analysed in § 3.6 in terms of phase interference of the dominant frequencies from hydrodynamic instabilities. Then, the main findings and conclusions are presented in § 4.
2. Theoretical and numerical approaches
2.1. Large eddy simulation
Large eddy simulation is performed to solve the compressible Navier–Stokes equations in general curvilinear coordinates. The spatial discretisation of the governing equations employs a sixth-order accurate compact scheme for derivatives and interpolations on a staggered grid (Nagarajan, Lele & Ferziger Reference Nagarajan, Lele and Ferziger2003). The time integration is carried out by a hybrid implicit–explicit method. The implicit second-order scheme of Beam & Warming (Reference Beam and Warming1978) is applied in the near-wall region to overcome the stiffness problem due to a fine boundary layer grid, whereas a third-order Runge–Kutta scheme is used for time advancement of the equations in flow regions away from the solid boundary. For the communication across the different methods, information is exchanged at overlapping points. An explicit sub-grid-scale model is not applied. However, outside the boundary layer, a sixth-order compact filter (Lele Reference Lele1992) is applied to control high-wavenumber numerical instabilities arising from grid stretching and interpolation between staggered grids. The transfer function associated with such filters has been shown to provide an approximation to sub-grid-scale models (Mathew et al. Reference Mathew, Lechner, Foysi, Sesterhenn and Friedrich2003).
No-slip adiabatic wall boundary conditions are enforced along the airfoil surface and characteristic plus sponge boundary conditions are applied in the far-field locations to minimise wave reflections (Wolf Reference Wolf2011). Periodic boundary conditions are used in the spanwise direction. Length scales, velocity components, density, pressure and temperature are non-dimensionalised as 
$\boldsymbol {x} = \boldsymbol {x}^{*}/L^{*}$, 
$\boldsymbol {u}=\boldsymbol {u}^{*}/a_{\infty }^{*}$, 
$\rho =\rho ^{*}/{\rho }_{\infty }^{*}$, 
$p = p^{*}/{\rho }_{\infty }^{*}{a_{\infty }^{*}}^{2}$ and 
$T=T^{*}/[(\gamma -1)T_{\infty }^{*}]$, respectively. Here, 
$L^{*}$ is the airfoil chord, 
$a_{\infty }^{*}$ is the freestream speed of sound, 
$\rho _{\infty }^{*}$ is the freestream density, 
$T_{\infty }^{*}$ is the freestream temperature and 
$\gamma$ is the ratio of specific heats. Variables with superscript 
$^{*}$ are given in dimensional units. Herein, time and frequency (Strouhal number) are presented non-dimensionalised by freestream velocity as 
$t=t^{*} U_{\infty }^{*} / L^{*}$ and 
$St=f^{*} L^{*} / U_{\infty }^{*}$, respectively. The present numerical methodology has been extensively validated for various 2-D and 3-D simulations of compressible airfoil flows at different configurations (Wolf et al. Reference Wolf, Azevedo and Lele2012a,Reference Wolf, Lele, Jothiprasard and Cheungb; Wolf, Azevedo & Lele Reference Wolf, Azevedo and Lele2013; Ramos et al. Reference Ramos, Wolf, Yeh and Taira2019).
The O-grid employed for the current LES is shown in grey lines for every three points in figure 1(a). For a smooth O-grid generation, the original airfoil trailing edge is truncated at 98 % of the chord and it has a curvature radius 
$r/L^{*} = 0.4\,\%$. The reference length scale 
$L^{*}$ is the unit chord from the original NACA0012 airfoil. The leading edge is placed at 
$(x, y) = (0, 0)$ and the airfoil is pivoted about this point. The total wingspan has 
$0.4 L^{*}$, being larger compared with previously reported studies of airfoil flows for this range of Reynolds number (Jones et al. Reference Jones, Sandberg and Sandham2008; Ducoin, Loiseau & Robinet Reference Ducoin, Loiseau and Robinet2016). The domain extends 37 chord-lengths outwards and a circular sponge, given by 
$A[(r - r_0)/\Delta r]^{n}$, occupies the last 
$\Delta r = 10$ chords. In this function, 
$r$ is the radial distance, 
$r_0$ is the sponge starting position and its coefficients are 
$A = 20$ and 
$n=4$.
Computational grids near the airfoil for (a) the O-mesh employed in the LES and (b) the H-mesh used in the linear stability analysis. The grids are shown for every three points in both directions.
The mesh has a distribution of points in the streamwise, wall-normal and spanwise directions given by 
$n_x=660$, 
$n_y=600$ and 
$n_z=192$, respectively, which results in approximately 
$76 \times 10^{6}$ grid points. The ratio of grid points along the suction side relative to the pressure side is approximately 5:3. The wall-normal distance of the first grid point is 
$\Delta n = 0.0001$ and the stretching ratio is 1.5 %. As discussed by Desquesnes et al. (Reference Desquesnes, Terracol and Sagaut2007), Jones & Sandberg (Reference Jones and Sandberg2011) and Fosas de Pando et al. (Reference Fosas de Pando, Schmid and Sipp2014b), the important mechanisms for tonal noise generation arise from 2-D flow instabilities in the laminar region of the flow. Motivated by this observation, a grid refinement study in terms of mean and fluctuation properties was conducted for 2-D simulations and shown by the present authors in Ricciardi, Wolf & Taira (Reference Ricciardi, Wolf and Taira2021). However, in the present study, the flow is 3-D and transitions to turbulence on the airfoil suction side, near the trailing edge. Considering only the turbulent flow region, the estimated mesh resolution in terms of wall units is given by 
$\Delta x^{+} < 10$, 
$\Delta y^{+} \approx 0.3$ and 
$\Delta z^{+} < 5$, which follow the best practices for wall-resolved LES. The time step is 
$\Delta t = 1.5 \times 10^{-4}$ and 75 convective time units are employed for post-processing and analysis of results. The 3-D simulation starts from a 2-D flow superposed with random noise and more than 30 convective time units are discarded before collecting flow statistics.
2.2. Linear stability and resolvent analyses
The mechanisms of generation and amplification of flow instabilities can be explained by linear stability theory. Hence, bi-global linear stability and resolvent analyses are employed in this work to investigate the most unstable frequencies and their sensitivity to disturbances, in addition to their role in the tonal noise generation. Through a Reynolds decomposition, it is possible to split the unsteady flow 
$\boldsymbol {q}(\boldsymbol {x},t)$ in a mean base flow 
$\bar {\boldsymbol {q}}(\boldsymbol {x})$ plus a time-dependent fluctuation component 
$\boldsymbol {q}'(\boldsymbol {x},t)$. If the fluctuations are sufficiently small, the Navier–Stokes equations can be linearised about the (mean) base flow. Although it is possible to consider the turbulent mean flow 
$\bar {\boldsymbol {q}}$ as the base flow, the linear stability analysis would not hold because such state is not an equilibrium point of the Navier–Stokes equations. Nonetheless, the use of a time-averaged base flow may provide some insights as a model (Taira et al. Reference Taira, Brunton, Dawson, Rowley, Colonius, McKeon, Schmidt, Gordeyev, Theofilis and Ukeiley2017, Reference Taira, Hemati, Brunton, Sun, Duraisamy, Bagheri, Dawson and Yeh2020). With addition of a time-dependent external forcing 
$\boldsymbol {f}'(\boldsymbol {x},t)$, the LNS equations can be written as
\begin{equation} \frac{\partial \boldsymbol{q}'}{\partial t} = \boldsymbol{\mathcal{L}} \boldsymbol{q}' + \boldsymbol{\mathcal{B}} \boldsymbol{f}' ,\end{equation}
where 
$\boldsymbol {\mathcal {B}} = \boldsymbol {\mathcal {B}}(\boldsymbol {x})$ is an operator which serves as a spatial window where forcing is applied. The matrix 
$\boldsymbol {\mathcal {L}} = \boldsymbol {\mathcal {L}}(\bar {\boldsymbol {q}})$ is the bi-global LNS operator with the spanwise time-averaged flow properties 
$\bar {\boldsymbol {q}}(\boldsymbol {x})$ as base flow.
The evolution of linear disturbances by (2.1) can be performed by a direct analysis of the operator 
$\boldsymbol {\mathcal {L}}$ or time-integrating the disturbances in a linearised version of the computational fluid dynamics (CFD) code. In the former case, with a transformation
\begin{equation} \bullet'(x,y,z,t) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \hat{\bullet}(x,y,\beta,\omega) \, {\rm e}^{\mathrm{i} (\beta z-\omega t)} \, {\rm d} \beta \, {\rm d} \omega , \end{equation}
where the wavenumber 
$\beta \in \mathbb {R}$ and the frequency 
$\omega \in \mathbb {C}$, it is possible to write the forced LNS equations in discrete form as
\begin{equation} {-}\mathrm{i} \omega \hat{\boldsymbol{q}} = \boldsymbol{\mathsf{L}} \hat{\boldsymbol{q}} + \boldsymbol{\mathsf{B}} \hat{\boldsymbol{f}} . \end{equation}
In this case, the linear operator becomes a function of the spanwise wavenumber 
$\boldsymbol{\mathsf{L}} = \boldsymbol{\mathsf{L}}(\bar {\boldsymbol {q}},\beta )$ and 
$\boldsymbol{\mathsf{B}}$ is the discrete version of 
$\boldsymbol {\mathcal {B}}$. If forcing 
$\hat {\boldsymbol {f}}$ is absent, the LNS equations can be analysed separately for each wavenumber 
$\beta$ as an eigenvalue problem (linear stability analysis) as
\begin{equation} \boldsymbol{\mathsf{Q}} \boldsymbol{\varLambda} = \boldsymbol{\mathsf{L}} \boldsymbol{\mathsf{Q}} .\end{equation}
Here, 
$\boldsymbol{\mathsf{Q}}$ holds the eigenvectors of 
$\boldsymbol{\mathsf{L}}$ and the eigenvalues appear in the diagonal matrix 
$\boldsymbol {\varLambda } = -\mathrm {i} \omega \boldsymbol{\mathsf{I}}$, where the frequency and growth rate are the real and imaginary parts of 
$\omega$, respectively. Due to the velocity gradients appearing off-diagonal, the linear operator becomes a non-normal matrix such that 
$\boldsymbol{\mathsf{L}}^{\mathcal {H}} \boldsymbol{\mathsf{L}} \ne \boldsymbol{\mathsf{L}} \boldsymbol{\mathsf{L}}^{\mathcal {H}}$. As a consequence, the eigenvectors are also non-normal. In this equation, the superscript 
$\mathcal {H}$ denotes the complex conjugate transpose (Hermitian). The matrix 
$\boldsymbol{\mathsf{L}}^{\mathcal {H}}$ is the adjoint operator, being related to regions of sensitivity within the flow. Its eigenvalue decomposition yields
\begin{equation} \boldsymbol{\mathsf{Q}}^{{{\dagger}}} \boldsymbol{\varLambda}^{\mathcal{H}} = \boldsymbol{\mathsf{L}}^{\mathcal{H}} {\boldsymbol{\mathsf{Q}}^{{{\dagger}}}} .\end{equation}
The eigenvalues of the adjoint operator are the complex conjugate of those from direct analysis. The eigenvectors 
$\boldsymbol{\mathsf{Q}}^{{\dagger} }$ represent the region of flow sensitivity and they are equal to 
$\boldsymbol{\mathsf{Q}}^{-1}$ for normal systems. On the other hand, this is not true for non-normal systems. A thorough review on the significance of adjoint operators is presented by Luchini & Bottaro (Reference Luchini and Bottaro2014) and references therein.
In case forcing is applied at a frequency 
$\tilde {\omega }$, (2.1) yields
\begin{equation} \hat{\boldsymbol{q}} = (-\mathrm{i} \tilde{\omega} \boldsymbol{\mathsf{I}} - \boldsymbol{\mathsf{L}} )^{{-}1} \boldsymbol{\mathsf{B}} \hat{\boldsymbol{f}} = \boldsymbol{\mathsf{H}} \, \boldsymbol{\mathsf{B}} \hat{\boldsymbol{f}} , \end{equation}
where the matrix 
$\boldsymbol{\mathsf{H}} = \boldsymbol{\mathsf{H}}(\bar {\boldsymbol {q}},\beta,\tilde {\omega })$ is the resolvent operator (Reddy & Henningson Reference Reddy and Henningson1993; Schmid Reference Schmid2007; McKeon & Sharma Reference McKeon and Sharma2010; Taira et al. Reference Taira, Brunton, Dawson, Rowley, Colonius, McKeon, Schmidt, Gordeyev, Theofilis and Ukeiley2017, Reference Taira, Hemati, Brunton, Sun, Duraisamy, Bagheri, Dawson and Yeh2020). For non-normal operators, eigenvalue sensitivity and energy amplification are related to the induced 
$L_2$ norm of operator 
$\boldsymbol{\mathsf{H}}$, i.e. the leading singular value 
$\sigma$ obtained from singular value decomposition (SVD)
\begin{equation} \boldsymbol{\mathsf{H}} = \boldsymbol{\mathsf{U}} \boldsymbol{\varSigma} \boldsymbol{\mathsf{V}}^{*} .\end{equation}
In this equation, 
$\boldsymbol{\mathsf{V}}$ and 
$\boldsymbol{\mathsf{U}}$ are unitary matrices holding right and left singular vectors and 
$\boldsymbol {\varSigma }$ is a diagonal matrix containing the singular values 
$\sigma$, such that
\begin{equation} \| \boldsymbol{\mathsf{H}} \| = \max( \sigma ) = \sigma_1 . \end{equation}
The first column of 
$\boldsymbol{\mathsf{V}}$ contains the forcing term that produces the largest response in the flow (first column in matrix 
$\boldsymbol{\mathsf{U}}$) with the amplification ratio given by 
$\boldsymbol {\varSigma }$.
A transformation using matrix 
$\boldsymbol{\mathsf{W}}$ is applied to convert variables 
$\hat {\boldsymbol {q}}$ and 
$\hat {\boldsymbol {f}}$ with an appropriate energy norm prior to performing the SVD. In the case of compressible flows, the Chu norm which relates density, velocity and temperature (Chu Reference Chu1965) is typically used. A spatial window 
$\boldsymbol{\mathsf{C}} = \boldsymbol{\mathsf{C}}(\boldsymbol {x})$ is applied to limit the domain of analysis and the system response in Fourier domain 
$\hat {\boldsymbol {y}}$ is given by
\begin{equation} \hat{\boldsymbol{y}} = \boldsymbol{\mathsf{W}} \boldsymbol{\mathsf{C}} \hat{\boldsymbol{q}} .\end{equation}Combining (2.6) and (2.9) leads to the modified resolvent operator
\begin{equation} \boldsymbol{\mathsf{H}}_w = \boldsymbol{\mathsf{W}} \boldsymbol{\mathsf{C}} (-\mathrm{i} \tilde{\omega} \boldsymbol{\mathsf{I}} - \boldsymbol{\mathsf{L}} )^{{-}1} \boldsymbol{\mathsf{B}} \boldsymbol{\mathsf{W}}^{{-}1} .\end{equation}The amplification mechanisms of flow disturbances can be identified by using the eigenvalue decomposition of the LNS operator, (2.4), in the resolvent, yielding
\begin{equation} \boldsymbol{\mathsf{H}}_w = \boldsymbol{\mathsf{W}} \boldsymbol{\mathsf{C}} (-\mathrm{i} \tilde{\omega} \boldsymbol{\mathsf{I}} - \boldsymbol{\mathsf{Q}} \boldsymbol{\varLambda} \boldsymbol{\mathsf{Q}}^{{-}1} )^{{-}1} \boldsymbol{\mathsf{B}} \boldsymbol{\mathsf{W}}^{{-}1} = \boldsymbol{\mathsf{W}} \boldsymbol{\mathsf{C}} \boldsymbol{\mathsf{Q}} ( -\mathrm{i} \tilde{\omega} \boldsymbol{\mathsf{I}} - \boldsymbol{\varLambda} )^{{-}1} \boldsymbol{\mathsf{Q}}^{{-}1} \boldsymbol{\mathsf{B}} \boldsymbol{\mathsf{W}}^{{-}1} ,\end{equation}
where 
$\boldsymbol {\varLambda }$ and 
$\boldsymbol{\mathsf{Q}}$ are the eigenvalues and eigenvectors from the solution of (2.4). The bounds of the induced 
$L_2$ norm are
\begin{equation} \| ( -\mathrm{i} \tilde{\omega} \boldsymbol{\mathsf{I}} - \boldsymbol{\varLambda} )^{{-}1} \| \le \| \boldsymbol{\mathsf{H}}_w \| \le \| ( -\mathrm{i} \tilde{\omega} \boldsymbol{\mathsf{I}} - \boldsymbol{\varLambda} )^{{-}1} \| \| \boldsymbol{\mathsf{Q}} \| \| \boldsymbol{\mathsf{Q}}^{{-}1} \| , \end{equation}
where the lower bound is the case of an operator with orthonormal eigenvectors. In this case, the norm depends only of the resonances, following a 
${1}/{R}$ decay based on the distance 
$R$ in complex plane with respect to the eigenvalues. For non-normal systems, the norm also depends on the pseudoresonances, measured by the product of eigenvectors 
$\boldsymbol{\mathsf{Q}}$ and their inverse 
$\boldsymbol{\mathsf{Q}}^{-1}$. The weighting 
$\boldsymbol{\mathsf{W}}$ and the two windowing matrices 
$\boldsymbol{\mathsf{B}}$ and 
$\boldsymbol{\mathsf{C}}$ are included in the norm calculation (see Schmid & Henningson (Reference Schmid and Henningson2001) and Symon et al. (Reference Symon, Rosenberg, Dawson and McKeon2018) for more details).
To perform the global stability (spectral) and resolvent (pseudospectral) analyses, the discrete linear operator is computed using the second-order accurate methodology from Sun et al. (Reference Sun, Taira, Cattafesta and Ukeiley2017) and Yeh & Taira (Reference Yeh and Taira2019). The base flow is the time–spanwise-averaged solution from the LES, 
$\bar {\boldsymbol {q}} = [\bar {\rho }, \bar {u}, \bar {v}, \bar {w}, \bar {T}]$. For the far-field and airfoil surface, Dirichlet boundary conditions are set for 
$[\rho ' , u' , v' , w' ] = [0, 0, 0, 0]$ and a Neumann boundary condition is set for 
$T'$ such that the wall-normal derivative 
$\partial T' / \partial n = 0$. At the outlet boundary, the same Neumann boundary condition is set for all flow variables. With these boundary conditions and the mean base flow 
$\bar {\boldsymbol {q}}$, the linear operator is computed in its discrete form 
$L(\bar {\boldsymbol {q}},\beta )$ for a prescribed spanwise wavenumber 
$\beta$.
The mean flow 
$\bar {\boldsymbol {q}}$ obtained is interpolated from the O-grid in figure 1(a) to the H-mesh shown in figure 1(b) for the modal analysis. This is necessary to improve spatial accuracy of the linear operator downstream and upstream of the airfoil, where direct and adjoint eigenvectors are supported, respectively (Yeh & Taira Reference Yeh and Taira2019). The present 2-D H-mesh has an extent of 
$x \in [-5,+7], y \in [-5,+5]$ and it is composed of approximately 
$520 \times 10^{3}$ grid points. The mesh refinement targets a Strouhal number cut-off 
$St=8.0$. The sponge region comprises a circle with parabolic growth 
$f = a(r-r_0)^{2}$ for 
$r > r_0$. The sponge is centred at 
$(x, y) = (0.5, 0)$ with 
$a = 1.0$ and 
$r_0 = 1.5$. The mesh refinement and sponge placement are assessed by convergence of physically meaningful eigenvalues while suppressing spurious eigenvalues. If the mesh is coarse or the sponge is placed too far from the airfoil, it was observed that the meaningful eigenvectors also exhibit spatial support in the same region of the spurious eigenvectors.
For the pseudospectrum calculation, the dashed red rectangle in figure 1(b) depicts the region considered for the present resolvent analysis in terms of operators 
$\boldsymbol{\mathsf{B}}$ and 
$\boldsymbol{\mathsf{C}}$ in (2.10). This windowing is important not only from a physical point of view, where the noise sources close to the airfoil surface are most important, but it also eliminates spurious eigenvectors, which are mostly observed far away from the body. Finally, the SVD in (2.10) is performed with the randomised algorithm presented by Ribeiro, Yeh & Taira (Reference Ribeiro, Yeh and Taira2020).
2.3. Mean flow perturbation
The properties of the resolvent operator can be studied as an initial value problem where large transient energy amplification is observed, even in stable problems, due to the non-normality (Trefethen et al. Reference Trefethen, Trefethen, Reddy and Driscoll1993; McKeon & Sharma Reference McKeon and Sharma2010). The time evolution of perturbations with respect to the mean flow can be performed with a linearised version of the CFD code according to (2.1). However, here we employ an approach that uses the present LES code with the addition of a force term to the right-hand side of the equations (Touber & Sandham Reference Touber and Sandham2009; Bhaumik, Gaitonde & Waindim Reference Bhaumik, Gaitonde and Waindim2015; Ranjan, Unnikrishnan & Gaitonde Reference Ranjan, Unnikrishnan and Gaitonde2018). In this approach, known as mean flow perturbation, the Reynolds decomposition 
$\boldsymbol {q} = \bar {\boldsymbol {q}} + \boldsymbol {q}'$ is applied to the nonlinear Navier–Stokes equations which are expressed as
\begin{equation} \frac{\partial (\bar{\boldsymbol{q}} + \boldsymbol{q}')}{\partial t} = \mathcal{N}(\bar{\boldsymbol{q}} + \boldsymbol{q}') . \end{equation} The turbulent mean flow 
$\bar {\boldsymbol {q}}$ is not an equilibrium state of the Navier–Stokes equations and its time derivative 
$\partial _t \bar {\boldsymbol {q}}$ is non-null. Hence, in order to use it as a base flow, an additional term 
$\mathcal {R}$ must be added in the numerical procedure to keep the base flow stationary. A Taylor series expansion about the base state 
$\bar {\boldsymbol {q}}$ yields
\begin{equation} \frac{\partial (\bar{\boldsymbol{q}} + \boldsymbol{q}')}{\partial t} = \mathcal{N}(\bar{\boldsymbol{q}} + \boldsymbol{q}') = \mathcal{N}(\bar{\boldsymbol{q}}) + \frac{\partial \mathcal{N}(\bar{\boldsymbol{q}})}{\partial \boldsymbol{q}} \boldsymbol{q}' + \frac{1}{2} \frac{\partial^{2} \mathcal{N}(\bar{\boldsymbol{q}})}{\partial^{2} \boldsymbol{q}} \boldsymbol{q}'^{2} + \cdots+ \mathcal{R} . \end{equation}Grouping together the base flow terms results in
\begin{equation} \frac{\partial \boldsymbol{q}'}{\partial t} + \left( \frac{\partial \bar{\boldsymbol{q}}}{\partial t} - \mathcal{N}(\bar{\boldsymbol{q}}) - \mathcal{R} \right) = \frac{\partial \mathcal{N}(\bar{\boldsymbol{q}})}{\partial \boldsymbol{q}} \boldsymbol{q}' + \frac{1}{2} \frac{\partial^{2} \mathcal{N}(\bar{\boldsymbol{q}})}{\partial^{2} \boldsymbol{q}} \boldsymbol{q}'^{2} + \cdots, \end{equation}where the base flow is stationary when
\begin{equation} \left( \frac{\partial \bar{\boldsymbol{q}}}{\partial t} - \mathcal{N}(\bar{\boldsymbol{q}}) - \mathcal{R} \right) = 0 .\end{equation} Considering that the imposed perturbation amplitude 
$\boldsymbol {q}'$ is sufficiently small leads to a linearised solution of the Navier–Stokes equations as
\begin{equation} \frac{\partial \boldsymbol{q}'}{\partial t} \approx \frac{\partial \mathcal{N}(\bar{\boldsymbol{q}})}{\partial \boldsymbol{q}} \boldsymbol{q}' = \boldsymbol{\mathsf{L}} \boldsymbol{q}' , \end{equation}
where the Jacobian 
${\partial \mathcal {N}(\bar {\boldsymbol {q}})}/{\partial \boldsymbol {q}}$ is the linearised operator 
$\boldsymbol{\mathsf{L}}$. In contrast to the previous approach in § 2.2, this methodology is matrix-free and the linear operator is never explicitly computed. Despite this, the method still accounts for the modal interaction that leads to transient energy amplification within a non-normal analysis of the linear operator. The system response is obtained by the time evolution of the disturbances and the initial condition is set as an impulsive excitation that triggers the dominant response. In the current analysis, the grid and numerical methodology employed are the same used in the LES.
3. Results
The physical mechanisms responsible for the generation of secondary tones in airfoil flows are investigated by post-processing the nonlinear results from LES and also employing linear stability theory. The linear analysis provides insights on the onset and growth of disturbances in the flow in addition to its receptivity. On the other hand, LES results allow the investigation of transition, intermittency and phase interference effects in the context of the multiple secondary tones and the acoustic feedback loop.
3.1. Mean flow
The mean flow is presented with contours of 
$u$-velocity normalised by the freestream speed of sound in figure 2. Important flow features can be observed including a long recirculation bubble that extends over the suction side, from the detachment at 
$x = 0.34$ until its reattachment at 
$x = 0.82$. In the magnified view, a small bubble can be observed on the pressure side, followed by a recirculation region at the trailing edge. The blue contours enclosed by magenta lines indicate regions of reversed flow.
Contours of mean streamwise velocity 
$\bar {u}$ normalised by freestream speed of sound. The magenta dashed lines depict the reversed flow boundaries which include a wide separation bubble on the suction side besides a small bubble on the pressure side, near the trailing edge.
Distributions of root-mean-square (RMS) for kinetic energy 
$k$ and pressure 
$p$ are presented in figures 3(a) and 3(b), respectively. As it can be seen from the plots, fluctuations start amplifying along the bubble and the highest values appear at the reattachment region, on the suction side. The green and blue dashed lines in figures 3(a) and 3(b), respectively, depict the locations of maximum 
$k$ and 
$p$ along the shear layer forming on the suction side. In the following analyses, these lines will be used as reference locations for data extraction to track flow disturbances. Again, the solid magenta line shows the recirculation bubble.
RMS values of (a) kinetic energy and (b) pressure. The green and blue dashed lines depict locations of maximum fluctuations along the shear layer on the suction side, whereas the magenta solid line delimits the reversed flow region.
The negative mean pressure coefficient 
$-C_p$ is shown in figure 4(a) with blue and red solid lines for suction and pressure sides, respectively. To highlight the regions with intense fluctuations, the pressure coefficient is presented as 
$-C_p \mp C_{p_{{RMS}}}$, where the blue dashed and red dotted lines correspond to the suction and pressure sides, respectively. The peak value of pressure coefficient on the suction side is observed near the leading edge at 
$x = 0.02$. From this point onward, a pressure increase leads to a drop in 
$-C_p$ until reaching a plateau along the LSB. The plateau extends up to 
$x = 0.7$, where the pressure further increases towards the trailing edge, causing another drop in 
$-C_p$. Pressure fluctuations are shown in figure 4(d) for both suction and pressure sides. The RMS values of 
$C_p$ are low along the entire pressure side, and upstream of the bubble (
$x \lesssim 0.5$) on the suction side. However, they increase towards the peak at 
$x \approx 0.78$ on the suction side.
Mean and RMS distributions of (a,d) pressure and (b,e) friction coefficients along the airfoil surface, and (c, f) tangential velocity profiles computed in the wall-normal direction 
$\Delta n$.
The mean skin friction 
$C_f$ distribution over the airfoil is presented in figure 4(b). The RMS values of skin friction 
$C_{f_{{RMS}}}$, shown in figure 4(e), are also added to the mean values. A magnified view highlights the presence of the LSB which starts at 
$x = 0.34$ and reattaches at 
$x = 0.67$. The reattachment is followed by a second detachment which reaches a minimum at 
$x = 0.75$ and reattaches at 
$x = 0.82$. The nature of the double detachment profile is discussed in details by Duck et al. (Reference Duck, Ruban, Theofilis, Hein and Dallmann2000). The RMS peaks at 
$x \approx 0.76$ due to the high value of 
$k$ shown in figure 3(a).
Mean tangential velocity profiles extracted along the wall-normal direction 
$\Delta n$ are presented in figure 4(c) at different locations along the chord. At the detachment point, 
$x = 0.34$, the red curve indicates the zero wall-normal derivative in the velocity profile. At 
$x = 0.73$, the velocity profile shown as a grey line exhibits the maximum reversed flow of 
$-0.13 U_\infty$. The orange line shows a velocity profile at 
$x = 0.95$, where 3-D effects are important due to turbulent transition. In figure 4( f), RMS values of the tangential velocity are presented and a single peak is observed for the profile computed at 
$x = 0.34$. Velocity fluctuations increase inside the bubble at 
$x = 0.64$ and become dominant at 
$x = 0.73$. In both these locations, the velocity fluctuations present triple peak profiles as also observed by Nash et al. (Reference Nash, Lowson and McAlpine1999), Desquesnes et al. (Reference Desquesnes, Terracol and Sagaut2007) and Duck et al. (Reference Duck, Ruban, Theofilis, Hein and Dallmann2000). Near the trailing edge at 
$x = 0.95$, where a turbulent regime may occur, mixing results in a smoother velocity fluctuation profile.
3.2. Vortex dynamics
Instantaneous 3-D flow fields over the airfoil are presented in figure 5 with iso-surfaces of Q-criterion coloured by 
$u$-velocity. A magenta shade along the airfoil surface depicts the region of reversed flow along the separation bubble. This figure shows 2-D rolls that are visible along the separation bubble and amplify on the suction side, leading to vortex shedding. Moreover, different flow regimes are observed at the trailing edge, where coherent structures alternate with periods of turbulent packets. For example, figure 5(a) shows a single laminar-like roll reaching the trailing edge at 
$t = 20.82$ whereas figure 5(b) shows that, at 
$t = 21.78$, vortex breakdown leads to a turbulent regime near the trailing edge. Readers are referred to the movie of the flow field provided as supplementary movie 1 are available at https://doi.org/10.1017/jfm.2022.129 for detailed inspection of the present flow dynamics.
Instantaneous iso-surfaces of Q-criterion coloured by 
$u$-velocity show regimes with (a) coherent structures at the trailing edge (high spanwise coherence, 
$t=20.82$) and (b) smaller-scale turbulent eddies (low spanwise coherence, 
$t=21.78$).
The analysis of flow snapshots is important to understand the vortex dynamics over the airfoil suction side. Based on the acoustic analogy of Ffowcs-Williams & Hall (Reference Ffowcs-Williams and Hall1970), efficient sound generation is achieved when pressure fluctuations are aligned with an edge, such that it will be maximum for 2-D-like perturbations. In this sense, figure 6 shows the spanwise-averaged 
$\omega _z$ vorticity, where blue and red contours represent negative and positive values, respectively. A dark blue colour indicates strong spanwise coherence of a 2-D vortex, whereas light blue contours are representative of uncorrelated turbulence. A thin black line is used to identify individual vorticity packets, based on the iso-contour of an entropy measure given by 
${p}/{\rho ^{\gamma }} - ({p_{\infty }}/{\rho _{\infty }^{\gamma }}) = 0.1\,\%$. This separates the region where vorticity is important from where the flow is irrotational, outside the boundary layer. A magenta dashed line represents the border of averaged negative 
$\bar {u}$ velocity, as discussed in figure 2. The simulation time is shown in convective time units on the lower left corner of all plots. A movie of figure 6 is provided as supplementary movie 2 to aid the dynamic visualisation of the flow features.
Spanwise-averaged vorticity 
$\omega _z$ shows vortex pairing over the airfoil. The magenta dashed line represents the boundary of reversed flow (negative 
$\bar {u}$). The black lines connecting the sub-figures mark the evolution of vortices with high (
$\boldsymbol {-}$ and 
${-}\,{-}$) or low (
$\cdots$ and 
${-}\cdot {-}$) spanwise coherence.
As shown in figure 6, the flow dynamics is dominated by events on the suction side, where vortices are shed from the LSB. During the shedding process, there is a possibility that the 2-D laminar vortices undergo a process of vortex pairing. In case this process is not successful, vortex bursting leads to turbulent transition and low coherence along the span, represented by dotted lines. These turbulent packets may interact with other vortices and further reduce their coherence, as indicated by the dashed-dotted line. On the other hand, if the vortices merge or only a solitary vortex is shed from the bubble, the trend is for either structure to keep its high coherence up to the trailing edge. Both of these processes are represented by the dashed (merging) and solid (single vortex) lines in figure 6. A detailed discussion on the consequence and cause of this behaviour is presented in §§ 3.3 and 3.6, respectively. Similar dynamics of vortex shedding from LSBs are observed experimentally by Kurelek et al. (Reference Kurelek, Lambert and Yarusevych2016).
As shown in the experimental investigations of Pröbsting & Yarusevych (Reference Pröbsting and Yarusevych2015) and Pröbsting et al. (Reference Pröbsting, Scarano and Morris2015), for this Reynolds number and angle of attack, it is expected that the pressure side does not play an important role in terms of the flow dynamics. In order to confirm this observation, a temporal signal of pressure fluctuation is extracted at 
$x = 0.98$, marked by the blue dashed line in figure 3(b). A second signal is obtained for the velocity derivative in the wall-normal direction on the pressure side at 
$x = 0.95$. Both signals are presented in figure 7(a) as well as the time-averaged 
${\rm d} \bar {u}/{\rm d} n$, indicating that the flow is, on average, recirculating (
${\rm d} \bar {u}/{\rm d} n<0$). In the velocity derivatives, the wall-normal direction is positive pointing outwards from the airfoil surface. Snapshots of the flow dynamics are presented in 7(b), where the red–blue contours are the 
$\omega _z$ vorticity, similarly to the previous figures, the green hatched colours show instantaneous regions of negative 
$u$-velocity whereas the dashed magenta line indicates the boundaries of the time-averaged pressure side bubble. In this figure, it is possible to see that when a vortex arrive at the trailing edge from the suction side at 
$t = 20.59$ (negative 
$p'$ peak), the flow reattaches on the bottom side (positive 
${\rm d} u/{\rm d} n$ at the wall). As the vortex leaves the airfoil surface on the suction side at 
$t=20.73$, the pressure increases on the suction side (positive 
$p'$) and, as a consequence, the flow detaches on the pressure side (negative 
${\rm d} u/{\rm d} n$). Thus, the magenta dashed line in figure 7(a) shows that the recirculation bubble exists on an time-averaged sense. Moreover, these figures show that the pressure side dynamics seems to be a reaction to the suction side events in the present flow.
Intermittency of the pressure side bubble exhibits dependence on the suction side vortices as presented for (a) temporal signals and (b) flow snapshots. In (a), the pressure signal is extracted on the suction side at 
$x = 0.98$ and the velocity derivative is computed on the pressure side at 
$x = 0.95$. In (b), the hatched green region indicates instantaneous values of negative 
$u$-velocity. In both plots, the magenta dashed line indicates negative values in the time-averaged 
$u$-velocity.
3.3. Intermittency and sound generation
The acoustic analogy of Ffowcs-Williams & Hall (Reference Ffowcs-Williams and Hall1970) states that efficient acoustic scattering occurs due to flow fluctuations near a trailing edge. In this sense, the hydrodynamic near field with pressure fluctuations is presented in figures 8(a) and 8(c) with black–white contours, while the acoustic field is presented in figures 8(b) and 8(d). In the latter, the levels are presented 10 times lower compared with the former. The figures also display 
$\omega _z$ vorticity in red–blue colours. To highlight the spanwise coherence, regions of 
$|\omega _z| > 24$ are plotted as filled contours. The vortex cores are related to lower pressures (white background contours) whereas the gaps between the cores display higher pressures (black contours). A coherent vortex at the trailing edge is presented in figure 8(a) at a time instant 
$t=20.82$. The highly coherent hydrodynamic structure leads to efficient acoustic scattering at the trailing edge which, in turn, induces the emission of an intense acoustic wave. The wave is scattered with a 
${\rm \pi}$ phase opposition relative to the incident fluctuations and it reaches the leading edge at a retarded time 
$t = 21.31$, as shown in figure 8(b). When the flow is turbulent, at 
$t = 21.78$, a less-coherent vorticity packet reaches the trailing edge as shown in figure 8(c). In this case, a weaker acoustic wave is emitted from the trailing edge, reaching the leading edge at a retarded time of 
$t = 22.27$, as exhibited in figure 8(d). A movie of figure 8 is submitted as supplementary movie 3 to show the relation between the acoustic field and vortex dynamics.
Acoustic radiation from the airfoil. (a) A coherent structure near the trailing edge at 
$t = 20.82$ and (b) its subsequent intense acoustic wave reaching the leading edge, at a retarded time 
$t = 21.31$. (c) An uncorrelated turbulent packet near the trailing edge at 
$t = 21.78$ and (d) its subsequent weaker acoustic wave reaching the leading edge, at a retarded time 
$t = 22.27$. Note that the levels from (b) and (d) are 10 times lower than those from figures (a) and (c).
Temporal signals are acquired at different positions along the blue dashed line in figure 3(b) in terms of spanwise-averaged non-dimensional pressure. The signal at the trailing edge is presented in figure 9(a) and displays intermittent deep valleys related to the passage of coherent structures as shown in figure 8(a). The signal also shows low-amplitude oscillations related to the turbulent packets that carry a negative pressure envelope, as described in figure 8(b). Such behaviour can be characterised as an amplitude modulation of the signal. The same trends are also observed in the acoustic pressure signal, which is extracted one chord above the trailing edge, at 
$(x,y) = (1c,1c)$, and presented in figure 9(c). It is possible to see that the low-pressure valleys from figure 9(a) result in high-pressure peaks that reach the observer position at a retarded time.
Time-signal and Fourier transform of spanwise-averaged pressure at different locations (a) and (b) inside the boundary layer along the chord and (c) and (d) in the acoustic field.
The magnitude of the Fourier coefficients is presented in figure 9(b) and 9(d) for the hydrodynamic and acoustic fields, respectively, displaying equidistant multiple tones. The Fourier transforms employ a frequency domain averaging with an overlap of 
$67\,\%$ amongst 4 bins covering the entire period. The frequency resolution is 
$\Delta f = 0.03$ and a Hanning window function with energy correction is applied. All tonal peaks are integer multiples of the lowest-frequency tone, 
$St \approx 0.5$. As it can be seen in figure 9(b), the dominant peak at the trailing edge is at 
$St \approx 3.5$ for probes located close to the leading edge (
$x = 0.05$), at the detachment point (
$x = 0.34$) and the trailing edge (
$x = 0.98$). This frequency is related to the passage of low-pressure packets, either turbulent or coherent. At 
$x = 0.64$, in the shear-layer roll-up region, the dominant peak moves to higher frequencies which may be related to the smaller-scale instabilities yet to be merged into larger vortices. In the acoustic field, figure 9(d) presents results extracted at one and five chords above the trailing edge. The same frequencies are observed where 
$St = 3.5$ is still the most prominent peak.
As can be seen from figure 9(a), coherent structures often reach the trailing edge with a slow timescale of 
$\Delta t \approx 2.0$ convective time units, which corresponds to a frequency 
$St = 0.5$. However, due to the chaotic motion, this process is not perfectly periodic and multiple vortices may reach the trailing edge within this period, as indicated by the time differences 
$\Delta t = 0.22$, 0.25, 0.29, 0.33 and 0.40 in the figure. In this faster time scales, the different intervals are related to frequencies 
$St = 4.5$, 4.0, 3.5, 3.0 and 2.5, respectively, which are integer multiples of the low-frequency tone at 
$St = 0.5$.
To better highlight this behaviour, a time-frequency analysis is also performed with the continuous wavelet transform (CWT) using the complex Morlet wavelet (Farge Reference Farge1992), as shown in figure 10. The present wavelet function is the product of a monochromatic complex exponential with a Gaussian envelope. Its spectral response is an exponentially decaying band-pass filter centred around the frequency of the complex exponential. The filter width depends on the Gaussian standard deviation. The scalogram presents how the correlation of the wavelet with a reference signal changes in time based on the wavelet frequency. Hence, higher values indicate strong fluctuations of the signal at a particular frequency and time. Results are obtained for the same probe locations considered in figure 9(b). The scalograms from figures 10(a) and 10(c) exhibit longer periods of intense pressure fluctuations at 
$St = 3.5$, whereas in figure 10(b), a broader range of frequencies is excited more intensely. As it can be seen from the probe at 
$x=0.64$, frequencies higher than 
$St = 3.5$ are excited, particularly 
$St = 4.5$, 
$5.0$ and 
$5.5$. One should note that the colour levels are different for each location, with lower and higher amplitudes for probes at 
$x=0.34$ and 
$0.98$, respectively. These figures show that the secondary tones present intermittent behaviour, with each frequency experiencing periods of silence (blue contours) and loud noise emission (green, yellow and red contours). Another important aspect captured by the scalograms is the amplitude modulation of the signals. Such behaviour has an important role in the generation of secondary tones, as discussed by Desquesnes et al. (Reference Desquesnes, Terracol and Sagaut2007), Pröbsting et al. (Reference Pröbsting, Serpieri and Scarano2014) and Ricciardi et al. (Reference Ricciardi, Arias-Ramirez and Wolf2020).
Time-frequency analysis of pressure fluctuations using a CWT at (a) the detachment point 
$x = 0.34$, (b) shear-layer roll-up 
$x = 0.64$ and (c) trailing edge 
$x = 0.98$.
3.4. Linear analysis and amplification mechanisms
In order to investigate the growth of velocity and pressure fluctuations at the tonal frequencies, figure 11 presents Fourier transforms of the temporal signals extracted along the airfoil suction side, following the green (velocities) and blue (pressure) dashed lines from figures 3(a) and 3(b), respectively. In this analysis, the temporal signals are computed for the spanwise-averaged normal and tangential velocity components 
$u_n$ and 
$u_t$, respectively, in addition to pressure 
$p$. Results are plotted together with the RMS values (solid black line) extracted from figure 3. Both 
$\hat {u}_t$ and 
$\hat {p}$ show strong amplification downstream of 
$x = 0.5$, which coincides with the location of the recirculation bubble. Furthermore, a plateau is observed upstream of 
$x = 0.5$ for pressure fluctuations, which indicates that acoustic waves dominate over hydrodynamic disturbances along this region. Pressure fluctuations at 
$St \approx 3.5$ exhibit higher amplitudes along the entire chord and similar observations can be made for the velocity components. As discussed by Kurelek et al. (Reference Kurelek, Yarusevych and Kotsonis2019), it is possible that these fluctuations act as a sub-harmonic forcing that stimulates the vortex merging process.
Magnitude of Fourier transform coefficients for (a) normal velocity 
$\hat {u}_n (St,x)$, (b) tangential velocity 
$\hat {u}_t (St,x)$ and (c) pressure 
$\hat {p} (St,x)$ computed along the green (velocities) and blue (pressure) lines from figures 3(a) and 3(b), respectively. Each colour represents a tonal frequency, identified in the colourbar, whereas the RMS is presented as the black line.
The mechanisms of generation and amplification of flow instabilities is explained by linear stability theory. In this regard, bi-global linear stability and resolvent analyses are employed to understand the most sensitive frequencies in terms of optimal forcing and response characteristics for 2-D perturbations, i.e. only for spanwise wavenumber 
$\beta = 0$. This is justified because 2-D structures are expected to radiate noise more efficiently (Sano et al. Reference Sano, Abreu, Cavalieri and Wolf2019). In the LNS operator, the off-diagonal terms, related to velocity gradients and shear, lead to a non-normal operator. Hence, orthogonality of the stability modes is not expected and interaction among modes within the system may lead to transient response that affect its short-term dynamics (Schmid Reference Schmid2007). This behaviour is properly captured by the resolvent analysis, which provides the information on eigenvalue sensitivity and transient energy amplification rate.
Following the theoretical and numerical methodology described in § 2.2, the eigenvalue spectrum from a bi-global stability analysis is presented in figure 12. The black circles represent the physically meaningful eigenvalues whereas the empty grey symbols are spurious eigenvalues. The latter could be identified due to their large displacement in the complex plane during the grid/sponge convergence study (not shown) as well as by visualisation of the associated eigenvectors (modes). Sensitivity of the operator to external disturbances is computed by the resolvent analysis along both real (frequency) and imaginary (growth rate) axes and it is measured by the pseudospectrum, presented as the contour plot in figure 12. The region of analysis in the resolvent operator is limited by means of matrices 
$\boldsymbol{\mathsf{B}}$ and 
$\boldsymbol{\mathsf{C}}$ in (2.10). Only the near-field grid points are considered, i.e. those inside the dashed red rectangle in figure 1. This is physically justified because the most intense quadrupole noise sources which are responsible for the incident field in the acoustic scattering process are close to the trailing edge (Wolf et al. Reference Wolf, Azevedo and Lele2012a). As an additional positive side effect, this improved results by reducing spatial support of spurious eigenvectors.
Spectrum (symbols) and pseudospectrum (contours in logarithmic scale) of the linearised compressible Navier–Stokes equations for 2-D spanwise perturbations 
$(\beta = 0)$. Higher frequencies present larger sensitivity as shown in light blue and green colours.
Results from the bi-global stability analysis in terms of the eigenspectrum are presented in figure 12. It is possible to see multiple eigenvalues at frequencies close to those from the LES, previously discussed in figure 9(b,d). For instance, the maximum relative deviation is less than 5 %, at 
$St \approx 2.0$, and around 2 % for the most unstable frequency 
$St \approx 5.0$. Thus, the vortex dynamics and the multiple tones seem to be triggered by linear mechanisms. The pseudospectrum, obtained from the resolvent analysis, is shown as a contour plot in logarithmic scale of the leading singular value from the SVD according to (2.7). The input–output amplification ratio depends on both resonance and pseudoresonance (see (2.12)) such that peaks appear when approaching an eigenvalue. The pseudoresonance levels depend on the mutual excitation of eigenvectors with more sensitive eigenvalues presenting higher peaks. In this sense, it is possible to see in figure 12 that higher frequencies are more sensitive compared with the lower frequencies. Moreover, the stable eigenvalues may affect the flow dynamics in terms of transient energy amplification given their sensitivities to disturbances.
The most unstable frequencies of 
$St \approx 4.5$, 
$5.0$ and 
$5.5$ are observed in the LES on the shear-layer roll-up region at 
$x=0.64$, as illustrated by figure 10(b). The interaction among multiple frequencies of the flow instabilities and the vortex merging reduce the shedding frequency in the LES, where the dominant tonal peak becomes 
$St \approx 3.5$. This nonlinear phenomenon is captured by the LES but not by the linear analysis. A more detailed discussion regarding the necessary conditions for vortex merging is presented in § 3.6.
The modal shape of each frequency is crucial to understand the system dynamics and to identify amplification regions within the flow. In this sense, figure 13(a)–13(c) present the real part of the global stability modes for the most unstable frequency 
$St = 5.08$ in terms of 
$u'$, 
$v'$ and 
$p'$, respectively. The modes are normalised with respect to the maximum value for each variable and the magnitude doubles every level to better represent lower values. The modal shapes computed for other relevant frequencies in the spectrum (resonances) show similar patterns compared with that described here and differences are only related to their characteristic wavelengths. Along the separation bubble, the modal structures of 
$u'$ and 
$v'$ are tilted against the mean shear, suggesting an Orr mechanism for transient growth of energy and amplification of disturbances (Schmid & Henningson Reference Schmid and Henningson2001). After flow reattachment, the structures become more aligned with the wall-normal direction. The pressure disturbances 
$p'$ display upstream travelling acoustic waves radiated from the trailing edge due to scattering mechanism. The modes displayed in the figure originate as shear layer instabilities on the airfoil suction side and extends along the wake. Hence, wake–boundary layer coupling mechanisms cannot be discarded.
Global stability modes of the most unstable frequency 
$St = 5.08$ for primitive variables (a) 
$u'$, (b) 
$v'$ and (c) 
$p'$.
The dominant resolvent modes are computed along the neutral stability axis 
$(\omega _i = 0.0)$ for 
$St = 5.08$ in terms of kinetic energy 
$k$ and presented in figure 14 normalised with respect to their maximum value. The response mode, shown in white–purple contours, depicts the location where disturbances are more energetically relevant. Here, its modal shape magnitude is presented such that it doubles every two levels. It is possible to see negligible levels from the leading edge until the recirculation bubble. Then, for 
$x > 0.4$, its magnitude increases quickly, indicating that the suction side bubble is an amplifier of disturbances. The forcing mode, upper-left plot with white–blue contours, highlights the most receptive location, where minimal disturbances result in larger transient growth by the system. This mode is related to the adjoint operator of the LNS equations and it is introduced in the context of airfoil secondary tones by Fosas de Pando et al. (Reference Fosas de Pando, Schmid and Sipp2017). In the present results, the leading-edge region on the suction side shows the higher values of forcing modes. Furthermore, these modes are absent on the bottom side which corroborates to the absence of pressure side-driven events for this flow configuration.
Response (white–purple contours) and forcing (upper-left plot with white–blue contours) modes of kinetic energy 
$k$ from the resolvent analysis along the neutral stability axis for 
$St = 5.08$.
3.5. Acoustic feedback loop mechanism
As discussed by several authors (Tam Reference Tam1974; Lowson et al. Reference Lowson, Fiddes and Nash1994; Desquesnes et al. Reference Desquesnes, Terracol and Sagaut2007; Jones & Sandberg Reference Jones and Sandberg2011; Fosas de Pando et al. Reference Fosas de Pando, Schmid and Sipp2014b), airfoil flows at transitional Reynolds numbers are subject to an acoustic feedback loop mechanism that is related to the multiple tones. In this feedback process, hydrodynamic perturbations are excited in the vicinity of the leading edge and convected downstream. The disturbances amplify along the recirculation bubble on the airfoil surface and, upon reaching the trailing edge, scatter acoustic waves that travel upstream and trigger new disturbances. In this context, figure 15(a)–15(h) present results from the resolvent analysis for different frequencies. The real part of the pressure forcing modes is shown in red–blue contours whereas black–white contours display the real component of the leading pressure response modes. The modes have been normalised with respect to their maximum values. To highlight the most receptive locations for each individual frequency, the yellow spots in the forcing modes are shown only for regions where the magnitude reaches higher values than 98 % of the maximum. As presented in figure 15, the locations of maximum receptivity change with frequency and vary along the interval 
$0.10 < x < 0.18$.
Real component of pressure response (black–white contours) and forcing (blue–red contours) modes from resolvent analysis along the neutral stability axis. The yellow spots represent the most sensitive regions where the forcing modes achieve at least 98 % of the respective maximum values for each frequency.
Aspects such as the amplitude and phase between the acoustic waves and flow receptivity at different frequencies impact the onset of hydrodynamic disturbances. To understand such complex interrelation, the acoustic feedback is analysed using the mean flow perturbation technique to solve the LNS equations around the 3-D mean flow in time domain (Jones & Sandberg Reference Jones and Sandberg2011; Fosas de Pando et al. Reference Fosas de Pando, Schmid and Sipp2014b). The method assumes that the base flow is steady and that it is not modified by the linear perturbations. The dynamics of the pressure fluctuations is presented in figure 16. For the initial condition, a pressure disturbance is applied at the leading edge, from 
$x = 0.0$ to 
$0.005$ over a short duration of 
$\Delta t = 1 \times 10^{-4}$ to excite all frequencies (impulse response). This triggers a wavepacket that is advected by the mean flow along the suction side. Figures 16(a) and 16(b) present the growth of disturbances for 
$x > 0.4$, where spatial amplification described by the response modes becomes important, as shown in figure 14. When the wavepacket reaches the trailing edge, acoustic waves are generated and travel upstream as shown in figure 16(c). After the wavepacket reaches the wake, hydrodynamic disturbances are no longer observed on the suction side, as displayed by figure 16(d). Afterwards, disturbances in the wake slowly decay whereas a new wavepacket amplifies on the suction side (figure 16e) and ( f). Similarly to the low frequency observed in the LES calculation, the wavepackets have a period of 
$\Delta t = 2.0$ and the next cycle is presented in figure 16(g)–16( j). The combined response from the suction side wavepacket, wake instability and upstream travelling acoustic waves keeps the fluctuation field active, amplifying the globally unstable frequencies in the eigenspectrum of figure 12, which ultimately become dominant. Readers are referred to the movie provided as supplementary movie 4 presenting the evolution of the linearised flow dynamics.
Snapshots of the LNS equations presented in terms of pressure fluctuations show periodic wavepackets composed by the multiple frequencies from the eigenspectrum of figure 12.
The feedback loop mechanism is investigated here by analysing the generation of flow instabilities by upstream propagating acoustic waves, as presented in figure 17. In this figure, the entropy measure 
${p}/{\rho ^{\gamma }}$ is given by red–blue contours whereas pressure fluctuations are given by blue–green–red lines, where solid and dashed lines are positive and negative values, respectively. The entropy measure is used to partially filter out fluctuations associated to acoustic waves whereas pressure is used to highlight acoustic waves. In figure 17(a), it is possible to see an incoming acoustic wave at 
$t = 5.888$ that, upon reaching the leading edge at 
$t = 5.920$, triggers an entropy fluctuation due to the secondary diffraction. This latter phenomenon is typical of airfoil noise as discussed by Miotto, Wolf & de Santana (Reference Miotto, Wolf and de Santana2017). The entropy fluctuation at the leading edge, presented in red colour and highlighted by a dotted line box, is propagated downstream in figure 17(b–d) at 
$t = 5.920$, 
$5.952$ and 
$5.984$, respectively. This disturbance overlaps with an incoming entropy wave moving upstream, carried by an acoustic wave depicted with blue dashed lines. Then, in figure 17(e) and 17( f) at 
$t = 6.016$ and 
$6.048$, respectively, the downstream propagating entropy fluctuation encounters a second positive acoustic wave in the region where the forcing modes are maximum, i.e. for 
$0.10 < x < 0.18$. This upstream propagating acoustic wave carries a positive entropy fluctuation which impinges on that generated at the leading edge. Thus, although the onset of instabilities occurs very close to the leading edge, the interaction between upstream travelling acoustic waves and downstream propagating boundary layer instabilities cannot be discarded because it occurs along the entire sensitive region of the flow. Albeit the analysis is conducted for pressure and entropy, local velocity fluctuations also occur and contribute to the incoming perturbations seen by the bubble. This overall process of the feedback loop mechanism is better visualised in supplementary movie 5.
Snapshots of the LNS equations presented in terms of entropy (contour) and pressure (lines) fluctuations. A downstream propagating entropy wave (highlighted by dotted line box) is generated due to leading-edge secondary diffraction. This wave impinges an upstream propagating entropy fluctuation carried by an acoustic wave at the region of maximum receptivity.
Flow disturbances are tracked along the black dashed line in figure 17 (also shown as a blue dashed line in figure 3b) and presented as a function of time and space in figures 18(a) and 18(b). The former shows pressure fluctuations 
$p'$ whereas the latter displays disturbances of an entropy measure 
$({p}/{\rho ^{\gamma }})'$. In these plots, the values are normalised with respect to the maximum and presented in a logarithmic scale. The positive and negative values are presented as red and blue colours, respectively. The first wavepacket shown in figure 16(a,b) is observed travelling downstream with a local mean velocity 
$\bar {u}$ as bottom-up contours starting at 
$t = 0.0$. Upon reaching the trailing edge, at 
$t > 1.2$, the wavepacket generates acoustic waves propagating upstream with 
$\bar {u} - a$ speed. These waves are observed as top-down contours with lower magnitudes. The acoustic pressure has a phase opposition of 
${\rm \pi}$ compared with the hydrodynamic fluctuations on the trailing edge and, hence, an incident negative hydrodynamic pressure fluctuation (blue) creates a scattered positive acoustic wave (red).
Space–time map of linear solution obtained by the mean flow perturbation for (a) pressure and (b) entropy. The map is normalised with respect to the maximum value and presented in a logarithmic scale. The positive and negative fluctuations are indicated by red and blue colours, respectively.
Upstream of 
$x \approx 0.4$, the acoustic waves dominate and it is not possible to see the onset of instabilities. Hence, a measure of entropy 
$({p}/{\rho ^{\gamma }})'$ is shown in figure 18(b) to partially filter the upstream acoustic disturbances and highlight the hydrodynamic content. In this figure, the contours of hydrodynamic fluctuations from bottom-up are better visualised. In both plots, a green region covering the entire period is placed at 
$0.1 \le x \le 0.18$, illustrating the maximum forcing region from figure 15. Along this region, disturbances generated in the vicinity of the leading edge (
$x \approx 0.0$) travel downstream and appear to interact with the incoming acoustic waves, as depicted in figure 17(e) and 17( f). Such interaction results in a chessboard-like interference pattern where the entropy levels carried by upstream-propagating acoustic waves (inside the boundary layer) are of the same order of the entropy disturbances generated at the leading edge. We conjecture that the simultaneous combination of acoustic and hydrodynamic interactions at this location results in the maximum forcing. Finally, black lines are drawn on the plots using the convective velocity information for both pressure and entropy disturbances. These lines are used to track the feedback dynamics and show that, if the feedback indeed closes along the most receptive region, the period of 
$\Delta t \approx 2.0$ is satisfied.
The analysis of the feedback mechanism for the LES results is presented based on the RMS values of either the pressure 
$p$ or the entropy measure 
${p}/{\rho ^{\gamma }}$ averaged along the 
$z$-direction (spanwise direction). Similarly to the linear analysis, time signals are extracted along the blue line in figure 3(b) and the RMS values are presented in a logarithmic scale in figure 19(a) for pressure and figure 19(b) for the entropy measure. In both cases, it is possible to see high fluctuation values for 
$x > 0.4$ travelling downstream towards the trailing edge. For 
$x < 0.4$, smaller values of 
$p'$ propagate upstream towards the leading edge. Similarly to the linear analysis, the former and latter disturbances represent hydrodynamic and acoustic waves, respectively. The entropy measure partially filters the acoustic waves moving upstream. Based on linear analysis results, the onset of disturbances occurs very close to the leading edge in a position that coincides with the small hump observed in figure 4 for pressure and friction coefficients in terms of RMS values 
$C_{p_{{RMS}}}$ and 
$C_{f_{{RMS}}}$, respectively. Thus, the wave secondary diffraction at the leading edge excites boundary layer disturbances which are also observed in the LES. The entropy fluctuations generated at the leading edge appear to interact with the upstream-propagating acoustic waves in the region of maximum receptivity, similarly to the phenomenon observed in the linear analysis. In figure 19, black lines are used to track disturbances assuming that the feedback mechanism takes place at 
$0.1 < x < 0.18$, indicated by the green shaded rectangle which represents the region of maximum receptivity from figure 15.
RMS values along the spanwise direction computed as a function of chordwise position 
$x$ and time 
$t$ for (a) pressure 
$p'$ and (b) entropy 
$(p/p^{\gamma })'$.
As a final remark, the repeating patterns observed in the simulations demonstrate that acoustic feedback occurs in the present flow. Results indicate that the loop depends on the leading-edge secondary diffraction and the region of maximum sensitivity, located downstream of the leading edge and upstream of the separation bubble. However, the intermittent interaction of vortex merging and bursting on the suction side may eventually disrupt the cycle by altering the phase and magnitude of the tonal frequencies with a temporary impact on the feedback loop.
3.6. On vortex merging and phase interference
From the previous sections, it is observed that the acoustic feedback loop is self-sustained by strong acoustic waves scattered from coherent flow structures advected past the trailing edge. These upstream-propagating acoustic waves create new flow instabilities which are amplified by the separation bubble and undergo a process of vortex pairing. Thus, it is important to understand which circumstances lead to generation of coherent structures or small-scale eddies, as discussed in § 3.3. For this task, the dynamics of wavepackets in the linear simulation is first investigated, followed by a similar analysis in terms of the LES data.
The amplification of the wavepacket along the recirculation bubble results in a maximum value of the perturbations evaluated by the mean flow perturbation at 
$x = 0.72$. The temporal signal of 
$p'$ at this location is presented in figure 20(a) depicting three envelopes. These are evaluated using the Hilbert transform to obtain the imaginary part and, thus, the magnitude of a real signal. Then, using the CWT, it is possible to extract the phase information of the previous signal. To highlight the relevant periods of the dynamics, the phase is weighted by the magnitude of the CWT coefficients as presented in figure 20(b). Comparing the temporal signal with the weighted phase, it is possible to see that the wavepacket envelopes peak when there is phase alignment across the multiple frequencies of the eigenspectrum shown in figure 12. Hence, constructive interference of multiple modes is crucial for maximum energy amplification. When the waves are out of phase, the wavepacket magnitude drops.
Wavepacket dynamics from linear solution including (a) temporal signal and (b) its phase scalogram at point of maximum amplitude (
$x = 0.72$).
The phase information from the LES data is presented together with the spanwise-averaged pressure fluctuations in figures 21(a) and 21(b) at 
$x = 0.64$. This location is representative of the vortex pairing region. These figures show that high-pressure disturbances depend on phase alignment of frequencies 
$4 \lesssim St \lesssim 6$. This frequency range corresponds to that where the most amplified modes from the linear analysis are observed. If more (fewer) frequencies are aligned, a stronger (weaker) pressure signal is observed. If the disturbances are in phase, the result is vortex merging and strong spanwise coherence. On the other hand, if they are out of phase, the vortex pairing fails and coherence is reduced, leading to lower values of spanwise-averaged pressure fluctuations. The repeating patterns in the pressure signal show that the feedback mechanism has a period of 
$\Delta t \approx 2.0$, related to the lower tone in the spectrum at 
$St = 0.5$.
Effects of constructive/destructive interference across multiple frequencies presented as the (a) phase-weighted scalogram and its effect in the (b) spanwise-averaged pressure fluctuations.
Destructive and constructive interference of vortical structures during the pairing process leads to an amplitude modulation of the signal (Desquesnes et al. Reference Desquesnes, Terracol and Sagaut2007; Pröbsting et al. Reference Pröbsting, Serpieri and Scarano2014; Ricciardi et al. Reference Ricciardi, Arias-Ramirez and Wolf2020). The modulation mechanism from the 2-D simulations performed in these references differs from the current 3-D case, where turbulent transition occurs. Here, the destructive interference leads to vortex breakdown and smaller pressure fluctuations. On the other hand, constructive interference of multiple frequencies creates spikes in the pressure signal as shown for certain time instants in figure 21 and also in figures 9(a) and 9(c).
4. Conclusions
A study of trailing-edge tonal noise has been conducted through a combination of LES and linear stability theory. The focus of this investigation has been on the secondary tones and acoustic feedback loop which arise in airfoil flows at moderate Reynolds numbers. A compressible LES is performed for the flow over a NACA 0012 airfoil at 
$\alpha = 3^{\circ }$ immersed in a freestream flow with Mach number 
$M_{\infty } = 0.3$ and Reynolds number 
$Re = 5 \times 10^{4}$. Despite intermittent flow transition to turbulent regime, the noise spectrum depicts a main tone with multiple equidistant secondary tones, and all tonal peaks observed are integer multiples of the lowest-frequency peak.
The mean flow displays two recirculation regions, one on each side of the airfoil, where the separation bubble on the suction side is longer than that on the pressure side. The flow dynamics on the airfoil is dominated by events on the suction side, aft of the separation bubble, as indicated by the RMS values of kinetic energy and pressure. The mean friction coefficient shows a double detachment pattern, typical of large fluctuations downstream the recirculation bubble. Mean velocity profiles show that the maximum reversed velocity inside the separation bubble is 13 % of the freestream value.
Linear stability analysis of the mean flow using a bi-global operator shows marginally stable and unstable eigenvalues at frequencies close to those of the tonal peaks observed in the nonlinear simulation, indicating that the overall flow dynamics originates from linear mechanisms. Results from resolvent analysis of the linear operator are displayed as the pseudospectrum, showing that the higher frequencies are more susceptible to be excited by pseudoresonances. The response modes show that the suction side recirculation bubble acts as an amplifier of disturbances. On the other hand, the forcing modes show a high-sensitivity region within the flow, starting downstream of the leading edge and extending until upstream of the bubble.
In the LES, the amplification of fluctuations along the LSB results ultimately in vortex shedding on the suction side. The vortices undergo a pairing process that may lead to either merging or bursting. This process moves the dominant tonal peak in the LES to a lower frequency than that observed in the unstable eigenvalues of the linear spectrum. Another result of the vortex pairing is that either coherent vortices or turbulent spots are advected towards the trailing edge. The success of vortex merging aft of the bubble depends on the phase alignment between the unstable frequencies observed in the linear eigenspectrum. When such frequencies are in phase, it is likely that the vortices will successfully merge. On the other hand, destructive interference results in vortex bursting. The different regimes where coherent structures alternate with uncorrelated eddies at the trailing edge act as an amplitude modulation of the signal. The magnitude of the radiated acoustic waves depends on the coherence level of the vortical structures reaching the trailing edge. A time-frequency analysis using CWTs shows that the flow has two different time scales. One is related to the lowest tonal frequency, which drives the acoustic feedback loop and depends not only on the coherent structures but also on the uncorrelated turbulent packets. The other is a faster time scale related to the passage of multiple coherent structures at the trailing edge. It has been found that the intermittent dynamics of this faster time scale affects the magnitude of tonal peaks with the possibility of changing the instantaneous dominant tonal frequency.
The feedback loop mechanism has been investigated by the time integration of the LNS equations using the mean flow perturbation technique. The system response to an impulsive actuation at the leading edge displays a wavepacket advected towards the trailing edge with subsequent radiation of acoustic waves due to a scattering mechanism. The wavepacket is composed of multiple frequencies in the eigenspectrum and its magnitude is maximum when the phases of the unstable frequencies from the eigenspectrum are aligned. Visualisation of the linearised solution shows that entropy disturbances are triggered by the secondary diffraction and propagate downstream along the boundary layer. These disturbances play an important role since they interact with upstream-propagating waves at the region of maximum forcing. Thus, it is likely that the feedback loop closes at the most sensitive region, described by the forcing modes, in order to satisfy the low-frequency period of the wavepacket cycle.
The pressure RMS values are analysed along the airfoil chord and show an influence of the feedback mechanism, where specific patterns repeat in time. Intermittency due to transition to turbulence may shift the phase and magnitude of all frequencies, disrupting the feedback loop. Despite this, the flow reaches a new equilibrium and reestablishes the merging/bursting process, indicating that the multiple stability modes may be the self-sustaining mechanism for multiple tones and flow transition.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2022.129.
Acknowledgements
The authors thank CENAPAD-SP (Project 551), UCLA cluster Hoffman2 and LNCC-Cluster SDumont (Project SimTurb) for providing the computational resources used in this study. We also thank C.-A. Yeh, C. Skene and J.H.M. Ribeiro for insightful discussions on stability and resolvent analyses.
Funding
T.R.R. and W.R.W. acknowledge Fundação de Amparo à Pesquisa do Estado de São Paulo, FAPESP, for supporting the present work under research grants no. 2013/08293-7, 2018/11835-0, 2019/20437-0 and 2021/06448-0, and Conselho Nacional de Desenvolvimento Científico e Tecnológico, CNPq, for supporting this research under grants no. 407842/2018-7 and 304335/2018-5. K.T. is grateful for support from the US Office of Naval Research under grant N00014-19-1-2460.
Declaration of interests
The authors report no conflict of interest.
Open access
