2.1 Asymptotic regimes

2.1.1 Euler limit ( $\tilde{R}e\rightarrow \infty$ )

In the limit of infinite acoustic Reynolds number (i.e.  $\unicode[STIX]{x1D707}\unicode[STIX]{x1D714}\rightarrow 0$ ), $\text{det}(\unicode[STIX]{x1D647})=0$ simplifies to $\unicode[STIX]{x1D706}(\unicode[STIX]{x1D701}^{2}-\unicode[STIX]{x1D706}^{2})=0$ , giving $\unicode[STIX]{x1D706}=0$ and $\unicode[STIX]{x1D706}=\pm \unicode[STIX]{x1D701}$ , or equivalently $\unicode[STIX]{x1D701}=1/\mathscr{M}$ and $\unicode[STIX]{x1D701}=1/(\mathscr{M}\pm 1)$ . The eigenmodes are

(2.8a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70E}=\left\{(\unicode[STIX]{x1D701},\tilde{\unicode[STIX]{x1D741}})\in \mathbb{R}\times \mathbb{R}^{3}\,:\,\unicode[STIX]{x1D701}={\displaystyle \frac{1}{\mathscr{M}}},\tilde{\unicode[STIX]{x1D741}}=\unicode[STIX]{x1D719}\left[\begin{array}{@{}c@{}}1\\ 0\\ 0\end{array}\right]\right\}, & & \displaystyle\end{eqnarray}$$

(2.8b ) $$\begin{eqnarray}\displaystyle a^{\pm }=\left\{(\unicode[STIX]{x1D701},\tilde{\unicode[STIX]{x1D741}})\in \mathbb{R}\times \mathbb{R}^{3}\,:\,\unicode[STIX]{x1D701}={\displaystyle \frac{1}{\mathscr{M}_{\pm }}},\tilde{\unicode[STIX]{x1D741}}=\unicode[STIX]{x1D719}\left[\begin{array}{@{}c@{}}1\\ \pm 1\\ \unicode[STIX]{x03B3}\end{array}\right]\right\}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D719}\in \mathbb{R}_{\neq 0}$ is a non-zero parameter (which in practice is absorbed by the smallness parameter, $\unicode[STIX]{x1D700}$ ) and $\mathscr{M}_{\pm }\equiv \mathscr{M}\pm 1$ . These eigenmodes are the so-called Kovásznay eigenmodes (see Kovásznay Reference Kovásznay1953, for details). They correspond to entropy waves ( $\unicode[STIX]{x1D70E}$ ) characterised by pure density fluctuations propagating at the flow speed (i.e. with phase speed $u_{0}$ ), and two acoustic waves, one propagating downstream ( $a^{+}$ , with phase speed $u_{0}+c_{0}$ ) and one propagating upstream ( $a^{-}$ , with phase speed $u_{0}-c_{0}$ ). These eigenmodes form a basis that is extensively discussed in Alferez & Touber (Reference Alferez and Touber2017) in the context of arbitrary fluids. In this work, $\unicode[STIX]{x1D70E}$ , $a^{+}$ and $a^{-}$ act as useful markers in the $\text{Re}(\unicode[STIX]{x1D701})$ – $\text{Im}(\unicode[STIX]{x1D701})$ plane since they all lie on the real axis (i.e. these modes have no spatial damping).

2.1.2 Landau regime ( $\tilde{R}e\gg 1$ )

In the limit of large but finite acoustic Reynolds number (i.e.  $\unicode[STIX]{x1D707}\unicode[STIX]{x1D714}\ll \unicode[STIX]{x03B3}p_{0}$ ), the wavenumber ( $k$ , or its dimensionless form $\unicode[STIX]{x1D701}$ ) is complex, reflecting the damping effect of transport coefficients on the Euler modes. Analytical expressions for the departure from the Euler modes in the $\text{Re}(\unicode[STIX]{x1D701})$ – $\text{Im}(\unicode[STIX]{x1D701})$ plane are obtained by seeking a solution in the form $\unicode[STIX]{x1D701}=a+\text{i}b/\tilde{R}e$ where $a$ assumes the $\unicode[STIX]{x1D701}$ value obtained from the Euler modes. Injecting such form in (2.7), expanding in powers of $1/\tilde{R}e$ , and cancelling terms up to first order provide expressions for $b$ . A similar approach is used to find the eigenvectors (i.e. seeking a solution in the form $\tilde{\unicode[STIX]{x1D741}}=\tilde{\unicode[STIX]{x1D741}}_{0}+(\boldsymbol{a}+\text{i}\boldsymbol{b})/\tilde{R}e$ up to first order in $1/\tilde{R}e$ , where $\tilde{\unicode[STIX]{x1D741}}_{0}$ is the corresponding eigenvector in the Euler limit). The resulting eigenmodes are

(2.9a ) $$\begin{eqnarray}\displaystyle \mathscr{S}=\left\{(\unicode[STIX]{x1D701},\tilde{\unicode[STIX]{x1D741}})\in \mathbb{C}\times \mathbb{C}^{3}\,:\,\unicode[STIX]{x1D701}={\displaystyle \frac{1}{\mathscr{M}}}-\frac{\text{i}}{\tilde{R}e\mathit{Pr}\mathscr{M}^{3}},\tilde{\unicode[STIX]{x1D741}}=\unicode[STIX]{x1D719}\left[\begin{array}{@{}c@{}}1\\ \text{i}/[\tilde{R}e\mathit{Pr}\mathscr{M}]\\ 0\end{array}\right]\right\}, & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(2.9b ) $$\begin{eqnarray}\displaystyle \mathscr{L}^{\pm }=\left\{(\unicode[STIX]{x1D701},\tilde{\unicode[STIX]{x1D741}})\in \mathbb{C}\times \mathbb{C}^{3}\,:\,\unicode[STIX]{x1D701}={\displaystyle \frac{1}{\mathscr{M}_{\pm }}}-\frac{\text{i}\unicode[STIX]{x1D6FF}^{+}}{2\tilde{R}e\mathscr{M}_{\pm }^{3}},\tilde{\unicode[STIX]{x1D741}}=\unicode[STIX]{x1D719}\left[\begin{array}{@{}c@{}}1\mp {\displaystyle \frac{\text{i}(\unicode[STIX]{x03B3}-1)}{\tilde{R}e\mathit{Pr}\mathscr{M}_{\pm }}}\\ \pm 1+{\displaystyle \frac{\text{i}\unicode[STIX]{x1D6FF}^{-}}{2\tilde{R}e\mathscr{M}_{\pm }}}\\[12.0pt] \unicode[STIX]{x03B3}\end{array}\right]\right\}, & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D719}\in \mathbb{C}_{\neq 0}$ and $\unicode[STIX]{x1D6FF}^{\pm }\equiv 4/3+\unicode[STIX]{x1D712}\pm (\unicode[STIX]{x03B3}-1)/\mathit{Pr}$ (note that $\tilde{\unicode[STIX]{x1D741}}$ can be arbitrarily stretched and rotated on the complex plane by $\unicode[STIX]{x1D719}$ , it is set here to align with $\tilde{p}$ in (2.8b ) for $\unicode[STIX]{x1D719}=1$ ). The damping coefficient $\text{Im}(\unicode[STIX]{x1D701})=-\unicode[STIX]{x1D6FF}^{+}/(2\tilde{R}e\mathscr{M}_{\pm }^{3})$ for the acoustic modes $\mathscr{L}^{\pm }$ (with $\unicode[STIX]{x1D6FF}^{+}$ the non-dimensional Kirchoff diffusivity), corresponds to the well-known acoustic damping found by Landau & Lifshitz (Reference Landau and Lifshitz1987) (equation (79.6) on page 301 of the 2nd edition), which is corrected here to account for a non-zero base-flow velocity (i.e. when $u_{0}\neq 0$ ). The above result also offers a more complete picture since the damping is clearly identified as a first-order correction (in $1/\tilde{R}e$ ) to Kovásznay eigenmode basis (including that of the entropy mode $\mathscr{S}$ ), with explicit expressions for the coordinates ( $\unicode[STIX]{x1D70C}$ , $u$ , $p$ ) of each eigenmode in the primitive-variable basis. This will be useful when discussing the effect of $\unicode[STIX]{x1D707}_{b}$ on turbulence.

2.1.3 Stokes regime ( $\tilde{R}e\ll 1$ )

In the limit of small acoustic Reynolds number (i.e.  $\unicode[STIX]{x1D707}\unicode[STIX]{x1D714}\gg \unicode[STIX]{x03B3}p_{0}$ ), the wavenumber ( $k$ , or its dimensionless form $\unicode[STIX]{x1D701}$ ) is again complex but is anticipated to align with the (stable) diagonal of the $\text{Re}(\unicode[STIX]{x1D701})$ – $\text{Im}(\unicode[STIX]{x1D701})$ plane. This reflects the tendency towards a diffusion-dominated flow, and in particular the celebrated ‘second problem of Stokes’ (i.e. oscillating plate), which has solutions in the form $\sin (\unicode[STIX]{x1D714}t-Kx)\exp (-Kx)$ with $Kc_{0}/\unicode[STIX]{x1D714}=(\tilde{R}e/2)^{1/2}$ . Thus, one expects that in the limit of small acoustic Reynolds number, $\unicode[STIX]{x1D701}=\unicode[STIX]{x1D702}\unicode[STIX]{x1D6E4}\exp (-\text{i}\unicode[STIX]{x03C0}/4)$ where $\unicode[STIX]{x1D702}\in \mathbb{R}_{\neq 0}$ and $\unicode[STIX]{x1D6E4}\equiv \tilde{R}e^{1/2}$ . Injecting such form in (2.7) and expanding in powers of $\unicode[STIX]{x1D6E4}$ gives

(2.10) $$\begin{eqnarray}\displaystyle \text{det}(\unicode[STIX]{x1D647})\underset{\unicode[STIX]{x1D6E4}\rightarrow 0}{=}\left\{A+{\displaystyle \frac{(B-A)\unicode[STIX]{x1D702}\mathscr{M}}{\sqrt{2}}}\unicode[STIX]{x1D6E4}+O(\unicode[STIX]{x1D6E4}^{2})\right\}-\text{i}\left\{{\displaystyle \frac{(B-A)\unicode[STIX]{x1D702}\mathscr{M}}{\sqrt{2}}}\unicode[STIX]{x1D6E4}+O(\unicode[STIX]{x1D6E4}^{2})\right\}, & & \displaystyle\end{eqnarray}$$

