3.1. Reference state

To investigate the influence of elasticity on viscoelastic flow focusing, it is essential to establish a well-defined reference state derived from experimental observations. The liquid phase is a dilute polymer solution, consisting of polyethylene oxide (PEO) dissolved in a glycerol–water mixture. This solution is characterised by a density of $\rho _1 = 1.0186 \times 10^3\ \textrm{kg}\,\textrm{m}^{-3}$ , a dynamic viscosity of $\mu _1 = 1.16 \times 10^{-3}\ \textrm{Pa s}$ , a relaxation time of $\lambda = 4.72 \times 10^{-4}\ \textrm{s}$ and a centreline velocity of $U_0 = 2.4\ \textrm{m/s}$ . These fluid properties can be experimentally tuned by adjusting the molecular weight and concentration of PEO in the solvent. The surrounding gas is nitrogen, with a density of $\rho _2 = 1.29\ \textrm{kg}\,\textrm{m}^{-3}$ and a dynamic viscosity of $\mu _2 = 2 \times 10^{-5}\ \textrm{Pa s}$ . The surface tension at the liquid–gas interface is $\sigma = 61.16\ \textrm{mN}\,\textrm{m}$ and the radius of the liquid jet is $R_1 = 0.07\ \textrm{mm}$ . Based on these parameters, the characteristic dimensionless groups are calculated as $ Re = 150$ , $\textit{We} = 7$ , $El = 0.1$ , $Q = 0.0013$ , $N = 0.0172$ and $X = 0.9$ .

As for the basic velocity profile, $U_{\!s}$ is larger than unity and represents the interfacial velocity normalised by the jet centreline velocity, while $K$ characterises the slope of the velocity profile at the liquid–gas interface. In the present reference state, the values $U_{\!s} = 1.33$ and $K = 1.2$ are determined by fitting theoretical predictions to experimental phase boundaries (see § 3.5). These parameters define a hyperbolic-tangent base profile that represents the experimentally realised shear-layer structure.

Figure 3.

Saddle points in the complex $k$ plane for the axisymmetric and helical modes under the reference state. The intersecting spatial branches arise from the downstream-propagating branch $k^{+}$ and the upstream-propagating branch $k^{-}$ . ( $a$ ) Axisymmetric mode: the negative imaginary part of the saddle-point frequency, $\omega _{0i} = -0.409134 \lt 0$ , indicates that the instability is convective. ( $b$ ) Helical mode: the negative imaginary part of the saddle-point frequency, $\omega _{0i} = -0.430328 \lt 0$ , likewise indicates a convective instability.

To validate the suitability of this reference state for spatial linear instability analysis, it is necessary to confirm that the base flow is convectively unstable, since spatial instability analysis is only applicable to convectively unstable configurations. We therefore perform a spatio–temporal instability analysis for the present reference state.

The analysis is based on the dispersion relation

(3.1) \begin{equation} D(n,\omega ,k; Re,We,El,U_{\!s},K,Q,N,X)=0, \end{equation}

where both the axial wavenumber $k$ and the perturbation frequency $\omega$ are complex quantities, $k = k_r + \textrm{i}k_i$ and $\omega = \omega _r + \textrm{i}\omega _i$ . The instability type is determined by locating a saddle point $(k_0,\omega _0)$ satisfying the zero-group-velocity condition

(3.2) \begin{equation} \frac {\partial \omega }{\partial k}=0. \end{equation}

In the complex $k$ plane, the solutions of $D=0$ form two spatial branches that, as $\omega _i$ decreases, approach each other and eventually coalesce at a saddle point. According to the Briggs–Bers criterion (Briggs Reference Briggs1964; Bers Reference Bers1973; Huerre & Rossi Reference Huerre and Rossi1998), only collisions between a downstream-propagating branch $k^{+}(\omega )$ and an upstream-propagating branch $k^{-}(\omega )$ are physically relevant. The sign of the imaginary part of the saddle-point frequency then distinguishes the instability type: $\omega _{0i}\gt 0$ corresponds to absolute instability, whereas $\omega _{0i}\lt 0$ indicates convective instability.

For the present reference state, the saddle point for the axisymmetric mode ( $n=0$ ) is located at $k_0 = 5.229015 - 0.921874\,\textrm{i}$ with the corresponding frequency $\omega _0 = 2.379980 - 0.409134\,\textrm{i}$ , as shown in figure 3( $a$ ). The negative imaginary part of the saddle-point frequency, $\omega _{0i}\lt 0$ , indicates that the instability is convective: perturbations are amplified while being advected downstream, without giving rise to self-sustained upstream-propagating disturbances. For the helical mode ( $n=1$ ), the saddle point converges to $k_0 = 5.355895 - 0.918524\,\textrm{i}$ with the corresponding frequency $\omega _0 = 2.255092 - 0.430328\,\textrm{i}$ , as shown in figure 3( $b$ ). The negative value of $\omega _{0i}$ again confirms the convective instability. These results demonstrate that, under the present reference state, both the axisymmetric and helical modes are convectively unstable. The absence of absolute instability justifies the use of a spatial instability framework in the following analysis.

3.2. Effects of liquid elasticity

Elasticity is the most distinctive feature of viscoelastic fluids and its influence on the flow-focused jet is first examined. In viscoelastic co-flow systems, the degree of elasticity is governed by both the intrinsic relaxation properties of the polymer chains and their concentration relative to the Newtonian solvent. The relaxation process, quantified through the elasticity number $El$ , characterises the competition between elastic and viscous time scales. A larger $El$ indicates that elastic stresses relax more slowly, allowing elastic effects to dominate the flow dynamics. The second controlling factor is the solvent–polymer composition ratio, represented by $X$ . A decrease in $X$ corresponds to an increase in the polymer fraction, thereby enhancing elastic contributions while reducing viscous damping from the solvent. Together, these two parameters determine how strongly elasticity influences the jet instability. In the present work, we focus on dilute polymer solutions and isolate the role of elasticity by varying $El$ while holding $X = 0.9$ fixed.

