2.1. Symmetries

The above system has three types of symmetries associated with it, like in the Newtonian case (Chandler & Kerswell Reference Chandler and Kerswell2013). The shift-reflect symmetry maps

(2.8) \begin{align} \mathcal {S} [u, v, T_{xx}, T_{xy}, T_{yy}, p](x, y) \rightarrow [{-}u, v, T_{xx}, -T_{xy}, T_{yy}, p]({-}x, y+\pi ), \end{align}

where $\mathbf {u} = u \hat{\mathbf{x}} + v \hat{\mathbf{y}}$ and $\mathbf {T} = T_{xx} \hat{\mathbf{x}}\hat{\mathbf{x}} + T_{xy} (\hat{\mathbf{x}}\hat{\mathbf{y}} + \hat{\mathbf{y}}\hat{\mathbf{x}}) + T_{yy}\hat{\mathbf{y}}\hat{\mathbf{y}}$ . A reflection symmetry maps

(2.9) \begin{align} \mathcal {R}[u, v, T_{xx}, T_{xy}, T_{yy}, p](x, y) \rightarrow [u, -v, T_{xx}, -T_{xy}, T_{yy}, p](x, 2\pi -y). \end{align}

In addition to these two discrete symmetries, there is the continuous translational symmetry

(2.10) \begin{align} \mathcal {T}_s[u, v, T_{xx}, T_{xy}, T_{yy}, p](x, y) \rightarrow [u, v, T_{xx}, T_{xy}, T_{yy}, p](x+s, y) \end{align}

for $s \in [0, L_x)$ . Every solution is therefore associated with a set of solutions generated via these symmetries, and in particular the shift-reflect symmetry means that any solution moving in the positive $\hat{\mathbf{x}}$ direction has an associated solution moving in the negative $\hat{\mathbf{x}}$ direction.

2.2. Base flow

There is a one-dimensional base flow which depends only on $y$ which is

(2.11) \begin{align} \mathbf {U} = \cos y \, \hat{\mathbf{x}}, \quad P=0, \end{align}

(2.12) \begin{align} T_{xx} = \frac {W}{1+\varepsilon W} \left (1 - \frac {\cos 2y}{1+4\varepsilon W}\right ), \quad T_{xy} = \frac {-1}{1+\varepsilon W} \sin y \quad \textrm{and} \quad T_{yy} = 0. \\[6pt] \nonumber \end{align}

2.3. Linearising

To examine the linear stability of the base state, small perturbations proportional to $\textrm{e}^{ik(x-ct)}$ of all dependent variables are considered where $k\in \mathbb {R}$ is the wavenumber, and $c \in \mathbb {C}$ is an eigenvalue to be found. This results in variables of the form

(2.13) \begin{align} \mathbf {u} &= \mathbf {U}(y)+ \begin{bmatrix} u^{\prime}(y) \\ v^{\prime}(y) \end{bmatrix} \textrm{e}^{ik(x-ct)}, \quad p = p^{\prime}(y)\textrm{e}^{ik(x-ct)} \nonumber \\ \mathbf {T} &= \mathbf {T}(y) + \left[\begin{array}{l@{\quad}l} \tau _{xx}^{\prime}(y) & \tau _{xy}^{\prime}(y) \\[3pt] \tau _{xy}^{\prime}(y) & \tau _{yy}^{\prime}(y) \end{array} \right]\textrm{e}^{ik(x-ct)}, \end{align}

where a prime denotes a perturbative quantity that, due to the boundary conditions, must be $2\pi n$ periodic in $y$ . The imaginary part of the eigenvalue, $c_i$ , determines the linear stability of the system, with $\sigma := kc_i$ the growth rate. To reduce slightly the linearised equations which determine the evolution of the perturbations, we take the curl of (2.1) to eliminate the pressure, and use (2.4) to write all $\mathbf {C}$ in terms of $\mathbf {T}$ . Equations (2.1)–(2.4) then become

(2.14) \begin{align} ikRe \big[ (U-c)\big(D^2-k^2 \big) - D^2U\big]v^{\prime} = &-(1-\beta )\big [ -k^2 D\big(\tau _{xx}^{\prime} - {\tau _{yy}^{\prime}}\big) + ik \big(D^2 + {k^2}\big) \tau _{xy}^{\prime} \big ] \nonumber \\ &+ \beta \big( D^2 - {k^2}\big)^2 v^{\prime}, \end{align}

(2.15) \begin{align} \big [ik(U-c) + \frac {1}{W} - \varepsilon \big(D^2 - k^2 \big) \big ]\tau _{xx}^{\prime} & = -v^{\prime}DT_{xx} + 2ikT_{xx} u^{\prime} + 2T_{xy} Du^{\prime} + 2\tau _{xy}^{\prime}DU \nonumber\\ & \quad+ {\frac {2ik}{W}u^{\prime}} , \end{align}

(2.16) \begin{align} \big [ik(U-c) + \frac {1}{W}- \varepsilon \big(D^2 - k^2 \big)\big ]\tau _{xy}^{\prime} & = -v^{\prime}DT_{xy} + ikT_{xx} v^{\prime} + \tau _{yy}^{\prime}DU \nonumber \\ & + \frac {1}{W} \big(Du^{\prime} + {ikv^{\prime}}\big) , \end{align}