where $A\equiv \unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D702}^{4}-(\unicode[STIX]{x1D6FC}+\unicode[STIX]{x1D6FD})\unicode[STIX]{x1D702}^{2}+1$ , $B\equiv (\unicode[STIX]{x1D6FC}+\unicode[STIX]{x1D6FD})\unicode[STIX]{x1D702}^{2}-2$ , with $\unicode[STIX]{x1D6FC}\equiv 4/3+\unicode[STIX]{x1D712}$ , $\unicode[STIX]{x1D6FD}\equiv \unicode[STIX]{x03B3}/\mathit{Pr}$ (note that $A$ and $B$ are real). A necessary (and sufficient) condition for non-trivial solutions to (2.5) to exist (in the required form and up to first order in $\unicode[STIX]{x1D6E4}$ ) is $A=0$ and $B\mathscr{M}=0$ . Setting $A=0$ gives $\unicode[STIX]{x1D702}^{2}=1/\unicode[STIX]{x1D6FC}$ or $\unicode[STIX]{x1D702}^{2}=1/\unicode[STIX]{x1D6FD}$ . If $\mathscr{M}\neq 0$ the second statement ( $B\mathscr{M}=0$ ) requires that $\unicode[STIX]{x1D702}^{2}=2/(\unicode[STIX]{x1D6FC}+\unicode[STIX]{x1D6FD})$ , which is compatible with the condition $A=0$ if $\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D6FD}$ . If $\mathscr{M}=0$ ( $u_{0}=0$ ), the condition $A=0$ suffices. In practice, setting either $\mathscr{M}=0$ or $\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D6FD}$ may be too restrictive. Instead, let $u_{\unicode[STIX]{x1D719}}\equiv \unicode[STIX]{x1D714}/\text{Re}(k)$ denote the phase speed and note that $\unicode[STIX]{x1D702}\mathscr{M}\unicode[STIX]{x1D6E4}/\sqrt{2}=u_{0}/u_{\unicode[STIX]{x1D719}}$ . Thus, the $O(\unicode[STIX]{x1D6E4})$ terms in (2.10) may be neglected when considering large phase speeds (compared to that of the base flow, i.e.  $u_{\unicode[STIX]{x1D719}}\gg u_{0}$ , independently of how large $u_{0}$ is), and setting $A=0$ becomes a sufficient condition for non-trivial solutions to exist. The eigenmodes corresponding to roots $\unicode[STIX]{x1D702}=\pm 1/\sqrt{\unicode[STIX]{x1D6FC}}$ and $\unicode[STIX]{x1D702}=\pm 1/\sqrt{\unicode[STIX]{x1D6FD}}$ (in the limit of $u_{0}/u_{\unicode[STIX]{x1D719}}\rightarrow 0$ ) are

(2.11a ) $$\begin{eqnarray}\displaystyle \mathscr{U}^{\pm }=\left\{(\unicode[STIX]{x1D701},\tilde{\unicode[STIX]{x1D741}})\in \mathbb{C}\times \mathbb{C}^{3}\,:\,\unicode[STIX]{x1D701}=\pm \sqrt{\frac{\tilde{R}e}{2\unicode[STIX]{x1D6FC}}}\left(1-\text{i}\right),\tilde{\unicode[STIX]{x1D741}}=\unicode[STIX]{x1D719}\left[\begin{array}{@{}c@{}}0\\ 1\\ 0\end{array}\right]\right\}, & & \displaystyle\end{eqnarray}$$

(2.11b ) $$\begin{eqnarray}\displaystyle \mathscr{P}^{\pm }=\left\{(\unicode[STIX]{x1D701},\tilde{\unicode[STIX]{x1D741}})\in \mathbb{C}\times \mathbb{C}^{3}\,:\,\unicode[STIX]{x1D701}=\pm \sqrt{\frac{\tilde{R}e}{2\unicode[STIX]{x1D6FD}}}\left(1-\text{i}\right),\tilde{\unicode[STIX]{x1D741}}=\unicode[STIX]{x1D719}\left[\begin{array}{@{}c@{}}0\\ 0\\ 1\end{array}\right]\right\}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D719}\in \mathbb{C}_{\neq 0}$ . They correspond to the diffusion of velocity ( $\mathscr{U}^{\pm }$ ) and pressure ( $\mathscr{P}^{\pm }$ ) fluctuations in ( $+$ ) and against ( $-$ ) the flow direction (‘ploughing’ space at arbitrarily fast phase speeds). For a given shear viscosity $\unicode[STIX]{x1D707}$ , these eigenmodes become increasingly relevant to small-wavelength waves (high-frequency ‘sound’). Note the scaling in $\tilde{R}e^{1/2}$ , implying that higher frequencies are less damped than lower frequencies. This is to contrast with the Landau regime, where the opposite is true.

Finally, if $\mathscr{M}\neq 0$ , injecting $\unicode[STIX]{x1D701}=(1+\text{i}b\tilde{R}e)/\mathscr{M}$ in (2.7) gives

(2.12) $$\begin{eqnarray}\displaystyle \displaystyle \mathscr{M}^{4}\text{det}(\unicode[STIX]{x1D647}) & = & \displaystyle \{4b(X^{2}-1)Z-2YX^{2}\mathscr{M}^{2}\}\nonumber\\ \displaystyle & & \displaystyle +\,\displaystyle \text{i}\left\{X[X^{2}(\mathscr{M}^{2}-Y)+Y]\mathscr{M}^{2}+[(1-X^{2})^{2}-4X^{2}]\frac{Z}{\tilde{R}e}\right\},\end{eqnarray}$$

where $X\equiv b\tilde{R}e$ , $Y\equiv 1+b(\unicode[STIX]{x1D6FC}+\unicode[STIX]{x1D6FD})$ and $Z\equiv \unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}b+1/\mathit{Pr}$ . Setting $Z=0$ (i.e.  $b=-1/(\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}\mathit{Pr})=-1/(\unicode[STIX]{x03B3}\unicode[STIX]{x1D6FC})$ ) then gives

(2.13) $$\begin{eqnarray}\displaystyle \text{det}(\unicode[STIX]{x1D647})\underset{\tilde{R}e\rightarrow 0}{=}\{O(\tilde{R}e^{2})\}+\text{i}\{O(\tilde{R}e)\}. & & \displaystyle\end{eqnarray}$$

Therefore, in the limit of vanishing acoustic Reynolds number ( $\tilde{R}e\rightarrow 0$ ) and in the presence of a base flow ( $u_{0}\neq 0$ ), the following eigenmode is also considered

(2.14) $$\begin{eqnarray}\displaystyle \displaystyle \mathscr{N}=\left\{(\unicode[STIX]{x1D701},\tilde{\unicode[STIX]{x1D741}})\in \mathbb{C}\times \mathbb{R}^{3}\,:\,\unicode[STIX]{x1D701}=\frac{1}{\mathscr{M}}-\frac{\text{i}\tilde{R}e}{\unicode[STIX]{x03B3}\unicode[STIX]{x1D6FC}\mathscr{M}},\tilde{\unicode[STIX]{x1D741}}=\unicode[STIX]{x1D719}\left[\begin{array}{@{}c@{}}1\\ 0\\ 1\end{array}\right]\right\}. & & \displaystyle\end{eqnarray}$$

Whilst reminiscent of the Euler entropy mode ( $\unicode[STIX]{x1D70E}$ ), mode $\mathscr{N}$ has equal coordinates in $\unicode[STIX]{x1D70C}$ and $p$ , i.e.  $|p-p_{0}|/p_{0}\sim |\unicode[STIX]{x1D70C}-\unicode[STIX]{x1D70C}_{0}|/\unicode[STIX]{x1D70C}_{0}\sim \unicode[STIX]{x1D700}\unicode[STIX]{x1D719}$ . Noting that $p(\unicode[STIX]{x1D70C}_{0}+\unicode[STIX]{x1D700}\unicode[STIX]{x1D70C}_{0},T_{0})-p(\unicode[STIX]{x1D70C}_{0},T_{0})\sim \unicode[STIX]{x1D700}\unicode[STIX]{x1D70C}_{0}(\unicode[STIX]{x2202}p/\unicode[STIX]{x2202}\unicode[STIX]{x1D70C})_{T=T_{0}}=\unicode[STIX]{x1D700}\unicode[STIX]{x1D70C}_{0}a_{0}^{2}=\unicode[STIX]{x1D700}p_{0}$ (where $a_{0}$ is Newton’s isothermal speed of sound), mode $\mathscr{N}$ corresponds to the transport by the base flow of isothermally compressed (or expanded) fluid particles (relative to the base state). It is referred here as Newton mode in reference to the isothermal compression of fluid particles. The Newton mode ( $\mathscr{N}$ ) will be shown to smoothly reconnect to the Euler entropy mode ( $\unicode[STIX]{x1D70E}$ ) when the bulk viscosity is large, and to play a primary role in promoting the solenoidal kinetic energy of small-scale turbulence.

2.1.4 Graphical summary of the asymptotic regimes

Figure 1. Graphical summary of the asymptotic limits of $\text{det}(\unicode[STIX]{x1D647})=0$ for $\mathscr{M}<1$ (a) and $\mathscr{M}>1$ (b) representing the damping factor $\text{Im}(\unicode[STIX]{x1D701})$ as a function of the phase Mach number $\mathscr{M}_{\unicode[STIX]{x1D719}}=1/\text{Re}(\unicode[STIX]{x1D701})$ . Arrows point in the direction of increasing $\tilde{R}e$ values. On the real axis ( $\mathscr{M}_{\unicode[STIX]{x1D719}}$ ) lie the Euler solutions ( $a^{\pm }$ , $\unicode[STIX]{x1D70E}$ ), which are reached by the Landau branches ( $\mathscr{L}^{\pm }$ , $\mathscr{S}$ ) when $\tilde{R}e\rightarrow \infty$ . In contrast, the Newton branch ( $\mathscr{N}$ ) and Stokes diffusion branches ( $\mathscr{U}^{\pm }$ , $\mathscr{P}^{\pm }$ ) depart from the real axis as $\tilde{R}e$ is increased (vertically down for $\mathscr{N}$ and according to $-1/\mathscr{M}_{\unicode[STIX]{x1D719}}$ for $\mathscr{U}$ and $\mathscr{P}$ ). The eight branches shown here act as guidelines to actual solutions to $\text{det}(\unicode[STIX]{x1D647})=0$ , which continuously transition from one Stokes branch to one Landau branch, as shown in figure 2.

Figure 1 provides a qualitative summary of the asymptotic limits discussed in the previous paragraphs. The vertical axis represents the damping factor $\text{Im}(\unicode[STIX]{x1D701})$ . The horizontal axis represents the phase Mach number, $\mathscr{M}_{\unicode[STIX]{x1D719}}=1/\text{Re}(\unicode[STIX]{x1D701})$ . In the limit of infinite acoustic Reynolds number, the dispersion relation converges to the Euler modes ( $a^{\pm }$ , $\unicode[STIX]{x1D70E}$ ). These limits are reached via the Landau modes ( $\mathscr{S}$ , $\mathscr{L}^{\pm }$ ) presented in § 2.1.2. In the limit of vanishing acoustic Reynolds number, the dispersion relation converges to the Newton ( $\mathscr{N}$ ) and Stokes diffusion ( $\mathscr{U}^{\pm }$ , $\mathscr{P}^{\pm }$ ) modes described in § 2.1.3. Both limits ( $\tilde{R}e\rightarrow \infty$ and $\tilde{R}e\rightarrow 0$ ) assume undamped solutions and thus lie on the real axis. In the Stokes diffusion case, such limit takes the form of infinite phase speeds (i.e.  $\mathscr{M}_{\unicode[STIX]{x1D719}}\rightarrow \pm \infty$ when $\tilde{R}e\rightarrow 0$ along the $\mathscr{U}^{\pm }$ and $\mathscr{P}^{\pm }$ branches).

On the $\mathscr{M}_{\unicode[STIX]{x1D719}}$ – $\text{Im}(\unicode[STIX]{x1D701})$ plane, both the Landau ( $\mathscr{S}$ , $\mathscr{L}^{\pm }$ ) and Newton ( $\mathscr{N}$ ) branches connect to the $\mathscr{M}_{\unicode[STIX]{x1D719}}$ axis following a vertical tangent. However, note the directional difference for the $\mathscr{N}$ branch, which departs from $\unicode[STIX]{x1D70E}$ (contrary to the $\mathscr{S}$ mode). The Stokes diffusion modes ( $\mathscr{U}^{\pm }$ , $\mathscr{P}^{\pm }$ ) follow a path with a $-1/\mathscr{M}_{\unicode[STIX]{x1D719}}$ scaling. For a given fluid and base flow, the acoustic Reynolds number ( $\tilde{R}e$ ) plays the role of a coordinate along each branch (corresponding to changing the forcing frequency). With this in mind, note from (2.11) that despite being superimposed, the same point on the $\mathscr{U}$ and $\mathscr{P}$ branches correspond to different $\tilde{R}e$ coordinates (due to the different stretching factors, $1/\sqrt{\unicode[STIX]{x1D6FC}}$ versus $1/\sqrt{\unicode[STIX]{x1D6FD}}$ ).