Figure 4.

Variations of the growth rate $-k_i$ and the phase speed $c_{\!p}$ with frequency $\omega$ for different elasticity numbers $El$ , at fixed parameters $ Re = 150$ , $\textit{We} = 7$ , $U_{\!s} = 1.33$ , $K = 1.2$ , $Q = 0.0013$ , $N = 0.0172$ and $X = 0.9$ : ( $a,c$ ) axisymmetric mode, ( $b,d$ ) helical mode. For $El = 0.1$ , the maximum growth rate $(-k_i)_{\textit{max}}$ and the most unstable phase speed $(c_{\!p})_{\textit{max}}$ are marked in ( $a$ ) and ( $c$ ), respectively. For convenience, the maximum growth rate $(-k_i)_{\textit{max}}$ and the most unstable phase speed $(c_{\!p})_{\textit{max}}$ are denoted as $-k_{i,\textit{max}}$ and $c_{\!p,\textit{max}}$ , respectively.

Figure 4 shows the spatial growth rate $-k_i$ and the phase speed $c_{\!p} = \omega /k_r$ as functions of $\omega$ for the axisymmetric and helical modes, highlighting the dependence of the dispersion relation on $El$ . In spatial instability theory, the frequency $\omega$ is prescribed to be real; therefore, $c_{\!p}=\omega /k_r$ represents the propagation speed of a wave of single real frequency $\omega$ . The spatial growth rates $-k_i$ are obtained by solving the quadratic eigenvalue problem (2.21) formulated in § 2.3. Consequently, each point on the dispersion curves corresponds to an eigenvalue of this problem. For a given real frequency $\omega$ , the criterion for identifying the spatial growth rate $-k_i$ is explained in Appendix B. By repeating this identification over the entire range of $\omega$ , the dispersion curves in figures 4( $a$ ) and 4(c) are constructed. Meanwhile, the eigenspectrum presented in Appendix D displays the full distribution of eigenvalues in the complex $k$ plane and their variation with elasticity. Further discussion of the rationale for the selection criterion, as well as of the structure of the eigenspectrum, can be found in Appendices B and D. Instability arises only within a finite frequency band, outside which disturbances decay exponentially. For each mode, the spatial growth rate reaches a maximum $-k_{i,\textit{max}}$ at the most unstable frequency $\omega _{\textit{max}}$ , corresponding to a resonance between the disturbance time scale and the dominant instability mechanism. This most amplified disturbance is therefore expected to dictate the jet dynamics in the unstable regime. The frequency $\omega _{\textit{max}}$ is associated with the axial wavenumber $k_{r,\textit{max}}$ , which sets the disturbance structure, while the wavelength $\lambda _{\textit{max}} = 2\pi /k_{r,\textit{max}}$ characterises the droplet size. The corresponding phase speed $c_{\!p,\textit{max}} = \omega _{\textit{max}}/k_{r,\textit{max}}$ defines the convective speed of the most unstable disturbance, thereby setting the characteristic propagation rate of perturbations along the jet. For both the axisymmetric mode (figure 4 $a$ ) and the helical mode (figure 4 $b$ ), $-k_{i,\textit{max}}$ initially increases slightly with $El$ and remains higher than that of the Newtonian case. At larger $El$ (e.g. $El=0.3$ ), however, the growth rate decreases sharply, ultimately rendering the viscoelastic jet more stable than its Newtonian counterpart. This trend is consistent with the temporal instability analysis of Ding et al. (Reference Ding, Mu, Si and Jian2022), who showed that jet instability reflects a competition between unrelaxed elastic stresses, which destabilise the flow, and nonlinear viscoelastic effects, which stabilise it. In addition, the right cutoff frequency shows a weak dependence on $El$ .

As shown in figures 4( $c$ ) and 4(d), the phase speed $c_{\!p}$ exhibits a non-monotonic dependence on frequency $\omega$ . At low frequencies, $c_{\!p}$ decreases with increasing $\omega$ for the axisymmetric mode. With further increase in frequency, the system enters a transitional regime in which $c_{\!p}$ varies only weakly with $\omega$ . At higher frequencies, $c_{\!p}$ increases monotonically with $\omega$ . In contrast, although the phase speed $c_{\!p}$ of the helical mode also decreases at low frequencies and increases again at high frequencies, it does not display a pronounced transitional regime as observed in the axisymmetric mode. At low frequencies, increasing the elasticity number $El$ reduces the phase speed, bringing it closer to the jet centreline velocity. At higher frequencies, however, increasing elasticity leads to a higher phase speed, indicating that elasticity enhances the propagation of high-frequency disturbances.

Figure 5( $a$ ) presents the maximum growth rate $-k_{i,\textit{max}}$ for the axisymmetric and helical modes as functions of $El$ . Consistent with figure 4, elasticity exerts a dual influence: it first enhances instability at small $El$ (e.g. $El \leqslant 0.01$ ) but suppresses it at larger $El$ . The mode with the larger $-k_{i,\textit{max}}$ is regarded as predominant, and the transition is clearly revealed by comparing the two curves, with the predominant mode at each $El$ marked by purple numbers. When $El$ is below approximately $0.08$ , the axisymmetric mode dominates, whereas increasing elasticity gradually shifts the predominance toward the helical mode. The critical elasticity number $El_{c} (\approx 0.08)$ , at which the two modes achieve equal spatial growth rates, signifies the onset of this mode transition. This highlights the role of elasticity in sustaining the axisymmetric mode at small $El$ and in promoting helical disturbances as $El$ increases.