(2.17) \begin{align} \big [ik(U-c) + \frac {1}{W}- \varepsilon \big(D^2 - k^2 \big)\big ]\tau _{yy}^{\prime} = 2ikT_{xy} v^{\prime} + \frac {2}{W}Dv^{\prime} , \end{align}

(2.18) \begin{align} iku^{\prime} + Dv^{\prime} = 0 ,\end{align}

where $D:= \textrm{d}/\textrm{d}y$ . The costly procedure of solving this eigenvalue problem over the whole domain $y\in [0,2\pi n]$ can be avoided by applying Floquet analysis just across one forcing wavelength $y\in [0,2\pi ]$ and including a modulation parameter $\mu$ to compensate. A mode with Floquet exponent $\mu$ has perturbations of the form $\phi ^{\prime} = \hat \phi (y) \textrm{e}^{i\mu y}$ , where $\phi ^{\prime}$ is the perturbation of any flow variable and $\hat \phi$ is $2\pi$ periodic. The resultant perturbation has periodicity $2\pi / \mu$ when $1/\mu \in \mathbb {N}$ , as all base flow quantities have the same periodicity as the forcing. Values of $1/\mu$ which factor into $n$ then satisfy the periodic boundary conditions over the large domain of $y\in [0,2\pi n]$ .

Figure 2.

The eigenvalue spectrum when $E=81$ , $Re=2$ , $\beta =0.95$ , $\varepsilon =0$ , $\mu =0$ and $k=0.2$ , with resolution $N_y=300$ (blue circles) and $N_y=400$ (red dots), where $N_y$ is the number of Chebyshev modes considered in the eigenvalue problem. The centre mode has unstable eigenvalues at $\text{c}=\pm 1.01016736+0.05925964i$ , while a stable continuous spectrum is seen with $c_i\lt 0$ . Insets show the polymer stress trace (colours) and streamfunction (contours) of an eigenmode $\phi$ , alongside symmetries of $\phi$ . While reflections $\mathcal {R}$ and translations $\mathcal {T}_s$ leave the eigenvalue $c=c_r+ic_i$ unchanged, shift-reflections $\mathcal {S}$ produce modes with eigenvalue $-c_r+ic_i$ that travel in the opposite direction.

3.1. Harmonic disturbances ( $n=1$ )

Kerswell & Page (Reference Kerswell and Page2024) show that there is an unstable eigenfunction when $n=1$ in ultra-dilute vKf that resembles the centre-mode eigenfunction in a channel. Here we go a step further and show that the $(Re,k)$ neutral curves show the distinctive centre-mode loops seen in channel flow (as shown in figure 10 of Khalid et al. (Reference Khalid, Chaudhary, Garg, Shankar and Subramanian2021a)), and they follow the same scaling relation for $E \ll 1$ . However, we also show that the behaviour for large $E$ is substantially different in this unbounded flow from that of the centre mode in channel flow, as the instability is not suppressed as $E\rightarrow \infty$ . Instead, the instability exists in inertialess vKf.

To consider the system with $n=1$ we take the linearised system with Floquet exponent $\mu = 0$ . We plot an example eigenvalue spectrum in figure 2, showing how an eigenmode and its eigenvalue $c=c_r+ic_i$ are affected by the symmetries discussed in § 2.1. The shift-reflect symmetry $\mathcal {S}$ generates a new eigenmode with eigenvalue $c=-c_r+ic_i$ that travels in the opposite direction. The symmetries of the centre mode in the $n=1$ case make the eigenmode invariant under reflection $\mathcal {R}$ . Lastly, translational symmetries $\mathcal {T}_s$ correspond to a phase shift of the mode. The stability of all such symmetries is the same, as the growth rate is unchanged under all operators.

Figure 3.

$(a{-}d)$ The centre-mode neutral curves in the $(Re, k)$ plane for $\beta =0.95$ , $\varepsilon =0$ , $\mu =0$ and $(a, b)$ $E = 0.3, 0.6, 1.0, 1.8, 3.1, 5.5, 9.5$ (light to dark) and $(c, d)$ $E=81, 107, 142, 187, 272, 359, 475$ (light to dark). Note that $(b,d)$ are scaled versions of $(a,c)$ , respectively, demonstrating that for small $E$ , $Re_{crit} \sim E^{-3/2}$ , while for large $E$ , $Re_{crit} \sim E^{-1}$ . We plot eigenfunctions with $k=0.2$ and (e) $E=3.1, Re=300$ , (f) $E=81, Re=2$ and (g) $E=81, Re=20$ . These correspond to the instability in the low- $E$ regime, the main loop in the high- $E$ regime and the secondary loop in the high- $E$ regime, respectively. Colours show the polymer stress trace field, while contours show the streamfunction. This figure demonstrates that the elastic instability seen at high $\beta$ and low $E$ is the centre mode, and that a different scaling regime exists at high $E$ .