Owing to the choice of plane-wave definition (see (2.4)), stable modes have a negative damping factor if propagating in the $x>0$ direction, and a positive damping factor if propagating in the $x<0$ direction. All asymptotic regimes give linearly stable modes. For supersonic base flows ( $\mathscr{M}>1$ ), the $\mathscr{L}^{-}$ branch moves from the $(\mathscr{M}_{\unicode[STIX]{x1D719}}<0,\text{Im}(\unicode[STIX]{x1D701})>0)$ quadrant to the $(\mathscr{M}_{\unicode[STIX]{x1D719}}>0,\text{Im}(\unicode[STIX]{x1D701})<0)$ quadrant. This is due to the $\mathscr{M}_{-}$ term in (2.9b ). Unlike for Euler flows, viscous and thermally conductive supersonic flows still sustain two modes ( $\mathscr{U}^{-}$ and $\mathscr{P}^{-}$ ) travelling upstream of the source point, which are less and less damped as the forcing frequency is increased.

2.2 Continuous ‘acoustic’ and ‘entropy’ paths

2.2.1 $\unicode[STIX]{x1D647}$ -matrix spectrum

The asymptotic branches discussed in the previous paragraphs are now combined with numerical solutions of $\text{det}(\unicode[STIX]{x1D647})=0$ to extract and interpret particular paths in the phase-speed ( $\mathscr{M}_{\unicode[STIX]{x1D719}}$ ) and damping-factor ( $\text{Im}(\unicode[STIX]{x1D701})$ ) plane. All numerical solutions presented here consider carbon dioxide in atmospheric conditions (i.e.  $\unicode[STIX]{x03B3}=1.29$ , $\mathit{Pr}=0.7$ ) flowing at $\mathscr{M}=0.5$ . However, the asymptotic solutions given above can readily be used to anticipate the effect of changing the specific heat ratio, the Prandtl number and the base-flow speed on the results presented here.

Figure 2. Damping factor, $\text{Im}(\unicode[STIX]{x1D701})$ , and phase Mach number, $\mathscr{M}_{\unicode[STIX]{x1D719}}$ , of the eigenmodes obtained from numerical solutions of $\text{det}(\unicode[STIX]{x1D647})=0$ for a range of acoustic Reynolds numbers ( $\tilde{R}e\in [10^{-3},10^{5}]$ ) and bulk- to shear-viscosity ratios ( $\unicode[STIX]{x1D712}\in \{0,10,10^{2},10^{3}\}$ , indicated inside each panel). The base flow is set to $\mathscr{M}=0.5$ , together with $\unicode[STIX]{x03B3}=1.29$ and $\mathit{Pr}=0.7$ . The acoustic ( $\mathscr{A}^{\pm }$ ) and entropy ( $\mathscr{E}$ ) paths are identified in red and black, respectively (branches not taking part in these paths are shown in orange). The Euler modes ( $\unicode[STIX]{x1D70E}$ , $a^{\pm }$ ) are indicated along the real axis for reference. A marker (white circle) is shown on the acoustic and entropy paths to indicate the $\tilde{R}e=10$ coordinates, where the line is made thicker on the side corresponding to $\tilde{R}e>10$ . Arrows on each path also give the direction of increasing $\tilde{R}e$ values. A mixture of log–log, linear and semi-log scales is used. The log–log and linear quadrants are marked with thin grey lines.

Figure 2 gives the dispersion relations (i.e. roots of (2.7), obtained numerically) for four different $\unicode[STIX]{x1D707}_{b}/\unicode[STIX]{x1D707}$ ratios ( $\unicode[STIX]{x1D712}\in \{0,10,10^{2},10^{3}\}$ ) over a range of acoustic Reynolds numbers ( $\tilde{R}e\in [10^{-3},10^{5}]$ ). The damping factor, $\text{Im}(\unicode[STIX]{x1D701})$ , is plotted against the phase Mach number, $\mathscr{M}_{\unicode[STIX]{x1D719}}=1/\text{Re}(\unicode[STIX]{x1D701})$ . Since $\mathscr{M}$ , $\unicode[STIX]{x1D707}$ and $\unicode[STIX]{x1D705}$ are all chosen to be non-zero, (2.7) is a fifth-order polynomial, giving five (distinct) roots. Each root describes one branch in figure 2 when varying $\tilde{R}e$ (keeping everything else constant), giving a total of five branches. In agreement with the asymptotic analyses, all branches from the numerical solutions are linearly stable.

The branches are made thicker when $\tilde{R}e>10$ to identify the parts of the branches corresponding to higher acoustic Reynolds numbers. The low- and high- $\tilde{R}e$ ends of each branch are constrained by two of the asymptotic regimes discussed in figure 1. Solutions to $\text{det}(\unicode[STIX]{x1D647})=0$ thus form continuous paths connecting the Stokes and Landau regimes. In reference to the Euler mode they reach (when $\tilde{R}e\rightarrow \infty$ ), three particular paths in the $\mathscr{M}_{\unicode[STIX]{x1D719}}$ – $\text{Im}(\unicode[STIX]{x1D701})$ plane are considered.

1. (i) Two continuous acoustic paths ( $\mathscr{A}^{\pm }$ ): the Landau acoustic mode ( $\mathscr{L}$ ) connects to one of the Stokes diffusion mode ( $\mathscr{U}$ or $\mathscr{P}$ ) and forms an acoustic path ( $\mathscr{A}$ ), shown in red in figure 2. The Stokes regime in the (signed) log–log portions of the plot is characterised by a unitary slope. In the cases shown here, $\mathscr{P}$ is found to connect to $\mathscr{L}$ for $\unicode[STIX]{x1D712}=0$ , and $\mathscr{U}$ to $\mathscr{L}$ for $\unicode[STIX]{x1D712}\in \{10,10^{2},10^{3}\}$ . Remarkably, the shape of the acoustic path is not sensitive to the choice of $\unicode[STIX]{x1D712}$ (for given $\unicode[STIX]{x03B3}$ , $\mathit{Pr}$ and $\mathscr{M}$ values). Rather, $\unicode[STIX]{x1D712}$ plays the role of a stretching factor on the $\tilde{R}e$ coordinate along the acoustic path (this is visible from the $\unicode[STIX]{x1D6FF}^{+}$ and $1/\sqrt{\unicode[STIX]{x1D6FC}}$ factors in (2.9b ) and (2.11a )). This is illustrated by the displacement of the white circle (placed at $\tilde{R}e=10$ ) along the $\mathscr{A}^{\pm }$ paths from one $\unicode[STIX]{x1D712}$ case to another in figure 2 (the marker recedes into the Stokes regime as $\unicode[STIX]{x1D712}$ is increased).

2. (ii) One continuous entropy path ( $\mathscr{E}$ ): the Landau entropy branch ( $\mathscr{S}$ ) either connects to one of the Stokes diffusion branches ( $\mathscr{U}^{+}$ for $\unicode[STIX]{x1D712}=0$ , $\mathscr{P}^{+}$ for $\unicode[STIX]{x1D712}=10$ ) or the Newton branch ( $\mathscr{N}$ for $\unicode[STIX]{x1D712}\in \{10^{2},10^{3}\}$ ). The resulting continuous entropy paths are shown in black in figure 2. The entropy path for the high- $\unicode[STIX]{x1D712}$ cases shown here therefore describes a closed loop in $\mathscr{M}_{\unicode[STIX]{x1D719}}$ – $\text{Im}(\unicode[STIX]{x1D701})$ space (figure 2 c,d).

The (two) remaining branches, shown in orange in figure 2, diverge as $\tilde{R}e\rightarrow \infty$ . These modes vanish (they are infinitely damped) in the Euler limit. If the base flow is supersonic ( $\mathscr{M}>1$ , not shown here), then the remaining (orange) branch on the $\mathscr{M}_{\unicode[STIX]{x1D719}}>0$ side of figure 2 (the one not connected to either $\unicode[STIX]{x1D70E}$ or $a^{+}$ ) is the branch connecting to the upstream-propagating acoustic mode ( $a^{-}$ , which will now be on the $\mathscr{M}_{\unicode[STIX]{x1D719}}>0$ side), whilst the two upstream-going Stokes diffusion branches ( $\mathscr{U}^{-}$ and $\mathscr{P}^{-}$ ) are the ones diverging in the Euler limit and not forming a part of the acoustic or entropy paths.

Figure 3. Magnitude of the real (thick grey) and imaginary (blue) parts of the wavenumber $\unicode[STIX]{x1D701}$ along the $\mathscr{A}^{+}$ (a,b) and $\mathscr{E}$ (c,d) paths for $\unicode[STIX]{x1D712}\in \{0,10\}$ (a,c) and $\unicode[STIX]{x1D712}\in \{10^{2},10^{3}\}$ (b,d), obtained from solving $\text{det}(\unicode[STIX]{x1D647})=0$ . The solid and dashed black lines give the asymptotic solutions from §§ 2.1.2 and 2.1.3 (the scalings and the branch types are reproduced here for convenience).

The acoustic and entropy paths are parametrised by the acoustic Reynolds number. Figure 3 gives the evolution of the magnitude of the real and imaginary parts of the wavenumber $\unicode[STIX]{x1D701}$ along the $\mathscr{A}^{+}$ and $\mathscr{E}$ paths (for all four values of $\unicode[STIX]{x1D712}$ ) obtained from numerical solutions to $\text{det}(\unicode[STIX]{x1D647})=0$ . Superimposed are the different asymptotic solutions described in §§ 2.1.2 and 2.1.3. The asymptotic solutions are in remarkable agreement with the numerical results over the Stokes and Landau regimes. The transition between the two regimes is limited to about one decade in $\tilde{R}e$ . In all cases, the transition corresponds to a maximum in the damping factor. For the acoustic path, $\unicode[STIX]{x1D712}$ operates a shift in $\tilde{R}e$ of the curve without modifying its overall shape. This confirms the above observation that the acoustic path is unchanged (for given $\unicode[STIX]{x03B3}$ , $\mathit{Pr}$ and $\mathscr{M}$ values) and that $\unicode[STIX]{x1D712}$ provides an alternative way of exploring this path at fixed $\tilde{R}e$ . In contrast, the entropy path cannot be explored via $\unicode[STIX]{x1D712}$ alone: first, because the Landau regime ( $\mathscr{S}$ ) is independent of $\unicode[STIX]{x1D712}$ ; second, because the path itself changes with $\unicode[STIX]{x1D712}$ in the Stokes regime either because a different branch is selected or because the intercept with the Landau regime shifts (see figure 3 d) owing to the fixed position of the Landau branch $\mathscr{S}$ . Note that this last case implies that the maximum damping on the primitive variables associated with the entropy path will reduce as $\unicode[STIX]{x1D712}$ is increased, a rather counterintuitive effect of $\unicode[STIX]{x1D707}_{b}$ (i.e. more bulk viscosity leads to less damping).

2.2.2 Eigenvectors