Figure 5( $b$ ) shows the corresponding axial wavenumber $k_{r,\textit{max}}$ and phase speed $c_{\!p,\textit{max}}$ . The variation of $k_{r,\textit{max}}$ mirrors that of figure 5( $a$ ), increasing initially and then decreasing as $El$ grows. In contrast, the phase speed exhibits a mode-dependent trend: for the axisymmetric mode, $c_{\!p,\textit{max}}$ decreases with increasing $El$ , while for the helical mode, it shows only a slight decline. In other words, increasing elasticity drives the phase speed closer to the jet centreline velocity (i.e. $c_{\!p,\textit{max}}\!\to \!1$ ), indicating a tendency toward a ‘centre mode’. This feature is analogous to the centre-mode instability reported in viscoelastic Poiseuille flow (Garg et al. Reference Gañán-Calvo and Riesco-Chueca2018). Notably, this inward shift of the phase speed is consistent with the evolution of the eigenfunctions shown in figure 9, which also exhibit a progressive concentration of disturbance energy toward the jet core. The above trends of the growth rate and phase speed with respect to $El$ for the axisymmetric and helical modes are summarised in table 2.

Table 2.

Summary of the effects of elasticity number ( $El$ ) and Weber number ( $\textit{We}$ ) on the maximum growth rate $-k_{i,\textit{max}}$ , phase speed $c_{\!p,\textit{max}}$ and predominant instability mechanism for the axisymmetric and helical modes. As discussed in § 3.5, the CPI corresponds to surface-tension-dominated behaviour, the KHI is driven primarily by gas-pressure fluctuations or interfacial shear, and the ESI arises from energy transfer between the base-flow shear and the perturbation velocities.

Figure 5.

Variations of the maximum growth rate $-k_{i,\textit{max}}$ , the most unstable axial wavenumber $k_{r,\textit{max}}$ and the most unstable phase speed $c_{\!p,\textit{max}}$ with elasticity number $El$ , at fixed parameters $ Re = 150$ , $\textit{We} = 7$ , $U_{\!s} = 1.33$ , $K = 1.2$ , $Q = 0.0013$ , $N = 0.0172$ and $X = 0.9$ : ( $a$ ) the maximum growth rate $-k_{i,\textit{max}}$ , ( $b$ ) the most unstable axial wavenumber $k_{r,\textit{max}}$ (left axis) and the corresponding phase speed $c_{\!p,\textit{max}}$ (right axis). The numbers 0 and 1 indicate that the predominant mode under the current parameter set is the axisymmetric mode and the helical mode, respectively.

3.3. Effects of Weber number

To identify the parameter regimes in which the axisymmetric and helical modes dominate, we examine the influence of the Weber number ( $\textit{We}$ ) on the key instability characteristics: the maximum spatial growth rate $-k_{i,\textit{max}}$ , the most unstable axial wavenumber $k_{r,\textit{max}}$ and the corresponding phase speed $c_{\!p,\textit{max}}$ . Figure 6( $a$ ) shows the variation of $-k_{i,\textit{max}}$ for the axisymmetric and helical modes as $\textit{We}$ increases. As in figure 5, the predominant mode at each Weber number is indicated by the purple numbers, with $0$ and $1$ denoting the axisymmetric and helical modes, respectively. For the axisymmetric mode, $-k_{i,\textit{max}}$ decreases at low to moderate $\textit{We}$ (i.e. $\textit{We} \lt 6$ ), where the most unstable wavenumber satisfies $k_{r,\textit{max}}\!\lt \!1$ , indicating a long-wave regime (figure 6 $b$ ) in which surface tension promotes instability. As $\textit{We}$ increases further, the most unstable wavenumber shifts beyond unity ( $k_{r,\textit{max}} \! \gt \!1$ ) and $-k_{i,\textit{max}}$ rises again, marking a transition to the short-wave regime (figure 6 $b$ ), where surface tension instead exerts a stabilising influence. By contrast, the helical mode displays a monotonic increase of $-k_{i,\textit{max}}$ with $\textit{We}$ . This behaviour defines a critical Weber number $\textit{We}_c$ : for $\textit{We}\!\lt \!We_c$ , the axisymmetric mode dominates, whereas for $\textit{We}\!\gt \!We_c$ , the helical mode becomes predominant. The underlying predominant mode transition arises from the interplay among surface tension, interfacial pressure fluctuations, interfacial shear and energy transfer between disturbances and the base flow, and will be further clarified in § 3.4. Notably, this trend is qualitatively consistent with theoretical predictions and experimental observations of Newtonian flow focusing (Gañán-Calvo & Barrero Reference Gañán-Calvo and Riesco-Chueca1999; Rosell-Llompart & Gañán-Calvo Reference Rosell-Llompart and Gañán-Calvo2008; Si et al. Reference Si, Li, Yin and Yin2009, Reference Si, Li, Yin and Yin2010).

Figure 6.

Variations of the maximum growth rate $-k_{i,\textit{max}}$ , the maximum axial wavenumber $k_{r,\textit{max}}$ and the most unstable phase speed $c_{\!p,\textit{max}}$ with Weber number $\textit{We}$ , where $El = 0.01$ , $ Re = 150$ , $U_{\!s} = 1.33$ , $K = 1.2$ , $Q = 0.0013$ , $N = 0.0172$ and $X = 0.9$ : ( $a$ ) the maximum growth rate $-k_{i,\textit{max}}$ , ( $b$ ) the most unstable axial wavenumber $k_{r,\textit{max}}$ (left axis) and the corresponding phase speed $c_{\!p,\textit{max}}$ (right axis). The numbers 0 and 1 indicate that the predominant mode under the current parameter set is the axisymmetric mode and the helical mode, respectively.