Figures 3(a) and 3(b) show unscaled and scaled $(Re,k)$ neutral curves for small values of $E$ . The neutral curves take the form of loops at small $E$ , as is the case in channel flow. For these loops, $Re_{crit} \sim E^{-3/2}$ and $k_{crit} \sim E^{-1/2}$ , where $Re_{crit}$ denotes the smallest Reynolds number on the neutral curve and $k_{crit}$ denotes the wavenumber at that point. These match the scalings for the centre mode in channel flow (Khalid et al. Reference Khalid, Chaudhary, Garg, Shankar and Subramanian2021a ). An eigenfunction in this regime is plotted in figure 3(e), resembling that of the channel flow centre mode as seen in figure 5 of Khalid et al. (Reference Khalid, Chaudhary, Garg, Shankar and Subramanian2021a ). The scaling relation and the eigenfunction both suggest that this instability is the centre mode.

Figures 3(c) and 3(d) show the neutral curves for large $E$ . The behaviour here is substantially different from that of channel and pipe flow, where at sufficiently large $E$ the centre mode is suppressed (Khalid et al. Reference Khalid, Chaudhary, Garg, Shankar and Subramanian2021a ; Chaudhary et al. Reference Chaudhary, Garg, Subramanian and Shankar2021). Here in vKf, the instability exists at all $E$ (only up to $E=475$ is shown), and we see that the loops have $Re_{crit} \sim E^{-1}$ and $k_{crit} \sim E^0$ as $E \rightarrow \infty$ . Equivalently one can consider the scalings in terms of $Re$ and $W$ to obtain that as $Re \rightarrow 0$ for fixed $W$ , $W_{crit} \sim Re^0$ and $k_{crit} \sim Re^0$ , meaning these critical parameters are independent of $Re$ . While $Re_{crit}$ and $k_{crit}$ are on a loop with these scalings, a secondary loop also exists at larger $Re$ . These secondary loops collapse as $E$ increases. It is worth pointing out that in this regime, the largest streamwise wavenumber at which instability is seen is $k \approx 0.9$ , and hence no linear instability is seen when simulating a box with $L_x=2\pi$ and $n=1$ , as is true in the Newtonian case (Marchioro Reference Marchioro1986). The eigenfunctions for large $E$ are seen in figures 3(f) and 3(g), showing the main loop and the secondary loop, respectively. The eigenfunctions on the secondary loop visually resemble that of the main loop, though activity is more spread out across the flow and less clearly localised to $y=\pi$ . Both visibly resemble the centre-mode eigenfunction in the low- $E$ regime in figure 3(e). To further demonstrate that the high- $E$ instabilities are the centre mode, figure 4 shows that the neutral curve loops from the high- $E$ regime continuously track into the loops from the low- $E$ regime. This means that the centre mode is responsible for the instability across all values of $E$ .

Figure 4.

The centre-mode neutral curves in the $(Re, k)$ plane for $\beta =0.95$ , $\varepsilon =0$ , $\mu =0$ and $E = 9.5, 17, 41, 81$ (dark to light). Secondary loops exist for the $E=41, 81$ curves. This shows that the neutral curve loops from the low- $E$ regime ( $E\lt 9.5$ ) can be continuously tracked into the main loops in the high- $E$ regime ( $E\gt 81$ ).

The neutral curves in the $(Re, E)$ and $(Re, W)$ planes are shown in figure 5. These demonstrate that the centre-mode instability is linearly unstable at vanishing $Re$ in vKf over a range of concentrations $\beta$ . This contrasts with the cases of inertialess pipe flow, which is linearly stable (Chaudhary et al. Reference Chaudhary, Garg, Subramanian and Shankar2021), and inertialess channel flow, which is only linearly unstable for ultra-dilute fluids with $\beta \gt 0.9905$ (Khalid et al. Reference Khalid, Chaudhary, Garg, Shankar and Subramanian2021a ,Reference Khalid, Shankar and Subramanian b ).

The neutral curves also show a second instability, which exists at vanishing $W$ and is inertial in nature. This instability is seen in Newtonian Kolmogorov flow, with a purely imaginary eigenvalue $c$ (Meshalkin & Sinai Reference Meshalkin and Sinai1961) and we identify an instability as inertial here if it similarly has zero frequency. This instability is suppressed by elasticity, with there being a maximum $E$ at which it exists.

Figure 5.

The neutral curves across $k\in \mathbb {R}$ in the (a) $(Re, E)$ and (b) $(Re, W)$ planes when $\varepsilon =0$ , $\mu =0$ and $\beta =0.5, 0.8, 0.9, 0.95$ (light to dark). This demonstrates that the centre mode exists in the inertialess system across a range of $\beta$ . Eigenfunctions for parameters on the neutral curves are shown in (c) when $(\beta , Re, W, k)=(0.5, 0.5, 5.78, 0.47)$ (blue circle) and (d) when $(\beta , Re, W, k)=(0.95, 0.5, 28.7, 0.60)$ (black square). Colours show the polymer stress trace field, while contours show the streamfunction.

3.2. The absence of PDI