Figure 4. Eigenvector trajectories along the $\mathscr{A}^{+}$ (red) and $\mathscr{E}$ (black) paths for $\unicode[STIX]{x1D712}=0$ (a), $\unicode[STIX]{x1D712}=10$ (b) and $\unicode[STIX]{x1D712}\in \{10^{2},10^{3}\}$ (c) on the surface of one eighth of the unitary sphere in the Euclidean space formed by the dimensionless primitive variable. Since $\tilde{\unicode[STIX]{x1D741}}\in \mathbb{C}^{3}$ , the moduli of the components of $\tilde{\unicode[STIX]{x1D741}}$ (divided by the magnitude of $\tilde{\unicode[STIX]{x1D741}}$ ) are plotted instead (i.e.  $(|\tilde{\unicode[STIX]{x1D70C}}|,|\tilde{u} |,|\tilde{p}|)/\Vert \tilde{\unicode[STIX]{x1D741}}\Vert$ ). Each path is split into the Stokes (solid line) and Landau (dotted line) regimes with the boundary defined as the $\tilde{R}e$ value for which $\text{Im}(\unicode[STIX]{x1D701})$ in figure 3 reaches its maximum. All data are based on numerical solutions to (2.5). In the limit of $\tilde{R}e\rightarrow \infty$ the acoustic and entropy trajectories terminate on the $a^{+}$ and $\unicode[STIX]{x1D70E}$ points, respectively. Once $\unicode[STIX]{x1D712}$ is sufficiently high for $\mathscr{U}^{+}$ to connect to $\mathscr{L}^{+}$ and $\mathscr{N}$ to connect to $\mathscr{S}$ , the acoustic and entropy paths on the sphere vary very little when varying $\unicode[STIX]{x1D712}$ . In particular, the entropy path then lies near the equator (going from $\mathscr{N}$ to $\unicode[STIX]{x1D70E}$ ), and the acoustic path lies near the meridian connecting $\mathscr{U}^{+}$ to $a^{+}$ .

The relative amplitudes of the primitive variables ( $\unicode[STIX]{x1D70C}$ , $u$ , $p$ ) along the $\mathscr{A}^{+}$ and $\mathscr{E}$ paths on the surface of one eighth of the unit sphere are explored in figure 4 (see caption for details), enabling the visualisation of the smooth transition from one of the Stokes limit ((2.9) and (2.11)) to the Euler limit (2.8). For the acoustic path, starting from either a pure fluctuation in $p$ ( $\mathscr{P}^{+}$ ) or $u$ ( $\mathscr{U}^{+}$ ), most of the displacement on the sphere occurs in the Stokes regime. The effect of the damping factor on the fluctuating fields is therefore expected to be easier to read in the Landau regime (where the eigenvector is not changing much) than in the Stokes regime (where the fluctuating field will be influenced by both changes in the damping factor and variations in the eigenvector components). Apart from when it is connected to $\mathscr{U}^{+}$ ( $\unicode[STIX]{x1D712}=0$ case), the entropy path is seen to stay close to the equator of the sphere where $u=0$ . At high $\unicode[STIX]{x1D712}$ values the path becomes insensitive to the actual value of $\unicode[STIX]{x1D712}$ and sees its largest excursion in the Landau regime. Under such conditions, the fluctuating fields associated with the entropy path can therefore be directly anticipated from the properties of the damping factor. Specifically, pronounced density and pressure fluctuations carried by the base flow are expected to arise when $\tilde{R}e\rightarrow 0$ (where $\text{Im}(\unicode[STIX]{x1D701})\rightarrow 0$ ).

2.3 The role of bulk viscosity in shaping the fluctuating field

Figure 5. Spatial evolution of $u^{\prime }\equiv u_{1}/c_{0}$ (where $\tilde{x}\equiv x/\unicode[STIX]{x1D706}_{0}$ with $\unicode[STIX]{x1D706}_{0}\equiv 2\unicode[STIX]{x03C0}c_{0}/\unicode[STIX]{x1D714}$ is the distance in units of the Euler acoustic wavelength) based on (2.15) over a range of acoustic Reynolds numbers ( $\tilde{R}e$ ) at time $tc_{0}/\unicode[STIX]{x1D706}_{0}=0.1\,(\text{mod}\,2\unicode[STIX]{x03C0})$ (see Movie 1 for a full period), for $\unicode[STIX]{x1D712}\in \{0,10,10^{2},10^{3}\}$ , $\unicode[STIX]{x03B3}=1.29$ , $\mathit{Pr}=0.7$ and $\mathscr{M}=u_{0}/c_{0}=0.5$ . The lighting of the $u^{\prime }$ surface is such that positive density–pressure correlations are highlighted and coloured according to the dilatational rate (see colour scale for $\unicode[STIX]{x1D70C}^{\prime }\equiv \unicode[STIX]{x1D70C}_{1}/\unicode[STIX]{x1D70C}_{0}$ , $p^{\prime }\equiv p_{1}/p_{0}$ , $\unicode[STIX]{x1D703}^{\prime }\equiv \unicode[STIX]{x03B8}\unicode[STIX]{x1D706}_{0}/c_{0}$ ). A triangular marker indicates the $\tilde{R}e$ location where one asymptotic-mode amplitude ( $\mathscr{P}^{\pm }$ , $\mathscr{U}^{\pm }$ , $\mathscr{L}^{\pm }$ , $\mathscr{N}$ , $\mathscr{S}$ ) has dropped by 10 % from its $\tilde{x}=0$ at the domain edge. Downward-pointing triangular markers correspond to $\unicode[STIX]{x1D712}$ -independent modes ( $\mathscr{P}^{\pm }$ , $\mathscr{S}$ ). Note the quiet regions (valleys) corresponding to the transition from the Landau to the Stokes regimes (when varying $\unicode[STIX]{x1D707}\unicode[STIX]{x1D714}$ ). The valleys shift to higher $\tilde{R}e$ values as $\unicode[STIX]{x1D712}$ is increased but is increasingly disturbed by $\unicode[STIX]{x1D70C}^{\prime }$ from the fixed $\mathscr{S}$ branch. On the Stokes side of the valley, long (in $\unicode[STIX]{x1D706}_{0}$ units) coherent structures in $u^{\prime }$ – $\unicode[STIX]{x1D703}^{\prime }$ ( $\mathscr{U}$ ) and $p^{\prime }$ ( $\mathscr{P}$ ) structures combine downstream with the isothermal pressure–density fluctuations ( $\mathscr{N}$ mode).

The variation of both the dispersion relation and the eigenvectors along all five continuous paths ( $\mathscr{A}^{\pm }$ , $\mathscr{E}$ and the remaining two Stokes branches) are now recombined to investigate their net effect on the fluctuating fields (again in the context of carbon dioxide with $\unicode[STIX]{x03B3}=1.29$ , $\mathit{Pr}=0.7$ flowing at $\mathscr{M}=0.5$ ), defined as

(2.15) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D741}_{1}=\left[\begin{array}{@{}ccc@{}}\unicode[STIX]{x1D70C}_{0} & 0 & 0\\ 0 & c_{0} & 0\\ 0 & 0 & p_{0}\end{array}\right]\text{Im}\left\{\mathop{\sum }_{n=1}^{5}\tilde{\unicode[STIX]{x1D741}}_{n}\,\text{exp}\left[\text{i}\left(\unicode[STIX]{x1D714}t-k_{n}x+\unicode[STIX]{x1D711}_{n}\right)\right]\right\}, & & \displaystyle\end{eqnarray}$$

where the $(k_{n},\tilde{\unicode[STIX]{x1D741}}_{n})$ pairs are the five solutions to (2.5). The eigenvectors, $\tilde{\unicode[STIX]{x1D741}}_{n}$ , are made unitary (based on the Euclidean norm) and the phase lags $\unicode[STIX]{x1D711}_{n}$ are chosen so as to align all $\tilde{u}$ components in the complex plane (to interpret as a forcing at $x=0$ synchronised on the velocity fluctuation). Figure 5 gives the spatial evolution of the velocity fluctuations, combined with those of density, pressure and dilatational rate, for a range of acoustic Reynolds numbers, where space is made dimensionless using the acoustic wavelength, $\unicode[STIX]{x1D706}_{0}\equiv 2\unicode[STIX]{x03C0}c_{0}/\unicode[STIX]{x1D714}$ .

A prominent feature of the fluctuating field is the existence of a range of acoustic Reynolds numbers over which perturbations are strongly damped and become rapidly negligible away from the source. This optimal (maximum) damping (see $\text{Im}(\unicode[STIX]{x1D701})$ in figure 3) creates two distinctive ‘quiet valleys’ in the $\tilde{R}e$ – $\tilde{x}$ maps (upstream and downstream of the source point) shown in figure 5. The position of each valley is governed by the intersection between the Stokes velocity-diffusion ( $\mathscr{U}^{\pm }$ ) and the Landau acoustic-damping ( $\mathscr{L}^{\pm }$ ) modes. The $\tilde{R}e$ intercept, $\tilde{R}e_{\text{I}}$ , can be seen from the asymptotic solutions to scale linearly with $\unicode[STIX]{x1D712}$ , i.e.  $\tilde{R}e_{\text{I}}\sim \unicode[STIX]{x1D712}$ if $\unicode[STIX]{x1D712}\gg 1$ . This is indeed visible from the approximate shift of the valley by approximately one decade in $\tilde{R}e$ for each increment in $\unicode[STIX]{x1D712}$ of one order of magnitude. However, the shift of the $\tilde{R}e$ – $\tilde{x}$ map with $\unicode[STIX]{x1D712}$ is not self-similar owing to the fixed ( $\unicode[STIX]{x1D712}$ independent) $\mathscr{P}^{\pm }$ and $\mathscr{S}$ branches. As $\unicode[STIX]{x1D712}$ is increased, the Landau entropy mode, $\mathscr{S}$ , interferes with the valley (on the downstream side only) and the ‘quiet’ region bears its signature: more and more pronounced density fluctuations within the valley as $\unicode[STIX]{x1D712}$ increases. In the process, notice that $\unicode[STIX]{x1D712}$ provides a way to align the $\mathscr{L}^{+}$ -mode damping with that of $\mathscr{S}$ . In the configuration shown here, this alignment occurs for $\unicode[STIX]{x1D712}\sim 10^{2}$ .