Figure 6( $b$ ) presents the most unstable axial wavenumber $k_{r,\textit{max}}$ and the corresponding phase speed $c_{\!p,\textit{max}}$ for the axisymmetric and helical modes. The axial wavenumber $k_{r,\textit{max}}$ increases monotonically with $\textit{We}$ , indicating that higher $\textit{We}$ shifts the instability toward shorter wavelengths. This trend reflects the competition between inertia and capillarity: at large $\textit{We}$ , inertia dominates over surface tension, making the system more susceptible to short-wave disturbances. The phase speed, however, shows mode-dependent behaviour. For the axisymmetric mode, $c_{\!p,\textit{max}}$ increases with $\textit{We}$ , whereas for the helical mode, it decreases initially and then levels off. At sufficiently high $\textit{We}$ , the phase speeds of both modes converge toward the same limiting value, close to the interfacial velocity $U_{\!s}$ . This convergence highlights a key feature: under strong inertial conditions, the propagation of disturbances – irrespective of symmetry – is primarily governed by interfacial dynamics rather than by the bulk liquid or gas alone, consistent with the interpretation from the energy budget analysis (see figure 12), and can therefore be identified as an ‘interface mode’. The above features of the growth rate and phase speed with respect to $\textit{We}$ for the axisymmetric and helical modes are summarised in table 2.

3.4. Mechanisms of predominant mode transition

To elucidate the physical mechanisms responsible for the transition of the predominant modes, energy budget analysis was conducted. Figure 7 displays the energy budgets of the axisymmetric mode and the helical mode obtained from the spatial linear instability analysis. The contributions from PRL, NVL, SHL, TEL, NEL and NVG are found to be small in magnitude, each accounting for less than 5 % of the kinetic energy, KE, over the entire parameter range considered. Therefore, their influence on the jet instability is limited and they are omitted from the figure for clarity.

Figure 7.

Energy budget of different azimuthal modes under the reference state at $El = 0.1$ , $ Re = 150$ , $\textit{We} = 7$ , $U_{\!s} = 1.33$ , $K = 1.2$ , $Q = 0.0013$ , $N = 0.0172$ and $X = 0.9$ : ( $a$ ) axisymmetric mode, ( $b$ ) helical mode. The dashed lines denote the relative rate of change of energy at the most unstable frequency $\omega _{\textit{max}}$ . The definitions of the individual energy terms are given in (2.23). The physical interpretation of energy terms are given in table 1.

For the axisymmetric mode (figure 7 $a$ ), the surface-tension term (SUT) is positive in the low-frequency (long-wave) regime, showing that surface tension promotes instability, whereas at higher frequencies (short waves) SUT turns negative, indicating stabilisation. In contrast, for the helical mode (figure 7 $b$ ), SUT remains negative, i.e. consistently stabilising. The pressure-related term (PRG) is positive, indicating that interfacial pressure fluctuations enhance instability, while the dissipation term (DIS) is negative, reflecting the stabilising effect of viscous damping. The shear-related terms (SHG and ATG) are positive, confirming that shear drives destabilisation. The Reynolds-stress term (REY) is weak and slightly stabilising for the axisymmetric mode, but destabilising for the helical mode, at the most unstable frequency $\omega _{\textit{max}}$ . This contrast originates from different momentum-exchange regions: the axisymmetric mode dissipates energy near the interface, whereas the helical mode extracts energy from the liquid shear layer (see figure 10( $c$ )), as further discussed below. Regarding elasticity, the unrelaxed axial tension (NE) is positive and destabilising, whereas the polymeric elastic stress (DIP) is negative and stabilising, indicating that elasticity exerts competing effects on instability.

Figure 8 present the energy budget at the most unstable frequency $\omega _{\textit{max}}$ for varying $El$ and $\textit{We}$ . Since the predominant mode dictates the instability dynamics of the liquid jet, the energy budget is evaluated only for this mode, highlighted by the red numbers in the figure. Note that ATG and SHG both represent interfacial shear contributions; therefore, they are combined and shown as SHG. In figure 8( $a$ ), when elasticity is weak ( $El\! \lt \! 0.03$ ), the pressure-related term (PRG) provides the largest positive input, indicating that gas-pressure fluctuations are the predominant driver of the axisymmetric instability. Interfacial shear (SHG) also contributes positively, reinforcing the growth. The instability in this regime is therefore governed jointly by gas-pressure fluctuations and interfacial shear, consistent with a Kelvin–Helmholtz (KHI) mechanism. As $El$ increases to $0.1$ , however, the predominant mode shifts from the axisymmetric to the helical mode. During this transition, the elastic contributions (NE and DIP) both increase in absolute value but remain much smaller than the other energy terms, suggesting that elasticity does not directly dominate the energy exchange. Instead, elasticity modifies other energy pathways – most notably by enhancing Reynolds-stress transfer (REY). At sufficiently large $El$ ( $=\! 0.3$ ), REY may even exceed PRG, indicating that instability is then primarily driven by energy exchange between the liquid shear layer and perturbation velocities, representing a new predominant instability mechanism – elasticity-enhanced shear-driven instability (ESI) distinct from that in Newtonian jets. A further discussion of how elasticity governs this mode transition will be presented in a subsequent section.

Figure 8.

Energy budget of the predominant mode at the most unstable frequency: ( $a$ ) varying elasticity number $El$ , where $\textit{We} = 7$ ; ( $b$ ) varying Weber number $\textit{We}$ , where $El = 0.01$ . All other parameters are the same as in the reference case. Red numbers indicate the predominant mode.