The linear stability analysis in § 3.1 was performed with $\varepsilon =0$ , and hence PDI is absent. However, to run simulations, introducing a finite $\varepsilon$ could potentially introduce PDI, as is the case for wall-bounded flows. We check that PDI is not in this system by considering how the neutral curves in the $(W,k)$ plane are affected by introducing finite $\varepsilon =10^{-3}$ in figure 6. Wavenumbers up to $k=100$ were considered. The centre-mode loops that exists with vanishing $\varepsilon$ are adjusted slightly, but no new modes of instability are identified, meaning PDI was not found in vKf. This was true when $\beta =0.95$ , meaning the simulations in § 4 and 5 do not contain PDI. In addition, we checked for PDI in a more concentrated fluid ( $\beta =0.2$ ), where PDI was demonstrated to be particularly unstable in bounded flows (Lewy & Kerswell Reference Lewy and Kerswell2024). The lack of PDI in Kolmogorov flow is consistent with Lewy & Kerswell (Reference Lewy and Kerswell2024) which considers $\beta \ll 1$ and suggests that PDI requires boundaries to exist in such a fluid. This finding confirms that PDI is not necessarily seeded at regions of maximal base shear when polymer stress diffusion is present, as was the case in the wall-bounded rectilinear flows.

Figure 6.

The neutral curves for non-zero-frequency modes in the $(Re, k)$ plane when $\mu =0$ , $\varepsilon =0$ (blue solid lines) and finite $\varepsilon =10^{-3}$ (red dotted lines) when (a) $\beta =0.95$ and $E = 8, 16, 32, 256, 512$ (light to dark) and (b) $\beta =0.2$ and $E=0.5, 1, 2, 8, 64, 256$ (light to dark). Wavenumbers as high as $k=100$ were considered. These demonstrate that PDI was not identified in Kolmogorov flow, and that finite $\varepsilon$ generally stabilises the centre-mode instability.

Figure 7.

$(a)$ The most unstable Floquet modes in the $(Re, W)$ plane for $\beta =0.95$ , with $\varepsilon =0$ and instabilities over wavenumbers $k \in \mathbb {R}$ are considered, with colour denoting which Floquet mode is most unstable. Colours correspond to $\mu =0$ (blue), $\mu =1/2$ (orange), $\mu =1/3$ (green), $\mu =1/4$ (cyan), $\mu =1/5$ (red), $\mu =1/6$ (brown), $\mu =1/7$ (pink). $(b)$ The same on a log scale. Eigenfunctions are plotted with parameters $(c)$ $(W, Re, k, \mu ) = (1, 10, 0.5, 0)$ , $(d)$ $(W, Re, k, \mu ) = (20, 0.5, 0.5, 1/2)$ and $(e)$ $(W, Re, k, \mu ) = (600, 20, 0.01, 1/7)$ . $(f)$ The maximum growth rate $\sigma ^*$ of each Floquet mode (same colours as in [ $a$ ]) across all $k\in \mathbb {R}$ for $Re=0$ , $\beta =0.95$ and $\varepsilon =0$ as $W$ varies. These plots demonstrate that all elastic instabilities are the centre mode, which is generally most unstable when $\mu =1/2$ , while the inertial instability is most unstable when $\mu =0$ .

3.3. Modulated disturbances ( $n\gt 1$ )

In this section we consider linear stability when $n\gt 1$ by considering modes with Floquet exponents $\mu$ satisfying $1/\mu \in \mathbb {N}$ . Figure 7 shows the most unstable Floquet modes in the $(Re, W)$ plane for a viscosity ratio of $\beta =0.95$ . All wavenumbers $k\in \mathbb {R}$ are considered, as well as the first seven Floquet modes in figures 7(a) and 7(b), which focus on the behaviour of the system at small and large $(Re, W)$ , respectively. The orange regions demonstrate that it is the $\mu =1/2$ Floquet mode that generally makes the centre-mode instability most unstable, corresponding to perturbations with wavelengths of $4\pi$ in the $y$ direction or double the forcing wavelength. These are therefore ‘subharmonic’ disturbances as they have half the spatial frequency of the forcing. The blue region shows that the inertial (Newtonian) instability is most unstable when $\mu =0$ , when there is no modulation and perturbations have the same wavelength as the forcing, i.e. ‘harmonic’ disturbances. The neutral curve in figure 7(a) resembles figure 1 in Boffetta et al. (Reference Boffetta, Celani, Mazzino, Puliafito and Vergassola2005), which considered the stability of a system equivalent to $n=64$ and $k\in \mathbb {N}/64$ . They identified the distinct elastic and inertial instabilities, but our plot allows us to in addition identify the most unstable wavelength of perturbation. It is generally the subharmonic that determines the linear stability of the centre mode. The neutral curves of figure 7(b) are zoomed out and use a log scaling which reveals that there is a part of the parameter space where even lower-order harmonics are most unstable.

Figure 8.

$(a)$ The maximum growth rate $\sigma ^*$ and most unstable wavenumber $k^*$ as $\beta$ varies when $Re=0$ , $\varepsilon =0$ , $\mu =0$ and $W=40, 80, 160$ (light to dark). Asymptotics derived in Appendix A are shown by the black dotted lines. Most unstable eigenfunctions are shown for $W=160$ and (b) $\beta =0$ and (c) $\beta =0.95$ with colours showing the polymer stress trace field and contours showing the streamfunction. ( $d)$ Plots of $\sigma ^*$ and $k^*$ in the inertialess UCM fluid for various Floquet modes with $\beta =0$ , $Re=0$ , $\varepsilon =0$ and $\mu =0, 1/2, 1/3, \ldots , 1/7$ with colours as in figure 7. The asymptotic limits as $W\rightarrow \infty$ are shown by horizontal black dashed lines. When $\mu =0$ , $W\sigma ^* \rightarrow 0.784$ and $Wk^* \rightarrow 1.526$ , while when $\mu \gt 0$ , $W\sigma ^* \rightarrow 1.139$ and $Wk^* \rightarrow 1.764$ . The centre mode is therefore generic across $\beta$ , existing even in the UCM fluid, and $k^* \sim W^{-1}$ and $\sigma ^* \sim W^{-1}$ .