On the Stokes side of the valley (e.g. increasing $\unicode[STIX]{x1D707}\unicode[STIX]{x1D714}$ values), the transition into the $\mathscr{U}$ mode produces coherent structures on the dilatational rate extending many acoustic wavelengths ( $\unicode[STIX]{x1D706}_{0}$ ) away (see the coloured regions in figure 5 and Movie 1, which is available at https://doi.org/10.1017/jfm.2019.531). Based on the asymptotic results, the extent ( $\unicode[STIX]{x1D706}$ ) of these coherent structures scales according to $\unicode[STIX]{x1D706}/\unicode[STIX]{x1D706}_{0}\sim (\unicode[STIX]{x1D712}/\tilde{R}e)^{1/2}$ (for $\unicode[STIX]{x1D712}\gg 1$ ). This implies that high- $\unicode[STIX]{x1D712}$ fluid flows make the observation of a given large dilatational structure size (relative to the acoustic wavelength) possible at lower forcing frequencies (i.e. larger eddy sizes) than in the low- or zero- $\unicode[STIX]{x1D712}$ fluid. Superimposed to these coherent dilatational structures are coherent pressure–density fluctuations from the Newton mode, $\mathscr{N}$ , downstream of the source point. These structures will be shown to play a major role in shaping turbulence in high- $\unicode[STIX]{x1D712}$ fluids.

3.1 Computational strategy and initial conditions

Numerical integrations of the Navier–Stokes equations (2.1) are performed using the in-house solver CompReal from Touber & Alferez (Reference Touber and Alferez2019). The equations are integrated in time explicitly using a third-order total-variation-diminishing Runge–Kutta scheme (Shu & Osher Reference Shu and Osher1989). Spatial derivatives are evaluated with a fourth-order centred finite difference dispersion-relation-preserving scheme (Tam & Webb Reference Tam and Webb1993) using a thirteen-point stencil (a four points per wavelength cutoff is used for the dispersion optimisation process). The centred finite-difference scheme is stabilised using a centred thirteen-point eighth-order explicit discrete filter optimised in spectral space to minimise small-scale dissipation up to a cutoff wavenumber, following the methodology developed by Bogey & Bailly (Reference Bogey and Bailly2004) (a five points per wavelength cutoff is used). This strategy accurately propagates large-wavelength content whilst minimising numerical errors at small scales (both from a dispersion and dissipation point of view), enabling, in particular, the study of acoustic waves down to wavelengths five times the grid spacing. Solutions to the Euler equations (i.e. setting $\unicode[STIX]{x1D707}=\unicode[STIX]{x1D707}_{b}=\unicode[STIX]{x1D705}=0$ in (2.1)) are also considered. These form discontinuities (eddy shocklets) and a shock-capturing strategy based on explicit localised artificial diffusivity (Cook & Cabot Reference Cook and Cabot2004) is employed. Artificial bulk viscosity and thermal conductivity are added to the momentum and energy equations where shocks are detected following the approach by Kawai, Shankar & Lele (Reference Kawai, Shankar and Lele2010). Detected discontinuities are redistributed across approximately four grid cells with spurious noise levels of approximately 1 % of the ‘shock’ strength (when coupled with accurate finite-difference operators).

Simulations of decaying homogeneous turbulence are performed in a two-dimensional periodic domain of size $L\times L$ (where $L$ is an arbitrary length), discretised with $N\times N$ uniformly distributed points. The restriction to two-dimensional space is motivated by the associated reduction in computational cost (which, at best, scales linearly with $\unicode[STIX]{x1D712}$ ) but also by the associated removal of the vortex-stretching mechanism, which will prove useful when highlighting a mechanism by which vorticity is promoted by $\unicode[STIX]{x1D707}_{b}$ that is also present in three-dimensional space (see § 3.5). A pseudo-dimensionless velocity field, $\boldsymbol{u}_{o}$ , is derived by projecting its components drawn from a normal distribution onto a solenoidal field (based on the Helmholtz decomposition) and rescaled to the prescribed spectrum $E(\unicode[STIX]{x1D709})=\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D709}\exp [-(\unicode[STIX]{x1D709}/\unicode[STIX]{x1D709}_{0})^{2}/2]$ , where $\unicode[STIX]{x1D709}$ is the dimensionless wavenumber ( $\unicode[STIX]{x1D709}=L/\unicode[STIX]{x1D706}$ , where $\unicode[STIX]{x1D706}$ is the period of the wave of interest) and $\unicode[STIX]{x1D709}_{0}$ is a constant to prescribe. The constant $\unicode[STIX]{x1D6FD}$ is adjusted to make the pseudo velocity field unitary, i.e.  $\text{max}(\boldsymbol{u}_{o}\boldsymbol{\cdot }\boldsymbol{u}_{o})=1$ . A pseudo-dimensionless pressure field, $p_{o}$ , is then derived from $\boldsymbol{u}_{o}$ by solving the Poisson equation $\unicode[STIX]{x0394}p_{o}=-(\unicode[STIX]{x1D735}\boldsymbol{u}_{o})^{\intercal }\boldsymbol{ : }\unicode[STIX]{x1D735}\boldsymbol{u}_{o}$ (where ‘ $\boldsymbol{ : }$ ’ is the contracted product; the Laplacian, $\unicode[STIX]{x0394}$ , and gradient, $\unicode[STIX]{x1D735}$ , operators are rescaled by $L^{2}$ and $L$ , respectively). A target thermodynamic state $(\unicode[STIX]{x1D70C}_{0},p_{0})$ is selected (this choice is arbitrary in ideal gases) and the pseudo-fields rescaled according to $p=p_{0}+\unicode[STIX]{x1D70C}_{0}c_{0}^{2}\,p_{o}$ and $\boldsymbol{u}=c_{0}\boldsymbol{u}_{o}$ , where $c_{0}^{2}=\unicode[STIX]{x03B3}p_{0}/\unicode[STIX]{x1D70C}_{0}$ ( $c_{0}$ is the isentropic ideal-gas sound speed evaluated at the target thermodynamic state). The temperature field is obtained from the equation of state assuming $\unicode[STIX]{x1D70C}(\boldsymbol{x},t=0)=\unicode[STIX]{x1D70C}_{0}$ (where $\boldsymbol{x}$ is the position vector). Thus, the initial velocity fluctuations are at most sonic and the initial kinetic energy is solely supplied to the solenoidal field. Most of the initial energy is concentrated around length scale $\ell =L/\unicode[STIX]{x1D709}_{0}$ .

The initial condition is integrated in time using the Euler equations with shock capturing (E), the Navier–Stokes equations without (N) and with (B) bulk viscosity. Labels (E), (N) and (B) are used throughout the section to identify the set of equations used to provide the displayed result. Solutions using (N) and (B) are integrated only over a few energetic-eddy turnover times, whereas the Euler solution is integrated further.

3.2 Choice of grid, fluid and flow properties in light of the linear theory results

The strategy described above is applied to $N=4128$ , a choice which is representative of current high-end resolutions (i.e. number of point per characteristic length scale) achieved in three-dimensional direct numerical simulations of turbulence. The bulk (B) and no-bulk (N) viscosity cases are run at same Reynolds number, $\mathit{Re}\equiv \unicode[STIX]{x1D70C}_{0}U_{0}L/\unicode[STIX]{x1D707}=10^{4}\unicode[STIX]{x03C0}$ , with $U_{0}=c_{0}$ (i.e. the maximum speed at $t=0$ ). The bulk-to-shear viscosity ratio in case (B) is set to $\unicode[STIX]{x1D712}=4\times 10^{3}$ (recall that both $\unicode[STIX]{x1D707}$ and $\unicode[STIX]{x1D707}_{b}$ are kept constant and uniform).

Starting from the definition of the acoustic Reynolds number, $\tilde{R}e=\unicode[STIX]{x1D70C}_{0}c_{0}^{2}/(\unicode[STIX]{x1D707}\unicode[STIX]{x1D714})$ , and considering the dispersion relation $\unicode[STIX]{x1D714}=u_{\unicode[STIX]{x1D719}}\text{Re}(k)$ (with $\text{Re}(k)=2\unicode[STIX]{x03C0}\unicode[STIX]{x1D709}/L$ ), one arrives at

(3.1) $$\begin{eqnarray}\displaystyle \displaystyle \tilde{R}e=\frac{1}{2\unicode[STIX]{x03C0}\unicode[STIX]{x1D709}}\frac{\mathit{Re}}{\mathscr{M}_{\unicode[STIX]{x1D719}}\mathscr{M}_{0}}, & & \displaystyle\end{eqnarray}$$

where $\mathscr{M}_{0}\equiv U_{0}/c_{0}=1$ and $\mathscr{M}_{\unicode[STIX]{x1D719}}\equiv u_{\unicode[STIX]{x1D719}}/c_{0}$ . Assuming $\mathscr{M}_{\unicode[STIX]{x1D719}}\sim 1$ (e.g. acoustic waves in a zero mean flow), equation (3.1) gives $\tilde{R}e\sim \mathit{Re}/(2\unicode[STIX]{x03C0}\unicode[STIX]{x1D709})$ . The smallest computed (but not resolved) scale $\unicode[STIX]{x1D706}=L/N$ (giving $\unicode[STIX]{x1D709}=N$ ) provides a lower bound on the $\tilde{R}e$ scaling, i.e.  $\tilde{R}e\sim \mathit{Re}/(2\unicode[STIX]{x03C0}N)$ . Conversely, the largest computed scale $\unicode[STIX]{x1D706}=L$ (giving $\unicode[STIX]{x1D709}=1$ ) provides an upper bound, $\tilde{R}e\sim \mathit{Re}/(2\unicode[STIX]{x03C0})$ . For $\mathit{Re}=10^{4}\unicode[STIX]{x03C0}$ and $N=4128$ , the simulation is thus expected to explore approximately four decades of acoustic Reynolds numbers centred around $10^{2}$ , i.e.  $10^{0}\lessapprox \,\tilde{R}e\,\lessapprox 10^{4}$ .

The inherent finite range of resolved acoustic Reynolds numbers (from (3.1)) is put in regards to the linear theory results in figure 6 (and Movie 2) for both $\unicode[STIX]{x1D712}=0$ and $\unicode[STIX]{x1D712}=4\times 10^{3}$ (with $\unicode[STIX]{x03B3}=1.29$ , $\mathit{Pr}=0.7$ ). Base-flow speeds ranging from $\mathscr{M}=0$ to $1.2$ are shown to illustrate the sensitivity of the $\tilde{R}e$ – $\tilde{x}$ map to local flow-speed variations induced by large-scale eddies. The main idea here is that the quiet ‘valleys’ (§ 2.3) are positioned beyond the smallest resolved scales for the $\unicode[STIX]{x1D712}=0$ case but towards the largest resolved scales for the non-zero- $\unicode[STIX]{x1D712}$ case. Therefore, whilst the zero- $\unicode[STIX]{x1D712}$ flow should be under the influence of the Landau damping towards the small scales, the non-zero- $\unicode[STIX]{x1D712}$ flow should exhibit a behaviour consistent with that of the Stokes regime at these same length scales. Based on the linear theory, the Stokes regime is characterised by the emergence of fast (supersonic, see $\mathscr{U}^{\pm }$ modes) large-scale contraction/expansion waves (large relative to the Euler acoustic wavelength, $\unicode[STIX]{x1D706}_{0}$ , at the same forcing frequency) on top of which positively correlated pressure–density perturbations propagate at the local flow speed. As the flow Mach number increases, so does the wavelength of these $\unicode[STIX]{x1D70C}^{\prime }$ – $\,p^{\prime }$ (isothermal, see $\mathscr{N}$ mode) fluctuations, and the two structures form large (compared to $\unicode[STIX]{x1D706}_{0}$ ) regions of in-phase $\unicode[STIX]{x1D70C}^{\prime }$ – $p^{\prime }$ – $\unicode[STIX]{x1D703}^{\prime }$ fluctuations (see higher $\mathscr{M}$ cases for the $\unicode[STIX]{x1D712}=4\times 10^{3}$ case in Movie 2). Such structures largely differ from the $\unicode[STIX]{x1D70C}^{\prime }$ – $p^{\prime }$ – $\unicode[STIX]{x1D703}^{\prime }$ fluctuations induced by acoustic waves (where $\unicode[STIX]{x1D703}^{\prime }$ is always in phase opposition with respect to $\unicode[STIX]{x1D70C}^{\prime }$ and $p^{\prime }$ , see Movie 2, and the $\unicode[STIX]{x1D70C}^{\prime }$ – $\,p^{\prime }$ fluctuations are at constant entropy and not temperature). The range of Mach numbers given in figure 6 and Movie 2 match the range of local Mach numbers obtained in the present turbulence simulations. The particular phase alignment of the $\unicode[STIX]{x1D70C}^{\prime }$ – $p^{\prime }$ – $\unicode[STIX]{x1D703}^{\prime }$ fluctuations in the Stokes regime will be shown to play a crucial role in enhancing solenoidal structures at small scales.