In figure 8( $b$ ) the elasticity number is fixed at $El = 0.01$ , corresponding to weak elastic effects and behaviour close to that of a Newtonian fluid, thereby allowing the influence of $\textit{We}$ to be examined independently. For a very small Weber number, $\textit{We} = 2$ , SUT attains the largest positive value among all terms, indicating that surface-tension work dominates the instability and the mechanism is CPI, corresponding to a varicose (axisymmetric) mode. As $\textit{We}$ increases, both PRG and SHG gradually surpass SUT and become the leading positive contributions, signifying a transition to a KHI driven by gas-pressure fluctuations (PRG) and interfacial shear (SHG). As $\textit{We}$ increases toward a critical value around $6$ , SUT weakens and decreases, explaining the initial drop in the maximum spatial growth rate reported in § 3.3. With further increase in $\textit{We}$ , SUT changes sign and becomes negative, i.e. stabilising, while SHG and PRG intensify due to stronger interfacial shear and gas–inertia coupling, leading to renewed growth of the axisymmetric mode. At higher $\textit{We}$ (e.g. $\textit{We} = 9$ ), REY also rises as the predominant mode shifts to the helical mode, but this increase is mainly a by-product of the predominant mode transition rather than the emergence of a new predominant mechanism. In fact, the dominant source of instability remains aerodynamic, with PRG consistently maintaining the largest magnitude. Moreover, as $\textit{We}$ increases from $9$ to $12$ , the relative contribution of REY diminishes, while PRG and SHG continue to grow, underscoring that the instability is still governed primarily by aerodynamic effects. This redistribution of energetic pathways coincides with the transition of the predominant mode to the helical mode. A concise summary of the predominant instability mechanisms for the axisymmetric and helical modes identified across the parameter space is provided in table 2.

To clarify how elasticity influences the instability mechanisms, we complement the energy budget analysis with a structural examination of the disturbance eigenfunctions. While the budget terms indicate where growth is generated or dissipated, they do not fully account for the predominant mode transition on their own. We therefore inspect the eigenfunctions at the most unstable frequency $\omega _{\textit{max}}$ across varying $El$ , in order to expose the redistribution of amplitudes and the associated structural changes.

Figure 9.

Spatial structure of the disturbance eigenfunctions of the predominant mode at the most unstable frequency $\omega _{\textit{max}}$ for different elasticity numbers $El$ . Panel ( $a$ ) shows the axial disturbance velocity $|\hat {u}_1|$ exhibiting its radial distribution and streamwise evolution in the $(r,z)$ plane, while panel ( $b$ ) presents the corresponding radial disturbance velocity $|\hat {v}_1|$ with its radial and streamwise variations. All other parameters are identical to those of the reference state. Here, $|\boldsymbol{\cdot }|$ denotes the modulus.

Figure 9 illustrates the spatial structure of the disturbance eigenfunctions of the predominant mode at the most unstable frequency $\omega _{\textit{max}}$ , highlighting their radial distribution and streamwise evolution. As the disturbances evolve downstream along the $z$ direction, their amplitudes increase progressively, indicating spatial growth of the instability. When the elasticity number is small ( $El=0.01$ ), the amplitude of the predominant mode at $\omega _{\textit{max}}$ attains its maximum at the interface $r=1$ and is mainly concentrated in the vicinity of the interface, while its influence is relatively weak away from the interface. However, as the elasticity number increases to a larger value (e.g. $El=0.4$ ), the structure of the predominant mode at $\omega _{\textit{max}}$ gradually shifts from being interface localised to extending further into the jet core, reflecting a redistribution of disturbance energy from the interfacial region towards the interior of the jet. This trend parallels the observation in § 3.2 that the most unstable phase speed $c_{\!p,\textit{max}}$ approaches the centreline velocity as $El$ increases, indicating a gradual transition from an interface-dominated mode to a centre-mode structure.

Figure 10.

Disturbance eigenfunctions of the predominant mode at the most unstable frequency $\omega _{\textit{max}}$ for different elasticity numbers $El$ . All other parameters are identical to those of the reference state. The predominant mode is axisymmetric at small elasticity ( $El = 0.01$ and $0.07$ ), whereas it becomes helical at larger elasticity ( $El = 0.1$ and $0.4$ ). ( $a$ ) Axial velocity $\hat {u}_1$ for the predominant mode; ( $b$ ) radial velocity $\hat {v}_1$ for the predominant mode; ( $c$ ) real part of the cross-term $-\hat {u}_1^{*} \boldsymbol{\cdot }\hat {v}_1$ for the predominant mode at the most unstable frequency $\omega _{\textit{max}}$ for different elasticity numbers $El$ . Here, $(\boldsymbol{\cdot })^{*}$ denotes complex conjugation and $\textrm{Re}(\boldsymbol{\cdot })$ denotes the real part.

The disturbance eigenfunctions of the predominant mode display qualitatively similar radial structures at different axial positions $z$ . Therefore, in order to concentrate on the radial characteristics of the mode, we fix the axial coordinate at $z=0$ without loss of generality and present the corresponding radial profiles of the eigenfunctions at this location. Figures 10( $a$ ) and 10(b) show the axial and radial disturbance velocities of the predominant mode at the most unstable frequency $\omega _{\textit{max}}$ for different elasticity numbers $El$ . The predominant mode is axisymmetric for small elasticity numbers ( $El = 0.01$ and $0.07$ ), whereas it becomes helical at larger elasticity ( $El = 0.1$ and $0.4$ ). For the predominant axisymmetric mode at small elasticity ( $El = 0.01$ and $0.07$ ), the axial disturbance velocity $\hat {u}_1$ (see figure 10 a) is nearly uniform in the core and attains its maximum at the interface, characteristic of an interfacial mode. As elasticity increases and the predominant mode becomes helical ( $El = 0.1$ and $0.4$ ), the axial disturbance velocity develops a distinct peak in the jet interior, whose amplitude increases as $El$ is further raised. A similar restructuring occurs in the radial disturbance velocity $\hat {v}_1$ (see figure 10 b). For the predominant axisymmetric mode at small $El$ , $\hat {v}_1$ resembles the Newtonian limit and increases monotonically with radius. As elasticity increases and the predominant mode switches to the helical one, $\hat {v}_1$ develops a central maximum whose amplitude can even exceed the interfacial value, highlighting the strong core localisation of the disturbance.