The $\mu =0$ inertial instability is shown in figure 7(c). Its trace field is antisymmetric about lines of peak base velocity, and it has zero frequency. The $\mu =1/2$ centre mode is shown in figure 7(d), and is clearly a modulated version of the $\mu =0$ centre mode shown in figure 5(d). Similarly, we plot the preferred mode when lower-order harmonics of the centre mode are most unstable in figure 7(e). This suggests that all elastic instabilities are centre modes for all Floquet modes, not just when $\mu =0$ .

The dependence of the centre-mode growth rate on $\mu$ is shown in figure 7(f). We consider vanishing inertia ( $Re=0$ ), and consider the maximum growth rate across all wavenumbers, $\sigma ^* := \max_{k \in \mathbb {R}} \sigma$ , as $W$ varies for fixed $\mu$ ( $k^*$ is defined as the most unstable wavenumber). We see that the harmonic disturbances ( $\mu =0$ ) are the most stable of those plotted. The subharmonic with modulation $\mu =1/2$ is the most unstable, and then subsequent modulations increase in stability.

3.4. The inertialess centre mode in the UCM fluid

So far we have mainly limited our results to the specific choice of $\beta =0.95$ . At this dilute concentration, the elastic instabilities in vKf for both $n=1$ and $n\gt 1$ have been identified as the centre mode. We now demonstrate that the centre mode exists not only when $\beta \sim 1$ , but for all $\beta$ . It even exists when $\beta =0$ and the model reduces to the UCM fluid.

We plot in figure 8(a) the maximum growth rate $\sigma ^*$ and the most unstable wavenumber $k^*$ of the inertialess instability as $\beta$ varies for the $n=1$ system. This shows that there is a smooth continuation of the centre mode at $\beta =0.95$ to the instability seen at lower  $\beta$ , suggesting that the instability seen at low $\beta$ is also the centre mode. Eigenfunctions are plotted at both $\beta =0$ and $\beta =0.95$ in figures 8(b) and 8(c) and clearly resemble each other, again suggesting that the instability at $\beta =0$ is the centre mode. The centre mode is not suppressed by low $\beta$ in vKf.

Figure 8(d) shows the behaviour of the various Floquet modes in the inertialess UCM fluid. As $W\rightarrow \infty$ , all harmonics tend to the same growth rate that is more unstable than the $\mu =0$ mode. Of these harmonics, the $\mu =1/2$ mode is the most unstable, as was the case when $\beta =0.95$ with vanishing inertia (shown in figure 7 g). The subharmonic is therefore most unstable for both high and low $\beta$ .

These plots also demonstrate that the correct asymptotic scalings for the inertialess centre mode are $k^* \sim W^{-1}$ and $\sigma ^* \sim W^{-1}$ as $W\rightarrow \infty$ . This is shown across all $\beta$ when $\mu =0$ (see figure 8 a) and various $\mu$ when $\beta =0$ (see figure 8 d). The equations governing the inertialess asymptotic limit of $W \rightarrow \infty$ are identified in Appendix A, and these asymptotics are plotted in both figures 8(a) and 8(d), showing their validity. The centre mode is therefore present in a very simple fluid when $Re=\beta =\varepsilon =0$ and $W \rightarrow \infty$ , demonstrating that while elasticity is required for the instability to exist, inertia, viscosity and polymer diffusion are not.

3.5. Relaminarisation in the $W\rightarrow \infty$ limit

Returning to the Oldroyd-B model, the flow becomes linearly stable as $W\rightarrow \infty$ for any fixed domain length. This is due to the centre mode only being unstable to a pocket of wavenumbers that scale like $k \sim W^{-1}$ , as suggested by the asymptotic analysis performed in Appendix A. Hence, at sufficiently large $W$ , all unstable wavelengths are longer than the channel itself, meaning the system is not susceptible to the centre-mode instability.

It will be useful to define $L_x^{min} := 2\pi / k_{max}$ , which is the shortest domain in which instability can be found, i.e. the system is stable when $L_x \lt L_x^{min}$ . When $Re=0$ , figure 9 demonstrates that $L_x^{min} \sim W$ , and hence that as $W \rightarrow \infty$ , any finite $L_x$ eventually becomes stable. Figure 9 also shows that the order in which the Floquet modes stabilise as $L_x/W \rightarrow 0$ depends on the concentration $\beta$ . This limit is achieved for fixed $L_x$ as $W \rightarrow \infty$ . In this limit, when $\beta \gt 0.62$ the $\mu =1/2$ mode stabilises before the $\mu =0$ mode, while the opposite is true for $\beta \lt 0.62$ . This is consistent with figure 7(f) in which $\beta =0.95$ , and the $\mu =1/2$ mode stabilises before the $\mu =0$ mode for negligible inertia as $W\rightarrow \infty$ . These results remain qualitatively true when inertia is introduced, with $L_x^{min} \sim W$ when $Re=1, 10$ (not shown). Relaminarisation therefore occurs as $W\rightarrow \infty$ both with and without inertia.