Figure 6. Spatial ( $\tilde{x}$ ) evolution of $p^{\prime }$ over a range of acoustic Reynolds numbers ( $\tilde{R}e$ ) for $\unicode[STIX]{x1D712}=0$ (a) and $\unicode[STIX]{x1D712}=4\times 10^{4}$ (b) and $\mathscr{M}\in \{0,0.3,0.6,0.9,1.2\}$ (with $\unicode[STIX]{x03B3}=1.29$ , $\mathit{Pr}=0.7$ ). The lighting of the $p^{\prime }$ surface follows its signed magnitude and is coloured following $\unicode[STIX]{x1D703}^{\prime }$ – $\unicode[STIX]{x1D70C}^{\prime }$ correlations (see colour scale). The triangular markers position the 10 % drop in the asymptotic modes. See § 2.3 and figure 5 for details. Movie 2 gives a full forcing period. The light is deemed and the colouring off for $\tilde{R}e$ numbers a priori not captured by the current turbulence simulations. The $\unicode[STIX]{x1D712}>0$ simulation is designed to observe Stokes regime within the inertial range.

3.3 Selective removal and promotion of turbulence structures by $\unicode[STIX]{x1D707}_{b}$

Figure 7. Isotropic energy spectral density ( ${\mathcal{E}}$ ) of the total $(\text{t})$ , solenoidal $(\text{s})$ and dilatational $(\text{d})$ velocity fields of the Euler simulation over the integration time $t^{\ast }\equiv tU_{0}/L\in [0,60]$ . The initial solenoidal energy deposition is shown in black dash line. The vertical line labelled ‘Filt.’ and transparent white patch indicate the range of wavenumbers explicitly filtered by the numerical method (grid with $N=4128$ points). Both a direct enstrophy cascade ( ${\mathcal{E}}\sim \unicode[STIX]{x1D709}^{-3}$ for $\unicode[STIX]{x1D709}>\unicode[STIX]{x1D709}_{0}$ ) and inverse energy cascade ( ${\mathcal{E}}\sim \unicode[STIX]{x1D709}^{-5/3}$ for $\unicode[STIX]{x1D709}<\unicode[STIX]{x1D709}_{0}$ ) can be observed in the solenoidal field, with the direct cascade rapidly settling (times in light blue) compared to the inverse cascade (times in yellow). The dilatational field progressively assumes a Burgers-like spectrum ( ${\mathcal{E}}\sim \unicode[STIX]{x1D709}^{-2}$ ). Owing to the difference in the direct cascade and Burgers-like spectral exponents ( $-3$ versus $-2$ ), the solenoidal and dilatational energies becomes comparable at high wave numbers, where the classical exponents no longer hold. The total velocity displays a spectrum coinciding with the dominant contribution (i.e. mostly the solenoidal field here).

The initial condition described in § 3.1 is applied to $\unicode[STIX]{x1D709}_{0}=8$ (giving $\tilde{R}e\approx 6\times 10^{2}$ ). This choice enables the observation of a dual cascade on the solenoidal velocity field (see Mininni & Pouquet Reference Mininni and Pouquet2013) of the Euler flow, which is deliberately centred in the vicinity of the transition into the Stokes regime of the high-bulk fluid case (i.e.  ${\mathcal{E}}\sim \unicode[STIX]{x1D709}^{-3}$ for $\unicode[STIX]{x1D709}>\unicode[STIX]{x1D709}_{0}$ , and ${\mathcal{E}}\sim \unicode[STIX]{x1D709}^{-5/3}$ for $\unicode[STIX]{x1D709}<\unicode[STIX]{x1D709}_{0}$ , where ${\mathcal{E}}$ is the energy spectral density). Strictly speaking, the possibility of observing a double cascade is based on the existence of two quadratic invariants (i.e. the space-averaged kinetic energy and enstrophy, see Kraichnan Reference Kraichnan1967; Boffetta & Ecke Reference Boffetta and Ecke2012, for details) and is therefore restricted to incompressible (decaying) turbulence in the inviscid limit. Figure 7 gives the isotropic energy spectral density over 60 large-eddy turnover time units ( $U_{0}/L$ ) extracted from the velocity field $(\text{t})$ and its projection onto the solenoidal $(\text{s})$ and dilatational $(\text{d})$ components (following the Helmholtz decomposition). The solenoidal field is seen to follow the inviscid- and incompressible-flow theory, except for the noticeable decay in space-averaged kinetic energy (which is mostly converted here into internal energy by the artificial viscosity in eddy shocklets). The robustness of the incompressible-flow theory is rooted in the fact that most of the kinetic energy remains confined to the solenoidal modes (see Aluie Reference Aluie2011; Wang, Gotoh & Watanabe Reference Wang, Gotoh and Watanabe2017; Wang et al. Reference Wang, Wan, Chen and Chen2018, for examples), as can be seen from the dilatational kinetic energy, which assumes a Burgers-like spectrum (Gotoh Reference Gotoh1994) with energy levels well below that of the solenoidal field (except at high wavenumbers, here $\unicode[STIX]{x1D709}>3\times 10^{2}$ ).

The robust and rapid establishment (in a couple of $U_{0}/\ell$ time units) of the enstrophy cascade is exploited here to seed the kinetic energy at high wavenumbers (i.e.  $\unicode[STIX]{x1D709}>\unicode[STIX]{x1D709}_{0}$ , down to the filter scale) in a manner which is representative of (two-dimensional) turbulence. These scales are expected to be under the influence of the Stokes regime. Figure 8 and Movie 3 give the spatial (and time) variations of the thermodynamic and mechanical pressures, the density, the temperature, the dilatational rate and the vorticity. The colouring is such that correlations in density and thermodynamic pressure are highlighted (white/blue regions correspond to in-phase pressure–density increases/decreases). The temperature field is shown in regards to such correlation, based on the mechanical pressure. This is to discriminate isothermal compressions from the more usual isentropic compression, and therefore help to identify structures consistent with mode $\mathscr{N}$ . The dilatational rate and vorticity fields illustrate the spatial (and temporal) distribution of both the dilatational and solenoidal kinetic energies.

The Euler-based (E) flow is characterised by a vortex-filamentation process down to the resolved scales ( $-3$ enstrophy cascade), and the emission of acoustic waves over the same range of scales (which steepen to form shocks over time, $-2$ Burgers-like spectrum). Owing to the absence of transport coefficients, and except at shocks, fluid particles retain their initial entropy (set non-uniform, see § 3.1) and the temperature field captures both the filamentation process (advection by the local solenoidal flow) and the reversible local increases/decreases induced by passing acoustic waves (see Movie 3). Vortex cores are local pressure–density minima and induce stagnation points with local pressure–density maxima. The zero- $\unicode[STIX]{x1D712}$ (N) flow is similar to (E) and the two fields are in fact difficult to distinguish in the early times (see Movie 3). The departure from (E) is mostly visible at small scales, where the transport coefficients (set amongst the lowest possible for such computational grid) dissipate the kinetic energy: the vortex-filamentation process is limited to a physical scale by the shear viscosity ( $\unicode[STIX]{x1D707}$ ), and that of the temperature by thermal conductivity ( $\unicode[STIX]{x1D705}$ ); acoustic waves are increasingly damped towards small scales (expected Landau regime).

In contrast, the high- $\unicode[STIX]{x1D712}$ (B) flow is characterised by $\unicode[STIX]{x1D70C}^{\prime }$ – $p^{\prime }$ structures resembling those of the vortex filaments. A closer look reveals that such structures are not visible in the temperature field, which remains similar to that of (N) – $\unicode[STIX]{x1D707}_{b}$ does not act directly on the temperature field. These structures therefore share the properties of the Newton eigenmode ( $\mathscr{N}$ ). The $\mathscr{N}$ and $\mathscr{S}$ branches along the entropy path (see figure 3) intersect at $\tilde{R}e=[(4/3+\unicode[STIX]{x1D712})\unicode[STIX]{x03B3}/\mathit{Pr}]^{1/2}/\mathscr{M}$ , which for the highest $\mathscr{M}\sim 1.2$ values reached here translates to $\tilde{R}e\sim 70$ . Structures with sizes similar or smaller than that given for $\tilde{R}e=50$ in figure 8 are thus expected to be under the influence of the $\mathscr{N}$ mode (on the left of the quiet valley in figure 6). The dilatational field, although an order of magnitude smaller than for (N), is also found to follow the filamentary property of vorticity, and the vorticity field itself appears sharper than in (N). These remarks contradict the notion that more bulk viscosity necessarily imply that more small scales are damped. Acoustic waves are indeed absent in (B), see Movie 3, but other sharp structures emerge, which are borrowed from the Stokes regime (similarly to the notion that acoustic waves are borrowed from the Euler limit when dissipative processes become negligible).

Figure 8. Pressure ( $p^{\prime }\equiv (p-p_{0})/(\unicode[STIX]{x1D70C}_{0}c_{0}^{2})$ ) and density ( $\unicode[STIX]{x1D70C}^{\prime }\equiv (\unicode[STIX]{x1D70C}-\unicode[STIX]{x1D70C}_{0})/\unicode[STIX]{x1D70C}_{0}$ ) fluctuations (in blue–white) at $t^{\ast }\equiv tU_{0}/L=0.4$ for cases (E), (N) and (B), with $\tilde{R}e$ -scale estimates ((3.1) with $\mathscr{M}_{\unicode[STIX]{x1D719}}\mathscr{M}_{0}=1$ ). Narrow-band overlays give the mechanical pressure ( $p_{m}^{\prime }\equiv (p_{m}-p_{0})/(\unicode[STIX]{x1D70C}_{0}c_{0}^{2})$ ) and temperature ( $T^{\prime }\equiv (T-T_{0})/T_{0}$ ) fluctuations (in blue-red), the dilatational rate ( $\unicode[STIX]{x1D703}^{\prime }\equiv \unicode[STIX]{x03B8}L/c_{0}$ ) and (positive) vorticity ( $\unicode[STIX]{x1D714}^{\prime }\equiv (L/c_{0})(\unicode[STIX]{x1D735}\times \boldsymbol{u})\boldsymbol{\cdot }\boldsymbol{e}_{z}$ ) fields (in orange–white). Fields are tiled following the $L$ -periodicity (in both vertical and horizontal directions). Lighting on $p^{\prime }$ placed on the top-right corner is applied for depth effect (dip at vortex core). Movie 3 provides a time animation of the above for $t^{\ast }\in [0,0.5]$ . Key remarks and interpretations are given in the text.

Remarkably, the correlation between small-scale thermodynamic pressure structures and the dilatational field (expected from Stokes regime) gives the mechanical pressure field in (B) a topology that is similar to that of (E) and (N), see figure 8 and Movie 3, at the exception of acoustic-wave footprints, which are necessarily visible in the thermodynamic (and thus mechanical) pressure for (E) and (N) (but are absent in (B)). No footprint arising from Stokes regime is visible in $p_{m}$ for (B). The time rate of change of the specific kinetic energy of a fluid particle (derived from (2.1)) is governed by

(3.2) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}{\dot{e}}_{c}=-\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\left(p_{m}-{\textstyle \frac{1}{3}}\unicode[STIX]{x1D707}\unicode[STIX]{x03B8}\right)+\unicode[STIX]{x1D707}\unicode[STIX]{x0394}e_{c}, & & \displaystyle\end{eqnarray}$$