Among the energy budget terms, the Reynolds-stress contribution (REY) exhibits the most pronounced variation across the predominant mode transition. Figure 10 $(c)$ shows the cross-correlation $-\textrm{Re}(\hat {u}_1^{*} \boldsymbol{\cdot }\hat {v}_1 )$ of the predominant mode as a function of $El$ , which directly links to the Reynolds-stress exchange identified in the energy budget. For $El=0.01$ and $El=0.07$ , the predominant mode is axisymmetric, whereas at larger $El$ it shifts to the helical mode. Here, $(\boldsymbol{\cdot })^{*}$ denotes complex conjugation and $\textrm{Re}(\boldsymbol{\cdot })$ denotes the real part. In flow-focusing configurations the liquid-phase velocity gradient $\textrm{d}U_1/\textrm{d}r$ is positive, and thus, the sign of REY follows directly from $-\textrm{Re}(\hat {u}_1^{*} \boldsymbol{\cdot }\hat {v}_1 )$ . At low elasticity (e.g. $El \le 0.07$ ), the axisymmetric mode exhibits negative values near the interface, with only a limited region of positive values inside the liquid, indicating that REY acts as an energy sink; its growth therefore relies primarily on interfacial mechanisms such as interfacial shear (SHG) and gas-pressure fluctuations (PRG). As $El$ increases to $0.1$ , however, $-\textrm{Re}(\hat {u}_1^{*} \boldsymbol{\cdot }\hat {v}_1 )$ becomes positive within the liquid shear layer, enabling the helical mode to efficiently extract energy from the velocity gradient. Quantitatively, the distribution shifts markedly toward the jet interior as $El$ increases from $0.07$ to $0.1$ , confirming that the predominant region of momentum exchange migrates from the interface to the liquid core. This indicates that the liquid shear layer (REY) becomes the primary source of instability.

3.5. Phase diagram of predominant mode transition

To elucidate how the predominant instability mode varies with $El$ and $\textit{We}$ , we construct a transition map in the $\textit{We}$ – $El$ plane at the reference parameters and compare theoretical predictions with experimental results. Figure 11 compares the predicted boundaries of predominant mode transition obtained from the spatial and temporal instability analyses with the experimentally observed mode distribution. A detailed description of the experimental set-up, measurement procedures and parameter calibration are presented separately in § 3.6, while the focus here is on the quantitative comparison between theory and experiment. The spatial instability analysis reveals a clear partitioning of modes in the $\textit{We}$ – $El$ space: at relatively large $\textit{We}$ , the helical mode dominates, consistent with figure 6, whereas at smaller $\textit{We}$ the predominant mode depends on $El$ – the axisymmetric mode dominates at low $El$ , while increasing $El$ gradually shifts the dominance to the helical mode, in line with figure 5.

Figure 11.

Phase diagram of the predominant modes in the $\textit{We}$ – $El$ plane at the reference parameters ( $ Re = 150$ , $U_{\!s}=1.33$ , $K=1.2$ , $Q=0.0013$ , $N=0.0172$ and $X=0.9$ ). The diagram compares the predicted boundaries of predominant mode transition obtained from the spatial (solid line) and temporal (dashed line) instability analyses with the experimentally observed mode distribution (symbols). The experimental cases are labelled in the order of increasing elasticity number as (I)–(VI), and these labels are indicated next to the corresponding data points in the present figure. The same labels (I)–(VI) are also used in figure 14 to denote the corresponding experimental images, so that each experimental condition in figure 14 can be directly located in the $\textit{We}$ – $El$ plane here. Red open squares denote axisymmetric modes, corresponding to cases (I), (II) and (IV), whereas blue open circles denote helical modes, corresponding to cases (III), (V) and (VI).

In contrast, the temporal instability analysis exhibits notable deficiencies: at small $El$ , it fails to capture the correct transition boundary, instead producing an approximately inverse relation in the $\textit{We}$ – $El$ plane, manifested as a sustained axisymmetric dominance over a broad $\textit{We}$ range. This behaviour stands in clear contradiction to spatial instability analysis and experimental observations, which reveal the emergence of the helical mode once $\textit{We}$ exceeds a finite threshold $\textit{We}_c$ . Consequently, the boundary obtained from spatial growth-rate analysis agrees much more closely with the experiments. More quantitatively, when $El = 0.2$ , the temporal analysis predicts a critical Weber number that is approximately $33\,\%$ larger than that obtained from the spatial analysis, and at $\textit{We} = 7$ the critical elasticity number from the temporal analysis exceeds the spatial prediction by more than $200\,\%$ . For completeness, a detailed examination of the temporal framework is presented in Appendix E. These discrepancies are consistent with the known limitations of the temporal (real- $k$ , complex- $\omega$ ) formulation in strongly convective, non-parallel jet configurations such as flow focusing: temporal instability analysis emphasises local amplification in time at a fixed station, while the instability observed in flow focusing is predominantly convective and manifests as downstream spatial growth. Consequently, the spatial (real- $\omega$ , complex- $k$ ) analysis provides a more faithful prediction of predominant mode boundaries and agrees more closely with experiments. Physically, the transition map highlights the interplay of inertia, capillarity and elasticity, while clarifying the limitations of temporal instability theory in jet instabilities.

Figure 12.