Figure 9.

The minimum $L_x$ at which laminar flow becomes unstable to perturbations with $\mu =0$ (blue) or $\mu =1/2$ (orange), when $Re=0$ , $\varepsilon =0$ and $W=50, 100, 200, 500, 1000$ (light to dark). Black dashed lines correspond to the asymptotic limit described in Appendix A, and they cross over at $\beta =0.62$ . The only region which is stable as $W\rightarrow \infty$ is shaded in red. This confirms that as $W \rightarrow \infty$ , $L_x^{min} \sim W$ , meaning only very long channels are linearly unstable for large $W$ . For $\beta \lt 0.62$ , the $\mu =0$ harmonic stabilises before the $\mu =1/2$ subharmonic as $L_x/W\rightarrow 0$ , while the opposite is true when $\beta \gt 0.62$ .

Figure 10.

Bifurcation plots for $\beta =0.95$ , $Re=0.5$ , $\varepsilon =10^{-3}$ and $L_x=4\pi$ . These show the deviation of the volume-averaged (a) trace $\varSigma$ and (b) kinetic energy $K$ from the laminar state, which has trace and kinetic energy $\varSigma _0$ and $K_0$ . We show stable solutions in both the $n=2$ system and the $n=1$ system. Blue corresponds to travelling-wave solutions, red to equilibria and green to limit cycles. The polymer stress trace of the solution at each of the six symbols are shown in figure 12. Bifurcation points (BP) due to the linear instability are shown with black crosses at $W=33, 86$ when $n=1$ and black pluses at $W=15, 76$ when $n=2$ .

5.1. The centre mode is responsible for transition to turbulence

Motivated by the presence of the centre-mode arrowhead in many of the coherent structures in § 4, we now demonstrate that the centre-mode instability is the cause of turbulence in this system. In fact, in the absence of any other elastic instabilities, all the dynamics ultimately comes from the centre-mode instability.

We first produce a state diagram in figure 13(a), which shows what types of final states are reached when simulations are initialised with finite-amplitude disturbances as the Weissenberg number $W$ and channel length $L_x$ are varied. Each state was found by taking the laminar base flow and adding a random perturbation to the first 50 % of the Fourier modes in fields $u$ , $T_{xx}$ and $T_{xy}$ , with an amplitude of 10 % of the mean laminar base field. All remaining initial fields were constructed from these. We simulated each pair of parameters ( $W$ , $L_x$ ) five times, and identify the final states as either laminar (blue circle), a travelling wave (green square), a periodic orbit (yellow cross), a quasi-periodic orbit (black triangle) or chaotic (purple star) using the procedure discussed in Appendix B. In regions of linear stability, this protocol produced exclusively laminar solutions, so to demonstrate subcriticality we simulated an additional 10 runs as well. In these, random perturbations were added to the first 5 % of Fourier modes with an amplitude of 20 % of the laminar base field for five runs, and an amplitude of 10 % of the laminar base field was set in the other five. This focus on longer-wavelength perturbations was sufficient to excite subcritical solutions; however, not all subcritical solutions were found, as a subcritical travelling-wave solution that is known to exist at $W=13, L_x=4\pi$ (see figure 10) was not found here. For some parameters these runs identified multiple types of final states, which are indicated as overlaid symbols in figure 13. These simulations used a resolution of $(N_x, N_y)=(128, 192)$ Fourier modes with a time step of $\textrm{d}t=8\times 10^{-3}$ and were confirmed using $(N_x, N_y)=(256, 256)$ .

Non-laminar solutions were only seen in regions of parameter space close to where the centre mode is linearly unstable, suggesting that these states originate from the centre-mode bifurcation. We plot all distinct final states in figure 13(b). While some states look similar (e.g. those at $W=80$ , $L_x=8\pi$ ), kinetic energy time series differ (not shown).

For the parameters considered, the shortest channel in which chaotic behaviour was found was $L_x=6\pi$ , at $W$ corresponding to the low-shear region of where the centre mode is linearly unstable. To check that the chaos was self-sustaining and not transient, we ran all chaotic simulations for $10^5$ time units of $2\pi U_0 / L_0$ , and saw no collapse in the time series of $K$ .

To examine the transition scenario, the turbulent state at $W=20$ and $L_x=6\pi$ was taken as a starting point and then $W$ lowered adiabatically. Figure 14 shows how the flow transitions from one where the time series of $K$ is chaotic to one that is periodic via regions of quasi-periodicity as $W$ is decreased, before eventually becoming a constant when the state is a travelling wave. At $W=18$ , the onset of turbulence is marked by the presence of bursting events that temporarily increase $K$ seemingly at random, before returning to periodic behaviour. When $W$ is decreased further, periodic behaviour (e.g. at $W=17$ ) and quasi-periodicity (e.g. at $W=15$ ) are found, until a travelling-wave state is realised at $W=7$ . The periodic and travelling-wave solutions closely resemble the centre-mode arrowhead.