where $e_{c}\equiv \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{u}/2$ , ${\dot{e}}_{c}\equiv \unicode[STIX]{x2202}e_{c}/\unicode[STIX]{x2202}t+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})e_{c}$ , and assuming a uniform dynamic viscosity. With the mechanical pressure field deprived from both acoustic waves and Stokes structures, $p_{m}$ assumes a topology reminiscent of that of incompressible flows (i.e. low values at vortex cores, see figure 8 and Movie 3). Thus, particles orbiting within isolated vortices are expected to experience pressure forces, $\unicode[STIX]{x1D735}p_{m}$ , orthogonal to their own velocities, $\boldsymbol{u}$ , thereby producing no work, such that $e_{c}$ remains constant. Similarly, particles cruising in between vortices, such as a contrarotating vortices, would first experience a pressure gradient in their direction of motion (speeding it up) followed by the opposite process (sometime leading to a stagnation point).

This picture is to contrast with the no-bulk calculation which is very much under the influence of the work done by compression/expansion processes due to the passage of acoustic waves (leading to an increase in friction work and viscous diffusion). Thus, the kinetic energy of high- $\unicode[STIX]{x1D712}$ fluid flows is not directly drained through an excess dissipation in acoustic waves (as may have been intuited). Instead, no acoustic wave is formed (in Stokes regime) and the absence of the small-scale structures from the Stokes mode in the mechanical pressure confers incompressible-like properties to the kinetic energy. This comment applies to the translational kinetic energy of a fluid particle.

Figure 9. Vorticity enhancement along a vortex filament (see (3.4)). From left to right: $\unicode[STIX]{x1D6FA}$ (normalised by $U_{0}/L$ ) (a,e); $-\unicode[STIX]{x1D6FA}\unicode[STIX]{x03B8}$ (b,f), $\unicode[STIX]{x1D6F1}$ (c,g), $\unicode[STIX]{x1D708}\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}$ (d,h) (all normalised by $(U_{0}/L)^{2}$ ); for case $\unicode[STIX]{x1D712}=4\times 10^{3}$ (a,b,c,d) and case $\unicode[STIX]{x1D712}=0$ (e,f,g,h) $0.1\,U_{0}/L$ after initial release at $t=0$ on the $N=4128$ grid. The blue lines are streamlines with a thickness indicative of the velocity magnitude. Magenta contours give $\unicode[STIX]{x1D6FA}L/U_{0}=65$ . The grey scale is clipped according to the value at the bottom right corner of each image (with white for positive). In addition to the ‘filamentation’ process (responsible for the local growth in enstrophy), the vorticity-dilatation correlation ( $\unicode[STIX]{x1D6FA}\unicode[STIX]{x03B8}$ term) is seen to benefit from the Stokes-mode coherence to further contribute to the enstrophy production (in the no-bulk case, $\unicode[STIX]{x1D6FA}\unicode[STIX]{x03B8}$ is less structured owing to the presence of acoustic waves and $|\unicode[STIX]{x1D6FA}|$ is both enhanced and reduced within the filament). The intensity of the vorticity at small scales in the high-bulk fluid case is initially enhanced.

The picture for the rotational kinetic energy (enstrophy) is different. The time rate of change of the vorticity of a fluid particle (derived from (2.1)) is governed by

(3.3) $$\begin{eqnarray}\displaystyle \dot{\unicode[STIX]{x1D74E}}=(\unicode[STIX]{x1D74E}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}-\unicode[STIX]{x03B8}\unicode[STIX]{x1D74E}-\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}\times \boldsymbol{f}+\unicode[STIX]{x1D708}\unicode[STIX]{x0394}\unicode[STIX]{x1D74E}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D74E}\equiv \unicode[STIX]{x1D735}\times \boldsymbol{u}$ , $\dot{\unicode[STIX]{x1D74E}}\equiv \unicode[STIX]{x2202}\unicode[STIX]{x1D74E}/\unicode[STIX]{x2202}t+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\unicode[STIX]{x1D74E}$ , $\unicode[STIX]{x1D717}\equiv 1/\unicode[STIX]{x1D70C}$ , $\boldsymbol{f}\equiv \unicode[STIX]{x1D735}(p_{m}-\unicode[STIX]{x1D707}\unicode[STIX]{x03B8}/3)-\unicode[STIX]{x1D707}\unicode[STIX]{x0394}\boldsymbol{u}$ , $\unicode[STIX]{x1D708}\equiv \unicode[STIX]{x1D707}/\unicode[STIX]{x1D70C}$ , assuming uniform dynamic viscosity. The first term on the right-hand side of (3.3) is the vortex stretching term, which vanishes in two-dimensional space. The second term, $\unicode[STIX]{x03B8}\unicode[STIX]{x1D74E}$ , plays a similar role (i.e. both terms express the conservation of angular momentum in the absence of external forces) and is still present in two-dimensional space (this term is the fluid equivalent to the action of an ice skater drawing its arms inward when spinning). The third term is an extension to the baroclinic production term ( $-\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}\times \unicode[STIX]{x1D735}p$ ) and the last term corresponds to the diffusion of vorticity by the shear viscosity.

The contribution of each term on the right-hand side of (3.3) in the context of a two-dimensional vortex filament is given in figure 9. Considering the local Cartesian coordinates given in the figure, the projection of (3.3) on $\boldsymbol{e}_{z}$ gives

(3.4) $$\begin{eqnarray}\displaystyle \hspace{-8.0pt}\begin{array}{@{}ccccccccccc@{}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}}{\unicode[STIX]{x2202}t}} & + & u{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}}{\unicode[STIX]{x2202}x}} & + & v{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}}{\unicode[STIX]{x2202}y}} & = & -\unicode[STIX]{x1D6FA}\unicode[STIX]{x03B8} & +\,\,\unicode[STIX]{x1D6F1}\,\,+ & \unicode[STIX]{x1D708}{\displaystyle \frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D6FA}}{\unicode[STIX]{x2202}x^{2}}} & + & \unicode[STIX]{x1D708}{\displaystyle \frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D6FA}}{\unicode[STIX]{x2202}y^{2}}},\\[10.00002pt] O(\unicode[STIX]{x1D6F7}) & & \!O(\unicode[STIX]{x1D6F7}\unicode[STIX]{x1D6F6}\unicode[STIX]{x1D716}) & & \!O(\unicode[STIX]{x1D6F7})\! & & O(\unicode[STIX]{x1D6F7}) & & \!O(\unicode[STIX]{x1D6F7}\unicode[STIX]{x1D6F9}\unicode[STIX]{x1D716}^{2}) & & \!O(\unicode[STIX]{x1D6F7}\unicode[STIX]{x1D6F9})\end{array} & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FA}\equiv \unicode[STIX]{x1D74E}\boldsymbol{\cdot }\boldsymbol{e}_{z}$ , $\unicode[STIX]{x1D6F1}\equiv -[\unicode[STIX]{x1D735}\unicode[STIX]{x1D717}\times \boldsymbol{f}]\boldsymbol{\cdot }\boldsymbol{e}_{z}$ , $\unicode[STIX]{x1D6F7}\equiv U\dot{\unicode[STIX]{x1D6FF}}/\unicode[STIX]{x1D6FF}^{2}$ , $\unicode[STIX]{x1D6F6}\equiv U/\dot{\unicode[STIX]{x1D6FF}}$ , $\unicode[STIX]{x1D716}\equiv \unicode[STIX]{x1D6FF}/l$ , $\unicode[STIX]{x1D6F9}\equiv \unicode[STIX]{x1D707}/(\unicode[STIX]{x1D70C}_{0}\dot{\unicode[STIX]{x1D6FF}}\unicode[STIX]{x1D6FF})$ . The orders of magnitude are obtained considering $\unicode[STIX]{x1D6FA}\sim U/\unicode[STIX]{x1D6FF}$ , $\unicode[STIX]{x03B8}\sim \dot{\unicode[STIX]{x1D6FF}}/\unicode[STIX]{x1D6FF}$ , $u\sim U$ , $v\sim \dot{\unicode[STIX]{x1D6FF}}$ and a characteristic time $\unicode[STIX]{x1D6FF}/\dot{\unicode[STIX]{x1D6FF}}$ . Since the filamentation process is induced by the energy-rich eddies from the initial condition, the shearing and squeezing speeds are $U\sim U_{0}$ and $\dot{\unicode[STIX]{x1D6FF}}\sim U_{0}$ (hence $\unicode[STIX]{x1D6F6}\sim 1$ ). Moreover, let $\unicode[STIX]{x1D711}\equiv L/\unicode[STIX]{x1D6FF}$ , giving $\unicode[STIX]{x1D6F7}=\unicode[STIX]{x1D711}^{2}(U_{0}/L)^{2}$ and $\unicode[STIX]{x1D6F9}=\unicode[STIX]{x1D711}/\mathit{Re}$ .

Initially, the vortex filament is a blob of vorticity with $\unicode[STIX]{x1D6FF}\sim l\sim \ell$ , giving $\unicode[STIX]{x1D716}\sim 1$ , $\unicode[STIX]{x1D711}\sim \unicode[STIX]{x1D709}_{0}$ , and $\unicode[STIX]{x1D6F9}\sim \unicode[STIX]{x1D709}_{0}/\mathit{Re}\ll 1$ . Thus, the time rate of change of the vorticity of a fluid particle within the vortex filament is initially under the influence of the vorticity dilatation and production terms, i.e.  $\dot{\unicode[STIX]{x1D6FA}}\approx -\unicode[STIX]{x1D6FA}\unicode[STIX]{x03B8}+\unicode[STIX]{x1D6F1}$ . Assuming a constant compression rate $\unicode[STIX]{x03B8}<0$ , the homogeneous part of this last equation suggests that the vorticity of a fluid particle initially undergoes an exponential growth (for $\unicode[STIX]{x1D6FA}>0$ and $\unicode[STIX]{x03B8}<0$ as in the configuration of figure 9), i.e.  $\unicode[STIX]{x1D6FA}\sim \unicode[STIX]{x1D6FA}_{t=0}\exp (|\unicode[STIX]{x03B8}|t)$ (ignoring forcing from $\unicode[STIX]{x1D6F1}$ ).

At later times, when the filament geometry is such that $\unicode[STIX]{x1D6FF}\ll l$ , $l\sim \ell$ and $\unicode[STIX]{x1D716}\mathit{Re}\sim 1$ , giving $\unicode[STIX]{x1D711}=(L/l)\unicode[STIX]{x1D716}^{-1}\sim \unicode[STIX]{x1D716}^{-1}$ , and $\unicode[STIX]{x1D6F9}\sim (\unicode[STIX]{x1D716}\mathit{Re})^{-1}\sim 1$ , the time rate of change of vorticity of a fluid particle inside the filament ( $\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x2202}t+v\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x2202}y$ ) is under the competing influence of the enhancement due to the vorticity-dilatation correlation ( $-\unicode[STIX]{x1D6FA}\unicode[STIX]{x03B8}$ ) and viscous diffusion ( $\unicode[STIX]{x1D708}\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D6FA}/\unicode[STIX]{x2202}y^{2}$ ), again ignoring $\unicode[STIX]{x1D6F1}$ for now. This is essentially the time captured in figure 9. The filamentation ceases when the diffusion term eventually dominates every other term ( $\unicode[STIX]{x1D6F9}\gg 1$ ).