Phase diagram of primary instability mechanisms in the $\textit{We}$ – $El$ plane, with $ Re = 150$ , $U_{\!s} = 1.33$ , $K = 1.2$ , $Q = 0.0013$ , $N = 0.0172$ and $X = 0.9$ . Open symbols denote the axisymmetric mode ( $n=0$ ) and solid symbols denote the helical mode ( $n=1$ ).

Since the transition of the predominant mode reflects a competition among distinct instability mechanisms, we further perform an energy budget analysis of the predominant mode. Figure 12 presents the phase diagram of the primary instability mechanisms in the $\textit{We}$ – $El$ plane, corresponding to the operating conditions shown in figure 11. Open symbols denote the axisymmetric mode, while solid symbols denote the helical mode. These mechanisms are identified according to the dominant energy transfer term in the budget analysis. Specifically, SUT characterises the instability driven by surface tension, corresponding to the classical CPI. Here PRG and SHG are associated with the KHI, reflecting the destabilising effects of gas-pressure fluctuations and velocity shear at the liquid–gas interface, respectively. In contrast, REY denotes an ESI, arising from the energy exchange between the disturbance and the base flow through Reynolds stress, typically linked to shear-driven amplification within the liquid core. From a physical perspective, at low $\textit{We}$ ( $\leqslant 2$ ), surface tension governs the jet dynamics, and disturbances are amplified through the interplay between interface curvature and capillary forces, such that SUT dominates. As $\textit{We}$ increases, the influence of surface tension weakens while gas-pressure fluctuations and interfacial shear become more pronounced, leading to a transition toward PRG and SHG, characteristic of KHI. With a further increase in $\textit{We}$ , viscoelastic jets with small elasticity behave similarly to Newtonian jets, where KHI remains the dominant mechanism. At higher elasticity (e.g. $El = 0.3$ ), however, REY may emerge as the prevailing mechanism, indicating that instability is mainly driven by shear amplification inside the liquid, thus revealing a new predominant instability mechanism that differs fundamentally from that in Newtonian jets. Eventually, at sufficiently large $\textit{We}$ (e.g. $\textit{We}=10$ ), interfacial dynamics regain dominance and the instability mechanism shifts back to KHI.

3.6. Experimental set-up

This section supplements the experimental results presented in § 3.5 by providing a detailed description of the measurement procedures, flow-focusing set-up and data interpretation. Rather than repeating the results, it aims to specify the experimental conditions under which the data in figure 11 were obtained. The following describes the flow-focusing configuration and measurement procedures used to obtain these data. It should be noted that experimental observations show that, in the region where the instability emerges and develops, the jet diameter varies only weakly along the axial direction and no noticeable tilt of the jet core is observed. Therefore, the assumption of a constant-radius, axisymmetric base flow adopted in the theoretical analysis is justified in this region.

To present quantitative experimental validations with the theoretical results, a flow-focusing device is used to generate a coaxially gas-driven viscoelastic jet. In a flow-focusing process the liquid jet is driven by a co-flowing high-speed gas stream, presenting either the axisymmetric or the helical mode of instability as it evolves downstream. The experimental platform is constituted by a capillary tube (with inner radius $R$ and outer radius $R_0$ ) and a sealed chamber with an orifice (with diameter $D$ and thickness $T$ ) at the bottom, coaxial with the capillary tube, as shown by the sketch map of figure 13( $a$ ) (Wang et al. Reference Wang, Mu, Zhao, Wu, Xu, He and Si2024). The distance between the exit of the tube and the focusing orifice is denoted by $H$ . In our experiments, the geometric parameters are maintained constant with $R=525\, \unicode {x03BC}$ m, $R_0=630\, \unicode {x03BC}$ m, $H=700\, \unicode {x03BC}$ m, $D=800\, \unicode {x03BC}$ m and $T=700\, \unicode {x03BC}$ m, respectively. The sealed chamber is assembled by transparent glass windows on all directions for the convenience of observation. A syringe pump is utilised to supply the viscoelastic focused liquid (with density $\rho _1$ , viscosity $\mu _1$ , flow rate $Q_1$ , surface-tension coefficient $\sigma$ and relaxation time $\lambda$ ) into the capillary tube and compressed nitrogen gas (with density $\rho _2 = 1.29\ \textrm{kg}\,\textrm{m}^{-3}$ , viscosity $\mu _2 = 1.81 \times 10^{-5}\ \textrm{Pa s}$ ) is supplied by a gas cylinder connected to a flow meter, serving as the focusing gas. One probe of the differential pressure transducer is inserted into the chamber and used to measure the gas-pressure difference $\Delta p$ . The viscoelastic focused liquids employed in the experiments are PEO solutions dissolved in glycerol–water mixtures (with solvent viscosity $\mu _s$ and molecular weight of polymer $M_w$ ). The mass fraction of PEO is fixed at $c = 0.05 \,\%$ for all solutions, while the mass fraction of glycerol, $c_G$ , varies among different samples. Their physical properties are listed in detail in table 3.

The dimensionless parameters and characteristic scales employed in the experiment are consistent with those presented in § 2 and are given by

(3.3) \begin{equation} Re = \frac {\rho _1 U_0 R_1}{\mu _1},\quad We = \frac {\rho _1 U_0^2 R_1}{\sigma },\quad El = \frac {\lambda \mu _1}{\rho _1 R_1^{2}}. \end{equation}

It is well known that changing the mass fraction $c$ or molecular weight of the polymer $M_w$ in PEO solutions leads to changes in both the solution viscosity $\mu _1$ and the relaxation time $\lambda$ . Although different solution formulations are attempted to achieve approximately identical $\mu _1$ , it remains difficult to obtain a sufficient number of data points while keeping $ Re$ constant. However, the following characteristics are present in flow focusing and have been experimentally verified (Gañán-Calvo Reference Gañán-Calvo, Montanero, Martín-Banderas and Flores-Mosquera1998; Herrada et al. Reference Herrada, Gañán-Calvo, Ojeda-Monge, Bluth and Riesco-Chueca2008; Si et al. Reference Si, Li, Yin and Yin2009).