Figure 14.

Kinetic energy time series of solutions (left), with the trace field at $t=5000$ (right) as $W$ is lowered from $W=20$ . The trace field colours use the same scale as figure 18 for comparison purposes. Parameters are $\beta =0.95$ , $Re=0.5$ , $\varepsilon =10^{-3}$ , $n=2$ and $L_x=6\pi$ . Symbols, as in figure 13(a), show chaos (purple star), quasi-periodic orbits (black triangle), periodic orbit (yellow cross) and a travelling wave (green square). This shows that the turbulent state at $W=20$ is connected to states that strongly resemble the centre-mode arrowhead.

Figure 15.

Power spectra for the states from figure 14 with $W=15$ and $W=10$ . This shows that the state at $W=15$ is truly a quasi-periodic orbit with discrete incommensurate frequencies, while the state at $W=10$ has a broad band of frequencies.

We plot the frequency spectra of the quasi-periodic behaviour at $W=15$ and $W=10$ in figure 15, where

(5.1) \begin{align} S_K(\omega ) = \left |\frac {1}{T} \int _0^TK(t)\textrm{e}^{-i\omega t} \textrm{d}t\right |^2 \end{align}

is the Fourier transform of the kinetic energy $K$ over a long time $T=50\,000$ . While $W=15$ clearly shows discrete incommensurate frequencies, $W=10$ shows a low-level broad band of frequencies. We still describe the latter as quasi-periodic, however, as it contains a small number of dominant frequencies that are incommensurate. Figure 14 used $\varepsilon =10^{-3}$ , and was run with a resolution of $(N_x, N_y) = (196, 128)$ . The results were robust to a doubling of the resolution in both directions. In addition, we checked that $\varepsilon =10^{-4}$ ran at $(N_x, N_y) = (384, 384)$ produced qualitatively similar results, despite the arrowheads being thinner (see Appendix C). The bursting scenario is still identifiable.

Finally, it is worth remarking that the transition scenario identified here is not the only one possible in vKf. In the regime discussed with $n=2$ and $L_x=6\pi$ , a bifurcation triggers irregular bursting events that become more frequent with increasing $W$ . Elastic turbulence is then seen above a critical $W$ where the disordered state exists at all times. An alternative transition to turbulence has been recently identified when $n=1$ and $L_x=8\pi$ , in which a sequence of period-doubling bifurcations cascade into turbulence (Nichols, Guy & Thomases Reference Nichols, Guy and Thomases2024).

Figure 16.

$(a, b)$ Simulations when $W=15$ (periodic, blue), $W=20$ (low-dimensional chaos, orange) and $W=30$ (turbulent, green), with rows of $(a)$ showing the trace field evolving over time, with times marked on the time series of $K$ in $(b)$ . $(c)$ The compensated power spectra for $W=15$ (blue), $20$ (orange), $30$ (green), $40$ (red), $50$ (purple), suggesting a regime where $E_K \sim k^{-4}$ , with regular spectra shown in the inset. $(d)$ Contributions to the kinetic energy due to the base shear ( $\mathcal {P}$ ), elastic forces ( $\mathcal {E}_{elast}$ ) and viscous forces ( $\mathcal {E}_{visc}$ ) as $W$ varies. Each quantity is averaged over a long time ( $T=20\,000$ ), with error bars showing plus and minus one standard deviation. Other parameters are $\beta =0.95$ , $Re=0.5$ , $\varepsilon =10^{-3}$ , $n=2$ and $L_x=8\pi$ in all plots.

5.2. Elastic turbulence properties

An analysis of ET in vKf is now presented in light of the fact that it occurs due to the centre-mode arrowhead losing stability as $W$ increases. We consider the power spectrum, as well as how energy is produced and dissipated in ET in a domain of $L_x=8\pi$ where turbulence exists across a wider range of $W$ (as per figure 13). These are compared with those of the centre-mode arrowhead structure. Three final states are considered: a periodic state at $W=15$ , a weakly chaotic state at $W=20$ and a strongly chaotic at $W=30$ (see figure 16). The polymer stress trace field is plotted at various times in figure 16(a), and the kinetic energy time series $K$ normalised by the laminar kinetic energy $K_0$ in figure 16(b). The periodic state at $W=15$ consists of two arrowheads drifting in the same direction (though still moving relative to each other). The low-dimensional chaos at $W=20$ again has two arrowheads drifting in the same direction, never making contact with each other. However, the increase in $W$ allows the ‘tail’ of the arrowhead to lengthen to the extent it now interacts with its ‘head’. We verify that the dynamics is truly chaotic, rather than quasi-periodic using the classification process discussed in Appendix B. The chaotic state at $W=30$ is more complicated, but we note that two processes can be identified, and we show these in the lower two rows of figure 16(a). The third row shows how multiple arrowheads can join together to form a zonal shear that spans the entire streamwise direction. This is a bursting event in which the kinetic energy peaks. The zonal shear then collapses, splitting into multiple arrowheads. The fourth row shows how the collision, coalescence and subsequent splitting of arrowheads contribute to the turbulent flow.