From these considerations on the vorticity equation, the vorticity-dilatation term can clearly play a significant role in the enhancement of $\unicode[STIX]{x1D6FA}$ during the early stages of the filamentation process (which is the main driver for the enstrophy cascade, see chapter 10 in Davidson Reference Davidson2004, for example). The high-bulk flow benefits from the strong correlation between $\unicode[STIX]{x1D6FA}$ and $\unicode[STIX]{x03B8}$ (which is related to Stokes regime). There is no definite vorticity-dilatation correlation in the no-bulk flow (which is entirely in the Landau regime) since the dilatational field is mostly due to acoustic waves travelling through the vortex filaments, thereby both promoting and reducing $\unicode[STIX]{x1D6FA}$ with no net effect (see figure 9 f). Term $\unicode[STIX]{x1D6F1}$ does not contribute much to the balance inside the filament shown in figure 9 in the high-bulk fluid case, and whilst more active in the no-bulk case, its contribution does not lead to a net effect (both signs occur within the structure).

3.4 Spectral signature

Figure 10. Isotropic energy spectral density at $t^{\ast }\in [0,0.3]$ for the solenoidal ( ${\mathcal{E}}_{s}$ ) and dilatational ( ${\mathcal{E}}_{d}$ ) velocity fields, and isotropic ‘energy’ spectral density for pressure ( ${\mathcal{E}}_{p}$ ), density ( ${\mathcal{E}}_{\unicode[STIX]{x1D70C}}$ ) and temperature ( ${\mathcal{E}}_{T}$ ) for case (B), coloured according to $\unicode[STIX]{x1D6E5}\equiv \log ({\mathcal{E}}^{\unicode[STIX]{x1D712}\neq 0})-\log ({\mathcal{E}}^{\unicode[STIX]{x1D712}=0})$ , with a lighting based on the slope of $\log ({\mathcal{E}}^{\unicode[STIX]{x1D712}\neq 0})$ with the light source placed above (darker implies steeper). The Euler (E) spectrum is shown in orange only at $t^{\ast }=0.3$ . The $\tilde{R}e$ axis is computed from (3.1) with $\mathit{Re}=10^{4}\unicode[STIX]{x03C0}$ , $\mathscr{M}_{\unicode[STIX]{x1D719}}\mathscr{M}_{0}=1$ . A dark ‘ribbon’ is applied where the spectrum reaches scales explicitly altered by the filter. Two red ‘ribbons’ mark where the entropy ( $\mathscr{E}$ ) and acoustic ( $\mathscr{A}^{+}$ ) paths reach their maximum damping (see figure 3) for $\mathscr{M}=1$ .

The statistical importance of the structures discussed previously are now assessed in the energy spectra. Figure 10 gives the isotropic spectral density of the velocity fields projected onto solenoidal ( ${\mathcal{E}}_{s}$ ) and dilatational ( ${\mathcal{E}}_{d}$ ) components over $2.4\ell /U_{0}$ from the initial solenoidal kinetic-energy injection at $\unicode[STIX]{x1D709}=8$ (downward triangular marker in ${\mathcal{E}}_{s}$ plot). The enstrophy ‘cascade’ is responsible for the rapid establishment of the ${\mathcal{E}}_{s}\sim \unicode[STIX]{x1D709}^{-3}$ spectrum in (E), which is eroded (with no clear spectral exponent) at small scales in (N) owing to the work done by friction ( $\unicode[STIX]{x1D707}$ acts on both shearing and compression/expansion motions), irreversibly converting the kinetic energy into internal energy (heat here). Some of the initial solenoidal kinetic energy is transferred to the dilatational kinetic energy following the emission of acoustic waves (and formation of eddy shocklets). Since the smallest resolved scales in (N) fall within the Landau regime (see figure 6), acoustic waves are increasingly damped at small scales and ${\mathcal{E}}_{d}$ rolls off before (in terms of $\unicode[STIX]{x1D709}$ values) reaching the scales explicitly removed by the filter.

In contrast, the high- $\unicode[STIX]{x1D712}$ case (B) exhibits a drop (relative to (N)) of two orders of magnitude in the dilatational kinetic energy for $\unicode[STIX]{x1D709}\in [10,80]$ , corresponding to a range of $\tilde{R}e$ values near the estimated optimal damping along the acoustic paths $\mathscr{A}^{\pm }$ (see § 2.3 and $\mathscr{A}^{+}$ marker in figure 10), where the transfer of the supplied solenoidal kinetic energy to the dilatational kinetic energy (via acoustic waves) is expected to be impaired. Consistent with this observation, scales coinciding with the emergence of Stokes regime experience easier transfers into dilatational kinetic energy (see the hump in ${\mathcal{E}}_{d}$ at $\tilde{R}e\sim 10^{2}$ ). These transfers are no longer acoustic in nature (as discussed previously) but are associated with the Stokes eigenmode, $\mathscr{U}$ , which prevails here by design. Together with the Newton mode, $\mathscr{N}$ (present so long a local flow exists), the flow is now characterised by in-phase $\unicode[STIX]{x1D70C}$ – $p$ – $\unicode[STIX]{x03B8}$ fluctuations which act to enhance the vorticity of vortex filaments (via vorticity-dilatation correlation, as discussed). This process increases the solenoidal kinetic energy contained at small scales (see the pronounced hump in ${\mathcal{E}}_{s}$ for (B), with energy levels exceeding those of (N) by more than three orders of magnitude).

Moving to the thermodynamic variables (figure 10), the pressure spectrum ( ${\mathcal{E}}_{p}$ ) develops a prominent hump with levels around the $\tilde{R}e=10$ mark exceeding those of (E) and (N) by four and five orders of magnitude, respectively. These scales correspond to the filamentary $p^{\prime }$ structures in figure 8 and Movie 3. Unlike for cases (E) and (N), the amplitude of the ${\mathcal{E}}_{p}$ spectrum at high wavenumbers in (B) is not acoustic in nature (see $\mathscr{P}$ and $\mathscr{N}$ eigenmodes), as the hump strikingly emerges on the Stokes side of the entropy path, $\mathscr{E}$ (see the marker indicating the optimal-damping position along the entropy path in figure 10). A similar observation applies to the density spectrum ( ${\mathcal{E}}_{\unicode[STIX]{x1D70C}}$ ), though not exceeding the energy levels from (E), consistently with the absence of $\unicode[STIX]{x1D70C}^{\prime }$ contributions on modes $\mathscr{P}^{\pm }$ . Remarkably, the temperature spectra ( ${\mathcal{E}}_{T}$ ) are nearly identical between (B) and (N) for an extended time (this is made visible by the absence of yellow/red colours), as was anticipated from Movie 3. Isothermal compressions are consistent with the Newton-eigenmode properties, characterised by in-phase $\unicode[STIX]{x1D70C}$ – $p$ fluctuations advected at the local flow speed and least damped at low $\tilde{R}e$ (the damping factor scales linearly with $\tilde{R}e$ , see (2.14) for details). At later times, a small hump develops on ${\mathcal{E}}_{T}$ as a consequence of the excess solenoidal energy at high wavenumbers. Indeed, vorticity filaments produce heat by friction work ( $\unicode[STIX]{x1D707}$ ) and the generated heat is diffused by thermal conductivity ( $\unicode[STIX]{x1D705}$ ), conferring the small hump in the temperature spectrum (i.e. inherited from the solenoidal kinetic-energy spectrum).

The above observations lead to a word of caution about numerical simulations of high- $\unicode[STIX]{x1D712}$ fluid flows. The choice of grid resolution ( $N$ ) and Reynolds number ( $\mathit{Re}$ ) is traditionally set so that the energy at small scales is drained by the shear viscosity and not the filter or the regularisation method employed (such solution is said to be grid converged). This is the case for (N) in figure 10 (naturally, Euler simulations are never grid converged and always stabilised artificially, here with shock capturing and filtering). However, when solving high- $\unicode[STIX]{x1D712}$ fluid flows, the setup must also take into account the transfer of solenoidal kinetic energy to the smallest computed scales (via the Stokes regime, as discussed in the present work). Evidence of the increased work done by the filter in case (B) is for example visible from the bright plateau past the filter scale in ${\mathcal{E}}_{s}$ (figure 10). Thus, not only the simulation time step must be adjusted to satisfy the viscous stability constraint, the mesh size must also be adjusted to capture this physical energy transfer to small scales. Both constraints have increased the simulation cost (in terms of time step size) of the present study by a factor $10^{4}$ when compared to a similar study looking only at (N) and (E). Thus, whilst simulations of decaying turbulence on $(4128)^{3}$ grid points are feasible on national computing facilities, the required number of time steps to cover a few eddy turnover times is still prohibitive.

3.5 Implications to three-dimensional turbulence

The ability of bulk viscosity to shape small-scale turbulence structures was demonstrated here in two-dimensional space, and the mechanisms at play based on a one-dimensional linear-mode analysis. The significance of these arguments to more practical three-dimensional flows is on a par with arguments relating to the relevance of acoustic waves in three-dimensional turbulent flows, which are derived from one-dimensional linear modes of strictly non-dissipative fluids. In the same spirit, the transition to Stokes regime, whilst anticipated in one dimension, remains relevant to multiple dimensions. This is demonstrated here by the ability to predict its occurrence in two-dimensional flows, despite fluctuation levels of the order of 10 % of the base-flow properties. Thus, eigenmodes such as $\mathscr{U}^{\pm }$ , $\mathscr{P}^{\pm }$ and $\mathscr{N}$ will survive a further dimensional increase, just like eigenmodes $a^{\pm }$ , $\unicode[STIX]{x1D70E}$ , $\mathscr{S}$ and $\mathscr{L}^{\pm }$ are still observed in three-dimensional flows (see Donzis & Jagannathan Reference Donzis and Jagannathan2016; Pan & Johnsen Reference Pan and Johnsen2017; Wang et al. Reference Wang, Wan, Chen and Chen2018, for examples). Additional modes emerge when increasing the dimensional space, such as shear (vorticity) waves, the bedrock of turbulence. In the nonlinear regime, these modes do interact, and the resulting flow can still be interpreted in the eigenmode basis (and more generally, the linear part of the Navier–Stokes operator is known to predict turbulence structures such as streaks, see Blesbois et al. Reference Blesbois, Chernyshenko, Touber and Leschziner2013, and references therein). In this work, a particular and significant (at least in two dimensions) interaction was shown to promote small-scale solenoidal motions, namely, the promotion of shear waves by dilatational modes (through the product of $\unicode[STIX]{x1D714}^{\prime }$ with $\unicode[STIX]{x1D703}^{\prime }$ , which displays high levels of correlation in the Stokes regime over many acoustic wavelengths). Whether this particular mechanism dominates over vortex stretching (see (3.3)) remains an open question, but its existence is certain, and the present work provides practical scalings, based on asymptotic solutions, to assess whether a given flow is capable of capturing the Stokes regime (by computing branch intersects in figure 3 for example) and expected associated properties. Whilst bulk viscosity is shown to facilitate the observation of the Stokes regime, by shifting it towards the inertial range, the regime itself is also present in zero-bulk fluids, and the results shown here may be of interest to studies dedicated to dissipative scales (see Khurshid, Donzis & Sreenivasan Reference Khurshid, Donzis and Sreenivasan2018, for example). Early attempts in turbulent Couette flows (Szemberg O’Connor Reference Szemberg O’Connor2018) at $\unicode[STIX]{x1D712}=10^{3}$ , where the kinetic energy is supplied at large scales and on the solenoidal modes only, indicate that the dilatational structures discussed here in two dimensions do form within vortex sheets and promote the solenoidal kinetic energy at small scales (unpublished work).