Table 3.

Physical properties of the focused liquids.

Figure 13.

( $a$ ) Sketch of the flow-focusing system considering the viscoelastic focused phase and the focusing gas phase, where $Q_1$ represents the flow rates of the focused liquid and the gas flow rate is determined by the pressure difference $\Delta p$ , $R$ and $R_0$ denote the inner and outer radii of the capillary tube, $D$ and $T$ represent the diameter and thickness of the orifice, and $H$ denotes the distance between the capillary tube and the orifice, respectively. ( $b$ ) Jet of the $n = 0$ mode downstream of the orifice exit is shown in an enlarged view. Here $U_1$ and $U_{\!a}$ denote the jet and the focusing gas velocity, respectively, and $R_1$ represents the jet radius. Additionally, $U_0$ and $U_{\!a}$ denote the base-flow velocity at the jet centreline and gas flow velocity at $r=D/2$ , respectively, and $R_1$ represents the jet radius.

1. (i) The gas flow velocities upstream and downstream of the focused orifice are governed by the Bernoulli equation $\Delta p = \rho _2 U_{\!a}^2 /2$ , so the external focusing gas velocity $U_{\!a} = \sqrt {2 \Delta p /\rho _2}$ is obtained.

2. (ii) The gas-to-liquid momentum ratio is very close to unity in flow focusing, i.e. $(1/2) \rho _{1} U_{0}^{2}=(1/2) \rho _{2} U_{a}^{2}$ , which is determined by the simplified averaged momentum equation in the axial direction in which the viscous terms are neglected. Therefore, the jet velocity can be given by $U_0 = \sqrt {\rho _2 /\rho _1} \ U_{\!a}$ .

3. (iii) The axial velocity of the jet is approximately uniform, so that $Q_1 = \pi R_1^2 U_0$ .

From the above relationships, the jet radius is obtained as $R_{1}= (\rho _{1}/ ({2 \pi ^{2} \Delta p} ) )^{1 / 4} Q_{1}^{1 / 2}$ , which is in good quantitative agreement with experimental measurements (Gañán-Calvo Reference Gañán-Calvo, Montanero, Martín-Banderas and Flores-Mosquera1998; Si et al. Reference Si, Li, Yin and Yin2009). Accordingly, the dimensionless parameters $ Re$ , $\textit{We}$ and $El$ can be rewritten as

(3.4) \begin{align} Re &= \left ( \frac {2 \rho _{1}^{3}}{\pi ^{2} \mu _{1}^{4}} \right )^{1/4} \left (Q_{1}^{2} \Delta p \right )^{1/4}, \end{align}

(3.5) \begin{align} We &= \left ( \frac {8 \rho _{1}}{\pi ^{2} \sigma ^{4}} \right )^{1/4} \left (Q_{1}^{2} \Delta p \right )^{1/4} (\Delta p)^{1/2}, \end{align}

(3.6) \begin{align} El &= \left ( \frac {2 \pi ^{2} \mu _{1}^{2} \lambda ^{2}}{\rho _{1}^{3}} \right )^{1/2} \left (Q_{1}^{2} \Delta p \right )^{-1/2} \Delta p . \end{align}

Therefore, by controlling $ (Q_{1}^{2} \Delta p )$ constant and varying only $\Delta p$ in the experiments, $Re$ $( \approx 150)$ can be conveniently maintained unchanged, while both $\textit{We}$ $(4\sim 10)$ and $El$ $(10^{-3}\sim 1)$ can be simultaneously adjusted. Using this method, a set of data points in the $\textit{We}$ – $El$ plane can be obtained with a single viscoelastic solution, and these points form a sloped trajectory in figure 11.

Figure 14.

The transition of the jets from the $n = 0$ mode to the $n = 1$ mode is observed in the experiments with $ Re \approx 150$ . Scale bars of 400 μm are provided in the figure. ( $a$ ) For the polymer with a molecular weight of $M_w=100$ kg mol⁻¹, where elastic effects are relatively weak, the non-axisymmetric instability of the jet is primarily induced by an increase in the $\textit{We}$ number. ( $b$ ) For the polymer with a molecular weight of $M_w=600$ kg mol⁻¹, where elastic effects are stronger, the helical instability is mainly caused by an increase in the $El$ number. The corresponding experimental results are displayed in the $\textit{We}$ – $El$ plane using open symbols (see figure 11).

As predicted by theory, with increasing $\textit{We}$ and $El$ , the dynamic evolution of the jet transitioning from the axisymmetric mode to the helical mode is clearly demonstrated in figure 14. The corresponding experimental results are displayed in the $\textit{We}$ – $El$ plane using hollow symbols (see figure 11). When the solution exhibits weak elasticity, as shown in figure 14( $a$ ), the variation in $El$ achieved in the experiments is relatively small. With the increase in $\textit{We}$ , the dominant mode gradually transforms from the axisymmetric mode to the helical mode, which validates the theoretical analysis presented in figure 6. In contrast, when the elasticity of the solution becomes strong, as depicted in figure 14( $b$ ), the increase in $El$ has a more pronounced influence on the mode transition: even at comparable $\textit{We}$ values, a higher $El$ promotes the earlier emergence of the helical mode, which can be seen by comparing figure 14( $a$ ) at $\textit{We}=7.6$ with figure 14( $b$ ) at $\textit{We}=7.4$ . This observation is consistent with the trends shown in figure 5.