In figure 16(c) the compensated power spectrum of ET in vKf is plotted for various  $W$ which is consistent with $E_K\sim k^{-4}$ . This is the scaling found by Foggi Rota et al. (Reference Foggi Rota, Amor, Le Clainche and Rosti2024), Lellep et al. (Reference Lellep, Linkmann and Morozov2024) and Singh et al. (Reference Singh, Perlekar, Mitra and Rosti2024), and is close to that of Berti & Boffetta (Reference Berti and Boffetta2010) who found $E_K\sim k^{-3.8}$ . Interestingly the periodic state at $W=15$ , the weakly chaotic state at $W=20$ and the strongly chaotic states at $W=30$ all show similar scaling laws, suggesting that the centre-mode arrowhead contains the same small-scale structures present in ET. The region in which this scaling is valid appears to be dependent on $W$ , with the extent of this region falling with $W$ . Each of the simulations was run at a resolution of $(N_x, N_y) = (256, 256)$ , except for $W=50$ where $(N_x, N_y) = (512, 512)$ was used. This high-resolution run gave a power spectrum that was visually similar to that run at the lower resolution (not shown).

We now consider the total energy balance and see how the inputted energy due to the base shear compares with the viscous and elastic dissipation. The perturbative kinetic energy budget is derived by expanding the momentum (2.1) about the base state and dotting this with the velocity perturbation $\mathbf {u}^{\prime}$ . This produces

(5.2) \begin{align} Re \left (\frac {\partial u^{\prime}_i}{\partial t} + u^{\prime}_j \partial _j U_i + U_j \partial _j u^{\prime}_i + u^{\prime}_j \partial _j u^{\prime}_i \right ) u^{\prime}_i = \left (- \partial _i p^{\prime} + (1-\beta ) \partial _j \tau ^{\prime}_{ji} + \beta \nabla ^2 u^{\prime}_i \right ) u_i^{\prime} ,\end{align}

where quantities with a prime denote real perturbations from the laminar state, and Einstein summation is used over repeated indices. Incompressibility and the product rule allow us to rewrite this as

(5.3) \begin{align} & \frac {1}{2}\frac {\partial (u^{\prime}_i u^{\prime}_i)}{\partial t} + u^{\prime}_i u^{\prime}_j \partial _j U_i + \frac {1}{2} \partial _j ( U_j u^{\prime}_i u^{\prime}_i) + \frac {1}{2} \partial _j (u^{\prime}_j u^{\prime}_i u^{\prime}_i) = - \frac {1}{Re}\partial _i (p^{\prime} u_i^{\prime}) \nonumber\\ & \quad + \frac {(1-\beta )}{Re} \partial _j (\tau ^{\prime}_{ji} u_i^{\prime}) - \frac {(1-\beta )}{Re} \tau ^{\prime}_{ji}\partial _j u_i^{\prime} + \frac {\beta }{Re} \partial _j (u_i^{\prime} \partial _j u^{\prime}_i) - \frac {\beta }{Re} \partial _j u_i^{\prime} \partial _j u^{\prime}_i . \end{align}

Taking a volume average of (5.3) produces an energy budget of

(5.4) \begin{align} \frac {\textrm{d}}{\textrm{d}t}(\text{KE}) = \mathcal {P} + \mathcal {E}_{elast} + \mathcal {E}_{visc}, \end{align}

where

(5.5) \begin{align} \text{KE} = \frac {1}{2} \langle \mathbf {u^{\prime}} \cdot \mathbf {u^{\prime}} \rangle _V, \qquad \mathcal {P} = -\langle \mathbf {\nabla } \mathbf { U} : \mathbf {u^{\prime}} \mathbf {u^{\prime}} \rangle _V ,\end{align}

(5.6) \begin{align} \mathcal {E}_{elast} = -\frac {1-\beta }{Re}\bigl \langle \boldsymbol {\tau ^{\prime}} :{\mathbf {\nabla } \mathbf {u^{\prime}}} \bigr \rangle _V, \qquad \mathcal {E}_{visc} = -\frac {\beta }{Re}\bigl \langle {\nabla \mathbf {u}^{\prime}} : \nabla \mathbf {u^{\prime}} \bigr \rangle _V ,\end{align}

with $\langle \,\cdot \, \rangle _V$ denoting a volume average. This energy budget demonstrates that the kinetic energy can be energised or dissipated via the base shear through $\mathcal {P}$ , elastic forces through $\mathcal {E}_{elast}$ or viscous forces through $\mathcal {E}_{visc}$ . We plot in figure 16(d) the long-time averages of each energising and dissipative term as $W$ varies, noting that integrating (5.4) over a long time $T$ gives $\langle \mathcal {P}\rangle _T + \langle \mathcal {E}_{elast}\rangle _T + \langle \mathcal {E}_{visc}\rangle _T = 0$ . We see that the only energising process is elasticity, that the base shear contributes little to the kinetic energy and that viscous forces dissipate energy. This is consistent with Buza et al. (Reference Buza, Page and Kerswell2022b ), which considers the energy budget in channel flow of the centre-mode eigenfunctions at higher $Re$ . We see that the energy budget is qualitatively similar for flows which are chaotic ( $W=30$ ) and for those which are periodic ( $W=15$ ), but that as $W$ increases, elasticity provides more energy to the flow which is dissipated via viscosity.