2. Mathematical description

We consider an incompressible axisymmetric slender viscoelastic jet with radius $r(z,t)$ and mean axial velocity $u(z,t)$ . The solvent has constant density $\rho$ and viscosity $\mu _s$ , whereas the polymer is characterised by its viscosity $\mu _p$ and relaxation time $\lambda$ . The liquid–air surface tension coefficient is $\sigma$ , assumed constant as well. The description of the jet dynamics combines the leading-order mass and momentum equations, derived originally by Eggers & Dupont (Reference Eggers and Dupont1994), augmented by the extra stresses due to the presence of the polymer (Ardekani et al. Reference Ardekani, Sharma and McKinley2010), given by

(2.1) \begin{equation} \partial _t\big(r^2\big) + \partial _z\big(u r^2\big)=0, \end{equation}

and

(2.2) \begin{equation} \rho \left ( \partial _t u + u\,\partial _z u \right ) = \rho g - \sigma \partial _z \mathcal{C} + \frac {3 \mu _s }{r^2}\partial _z\!\left (r^2\,\partial _z u\right ) + \frac {1}{r^2}\partial _z\big(r^2\,T\big), \end{equation}

where $g$ is the gravitational acceleration, $\mathcal{C}$ is the full curvature,

(2.3) \begin{equation} \mathcal{C} = \frac {1}{r}\left (1+r_z^2\right )^{-1/2} \!\!\! - r_{\textit{zz}}\left (1+r_z^2\right )^{-3/2}, \end{equation}

and $T = \tau _{\textit{zz}} - \tau _{\textit{rr}}$ is the polymeric tensile stress, written in terms of the diagonal components of the polymeric extra-stress tensor $\boldsymbol{\tau }$ . The evolution of $\boldsymbol{\tau }$ is assumed to be dictated by the Oldroyd-B/Giesekus constitutive model

(2.4) \begin{equation} \boldsymbol{\tau } + \lambda \,\overset {\triangledown }{\boldsymbol{\tau }} + \frac {\alpha \lambda }{\mu _p}\,\boldsymbol{\tau }\boldsymbol{\cdot }\boldsymbol{\tau } = 2\mu _p\,\boldsymbol{D}, \end{equation}

with $\boldsymbol{D}=1/2(\boldsymbol{\nabla }\boldsymbol{v}+\boldsymbol{\nabla }\boldsymbol{v}^T)$ , $\overset {\triangledown }{\boldsymbol{\tau }} = {\mathrm{D}}/{ \mathrm{D}t} \boldsymbol{\tau }- (\boldsymbol{\nabla }\boldsymbol{v})^T \boldsymbol{\cdot }\boldsymbol{\tau } - \boldsymbol{\tau } \boldsymbol{\cdot }(\boldsymbol{\nabla }\boldsymbol{v} )$ the upper-convected derivative, $ ({\mathrm{D}}/{ \mathrm{D}t})$ the material derivative operator and $\alpha$ is the Giesekus mobility factor. Note that setting $\alpha =0$ recovers the Oldroyd-B model, which will be the choice for the rest of this paper. The diagonal components of (2.4) read

(2.5) \begin{align} \lambda \big (\partial _t\tau _{\textit{zz}} + u\partial _z\tau _{\textit{zz}} - 2\,u_z\,\tau _{\textit{zz}}\big ) + \tau _{\textit{zz}} + \frac {\alpha \lambda }{\mu _p}\,\tau _{\textit{zz}}^2 &= 2\mu _p\,u_z, \end{align}

(2.6) \begin{align} \lambda \big (\partial _t\tau _{\textit{rr}} + u\,\partial _z\tau _{\textit{rr}} + u_z\,\tau _{\textit{rr}}\big ) + \tau _{\textit{rr}} + \frac {\alpha \lambda }{\mu _p}\,\tau _{\textit{rr}}^2 &= -\mu _p \,u_z. \\[10pt] \nonumber \end{align}

Note that, whereas finite extensible nonlinear elastic models are more appropriate when the polymer approaches its finite extension limit, since the present study concerns the weakly elastic, early linear stage of the instability about a steady stretched base state, the Giesekus/Oldroyd-B closure is deemed as a baseline constitutive model, while finite extensibility is more appropriate in scenarios with stronger stretching and/or nonlinear post-critical dynamics.

The mathematical description is completed upon boundary conditions, namely,

(2.7) \begin{equation} r(0,t)-R \;=\; u(0,t)-\frac {Q}{\pi R^2} \;=\; \tau _{\textit{zz}}(0,t) \;=\; \tau _{\textit{rr}}(0,t) \;=\; 0, \end{equation}

where $R$ is the nozzle radius and $Q$ is the prescribed flow rate. The assumption of fully relaxed polymer at the exit, $\tau _{\textit{zz}}(0)=\tau _{\textit{rr}}(0)=0$ , provides a minimal and internally consistent closure for a free-jet model posed only for $z\geqslant 0$ . It is appropriate when the residence time in the injector is long compared with the relaxation time $\lambda$ and when upstream deformation rates are weak. At the same time, it is not expected to be universally accurate in experiments, since polymer solutions may exit the injector with a non-negligible elastic memory of the internal shear and extension history, as evidenced by the classical die-swell (Barus) effect. In particular, extrudate swell in dilute polymer solutions is driven primarily by shear-generated normal stresses, notably the first normal-stress difference, rather than by a stress-free state at the die exit (e.g. Mitsoulis Reference Mitsoulis2010), and elasticity is known to influence jetting–dripping thresholds in related jetting configurations (Montanero & Gañán-Calvo Reference Montanero and Gañán-Calvo2008).

In the present work we nevertheless adopt (2.7) as the reference inlet-stress condition. The reason is that any non-zero specification of $\tau _{\textit{zz}}(0)$ and $\tau _{\textit{rr}}(0)$ within a reduced free-surface jet model remains an effective closure unless the internal nozzle flow and die-exit region are modelled explicitly. Asymptotic inlet-stress relations derived from the outer jet equations, as in local slender analyses of viscoelastic jets under gravity (Alhushaybari & Uddin Reference Alhushaybari and Uddin2019), provide mathematically consistent alternatives, but they still do not encode the dependence on injector geometry, entrance effects or residence-time history that controls die swell in practice. Likewise, prescribing finite inlet stresses directly introduces additional closure parameters whose values are not fixed by the present reduced model alone. For these reasons, the relaxed condition (2.7) should be interpreted here as a baseline choice that isolates the role of viscoelasticity in the downstream stretched-jet dynamics, rather than as a universally realistic experimental boundary condition. To assess the robustness of the main conclusions, we examine in Appendix C the sensitivity of the base state and global-stability metrics to be defined below to alternative non-zero inlet-stress closures. A fully coupled nozzle–jet description, or a calibrated inlet condition based on independent measurements of die swell and upstream rheology, remains a natural extension of the present framework and is left for future work.

Lastly, initial conditions are also needed for the temporal marching of the model. We shall describe them below if needed.

2.1. Choice of scales and non-dimensionalisation

Following Rubio-Rubio et al. (Reference Rubio-Rubio, Sevilla and Gordillo2013), we take the gravito–capillary scales

(2.8) \begin{equation} \ell _\sigma = \left (\frac {\sigma }{\rho g}\right )^{1/2}, \quad U_\sigma =\sqrt {g\ell _\sigma }, \quad T_\sigma =\sqrt {\ell _\sigma /g}, \quad \tau _c = \mu _p U_\sigma /\ell _\sigma , \end{equation}

as relevant scales for length, velocity, time and polymer stress, respectively. With these choices, the dimensionless model reads

(2.9) \begin{align}&\qquad\qquad\qquad\qquad \partial _t\big(r^2\big) + \partial _z\big(u r^2\big)=0, \end{align}

(2.10) \begin{align}& \partial _t u + u\,\partial _z u = 1 - \partial _z \mathcal{C} + \frac {\varGamma }{r^2}\partial _z\!\left (r^2\,\partial _z u\right ) + \frac {\varGamma (1-\beta )}{3\beta r^2}\partial _z\big(r^2\,T\big), \end{align}

(2.11) \begin{align}&\quad\,\, \textit{De}\big (\partial _t\tau _{\textit{zz}} + u\,\partial _z\tau _{\textit{zz}} - 2\,u_z\,\tau _{\textit{zz}}\big ) + \tau _{\textit{zz}} + {\alpha \textit{De}}\,\tau _{\textit{zz}}^2 = 2 u_z, \end{align}

(2.12) \begin{align}&\quad\,\, \textit{De}\big (\partial _t\tau _{\textit{rr}} + u\,\partial _z\tau _{\textit{rr}} + u_z\,\tau _{\textit{rr}}\big ) + \tau _{\textit{rr}} + {\alpha \textit{De}} \,\tau _{\textit{rr}}^2 = -\,u_z, \end{align}

where each dependent variable should be understood as scaled with (2.8), i.e. $z \to z/\ell _\sigma$ , $u \to u/U_\sigma$ , $t \to t/T_\sigma$ , $\tau _{ij} \to \tau _{ij}/\tau _c$ . In these units, the dimensionless counterpart of the baseline boundary conditions (2.7) is given by

(2.13) \begin{equation} r(0,t)={\textit{Bo}}^{1/2}, \quad u(0,t) = \textit{We}^{1/2}{\textit{Bo}}^{-1/4}, \quad \tau _{\textit{zz}}(0,t) = \tau _{\textit{rr}}(0,t) = 0. \end{equation}

Alternative non-zero inlet-stress closures are examined separately in Appendix C. Suitable outlet conditions must also be prescribed downstream in order to numerically integrate these equations. A possible choice is a jet with zero axial curvature, i.e. $r_z = r_{\textit{zz}} = 0$ , at $z = L \gg 1$ , as done for instance by Shukla & Gallaire (Reference Shukla and Gallaire2020). Above, ${\textit{Bo}} = \rho g R^2/\sigma$ is the Bond number, $\textit{We} = \rho U^2 R/\sigma$ is the Weber number based on the mean exit speed $U = Q/(\pi R^2)$ , $\varGamma = {3\nu _s}/{\sqrt {g\ell _\sigma ^3}}$ is the Kapitza number based on the solvent kinematic viscosity $\nu _s = \mu _s/\rho$ , $\textit{De}=\lambda /T_\sigma = \lambda \sqrt {g/\ell _\sigma }$ is the Deborah number and $\beta = \mu _s/(\mu _s + \mu _p)$ is the solvent-to-total viscosity ratio. Finally, note that the Newtonian limit is recovered for $1-\beta = \textit{De} = 0$ . It is worth noting that the stress scale $\tau _c=\mu _p U_\sigma /\ell _\sigma$ measures the polymer stress relative to its own polymer-viscous scale, whereas its feedback on the momentum equation is weighted by the prefactor $(1-\beta )/\beta$ . Thus, for $\beta \lesssim 1$ the dimensionless stress fields may remain $O(1)$ while their dynamical effect on the jet becomes asymptotically weak; in the exact limit $\beta =1$ , the polymer contribution vanishes and the appropriate description is the Newtonian system obtained by dropping the polymer-stress equations altogether.

3. The linearised equations

We assume the usual normal-mode decomposition for each fluid variable $\boldsymbol{q} = [r, u, \tau _{\textit{zz}}, \tau _{\textit{rr}}]^{\mathrm{T}}$ about a steady base

(3.1) \begin{equation} \boldsymbol{q}(z,t) = \boldsymbol{q}_0(z) + \epsilon \hat {\boldsymbol{q}}(z)\mathrm{e}^{\omega t}, \end{equation}

where $\epsilon \ll 1$ is an arbitrarily small number, $\boldsymbol{q}_0(z)$ is the steady solution to (2.10)–(2.12), $\omega =\omega _r +\mathrm{i}\omega _i$ is a complex growth rate whose real and imaginary parts represent the temporal amplification rate and oscillation frequency, respectively, and $\hat {\boldsymbol{q}}$ are (direct) eigenfunctions. Substituting the ansatz (3.1) into (2.10)–(2.12) and collecting $O(\epsilon )$ terms yields a linear system that can be written in compact form as the generalised eigenvalue problem

(3.2) \begin{equation} \left ( \mathcal{A} - \omega \mathcal{B} \right ) \hat {\boldsymbol{q}} = \boldsymbol{0}, \end{equation}

where $\mathcal{A}$ and $\mathcal{B}$ are linear differential operators encompassing boundary conditions. We employ a standard Chebyshev collocation technique in order to discretise (3.2) (see Appendix B for details) which in turn yields the matrix pencil $({{\unicode{x1D63C}}} - \omega {{\unicode{x1D63D}}}) \boldsymbol{\tilde {q}} = 0$ , where ${\unicode{x1D63C}}, {\unicode{x1D63D}}$ and $\tilde {\boldsymbol{q}}$ represent the discretised counterparts, whose spectrum $\omega$ determines the global temporal stability of the spatially varying base flow. This full-curvature formulation recovers the Newtonian stretched-jet framework of Rubio-Rubio et al. (Reference Rubio-Rubio, Sevilla and Gordillo2013) in the appropriate limit and here serves as the basis for a direct–adjoint analysis of weakly viscoelastic mode selection. Because the base state varies significantly with $z$ under gravitational acceleration and capillary thinning, local parallel stability theory is insufficient to capture the onset of self-excited jet oscillations. Thus, a global approach is the correct strategy to predict stability thresholds and oscillation frequencies in the jetting regime. In the discretised setting, we define a weighted inner product consistent with the spectral collocation quadrature

(3.3) \begin{equation} \langle \boldsymbol{a},\boldsymbol{b}\rangle _Q \equiv \boldsymbol{a}^{\mathrm{H}}\,{\unicode{x1D64C}}\,\boldsymbol{b}, \end{equation}

where $(\boldsymbol{\cdot })^{\mathrm{H}}$ denotes the conjugate transpose and $\unicode{x1D64C}$ is a Hermitian positive–definite weight matrix assembled from the collocation weights (see, for instance, Trefethen (Reference Trefethen2000)). Direct eigenfunctions are normalised such that

(3.4) \begin{equation} \langle \tilde {\boldsymbol{q}} , \tilde {\boldsymbol{q}} \rangle _Q = 1. \end{equation}

To quantify the sensitivity of the global modes to localised modifications of the base flow and/or of the linear operator, we also introduce the corresponding adjoint (left) eigenfunctions. From a mathematical standpoint, the adjoint eigenfunction provides the receptivity kernel of the global oscillator, i.e. it weights external forcing and initial perturbations in the energy inner product, thereby filtering generic disturbances into the component that most efficiently excites the selected direct mode and pinpointing the region(s) of strongest feedback sensitivity. The (discrete) adjoint eigenfunctions $\tilde {\boldsymbol{q}}^{\dagger }$ are then defined as solutions of the $Q$ -adjoint generalised eigenvalue problem

(3.5) \begin{equation} \big ({\unicode{x1D63C}}^{\dagger }-\bar {\omega }\,{\unicode{x1D63D}}^{\dagger }\big )\tilde {\boldsymbol{q}}^{\dagger }=\boldsymbol{0}, \quad {\unicode{x1D63C}}^{\dagger }={\unicode{x1D64C}}^{-1}{\unicode{x1D63C}}^{\mathrm{H}}{\unicode{x1D64C}},\quad {\unicode{x1D63D}}^{\dagger }={\unicode{x1D64C}}^{-1}{\unicode{x1D63D}}^{\mathrm{H}}{\unicode{x1D64C}}, \end{equation}

and we impose the standard bi-orthogonality (normalisation) condition

(3.6) \begin{equation} \big \langle \tilde {\boldsymbol{q}}^{\dagger },{\unicode{x1D63D}}\,\tilde {\boldsymbol{q}} \big \rangle _Q = \tilde {\boldsymbol{q}}^{\dagger \mathrm{H}}{\unicode{x1D64C}}\,{\unicode{x1D63D}}\,\tilde {\boldsymbol{q}} =1, \end{equation}

which enables a direct projection of perturbations onto individual global modes.

Building on this framework, we characterise the spatial structure of the global instability using two complementary analyses based on the sensitivity of eigenvalues with respect to localised changes of the operator $\mathcal{A}$ . First, the wavemaker (or structural sensitivity) field measures the sensitivity of $\omega$ to a localised feedback perturbation of the operator and is typically constructed from the pointwise overlap of the direct and adjoint eigenfunctions; it highlights the region where the instability mechanism is most receptive to modifications of the local coupling between state variables, and is therefore a proxy for the core of the global feedback loop (Giannetti & Luchini Reference Giannetti and Luchini2007). Equivalently, the wavemaker identifies the spatial region where cause and effect overlap most strongly. The direct mode is large enough there for the local dynamics to matter, and the adjoint mode is also large enough there for the eigenvalue to react strongly to a perturbation. For that reason, it is often interpreted as the core of the global oscillator, namely the part of the flow where a small structural change has the greatest impact on the selected growth rate and frequency. It is defined as the real scalar

(3.7) \begin{equation} S(z) = \left \lVert {\unicode{x1D64C}} \tilde {\boldsymbol{q}}^\dagger \right \rVert \, \left \lVert \tilde {\boldsymbol{q}}\right \rVert , \end{equation}

which corresponds to the Frobenius norm of the structural-sensitivity tensor

(3.8) \begin{equation} {\unicode{x1D64E}} = {\unicode{x1D64C}} \tilde {\boldsymbol{q}}^\dagger \otimes \tilde {\boldsymbol{q}}. \end{equation}

Additionally, the momentum-restricted wavemaker is defined as

(3.9) \begin{equation} S_{\textit{mom}}(z) = \big \lVert {\unicode{x1D650}} \big( {\unicode{x1D64C}} \tilde {\boldsymbol{q}}^\dagger \big) \big \rVert \, \left \lVert {\unicode{x1D650}} \tilde {\boldsymbol{q}}\right \rVert , \end{equation}

where $\unicode{x1D650}$ is a projection matrix that extracts the subset of state components associated with the momentum equation. Numerically, the wavemaker is assembled node by node on the collocation grid. For each axial point $z_j$ , one extracts the local direct and adjoint state vectors and computes the corresponding Euclidean norm in state space. Thus $S(z_j)$ is the Frobenius norm of the local structural-sensitivity tensor, whereas $S_{\textit{mom}}(z_j)$ is the same quantity restricted to the momentum component through the projector $\unicode{x1D650}$ . The reader is referred to Giannetti & Luchini (Reference Giannetti and Luchini2007), Pralits, Brandt & Giannetti (Reference Pralits, Brandt and Giannetti2010) and Marquet & Lesshafft (Reference Marquet and Lesshafft2015) for implementation details.

Second, the endogeneity $E = E_r + \mathrm{i} E_i$ , defined as

(3.10) \begin{equation} E(z) = \tilde {\boldsymbol{q}}^\dagger \boldsymbol{\cdot }\left ({\unicode{x1D63C}} \tilde {\boldsymbol{q}} \right )\!, \end{equation}

where the dot indicates scalar multiplication of two vectors, provides a local budget for eigenvalue selection by quantifying the pointwise contribution of the unperturbed operator acting on the direct mode, weighted by the adjoint mode. In its discrete form, it yields a spatial density whose integral recovers the eigenvalue, namely,

(3.11) \begin{equation} \int _0^L E(z) \, \mathrm{d}z = \omega , \end{equation}

and whose real and imaginary parts can be interpreted as contributions to growth-rate and frequency selection, respectively (Marquet & Lesshafft Reference Marquet and Lesshafft2015). Accordingly, $\textrm{Re} (E)$ can be read as a local contribution to amplification or decay, while $\textrm{Im} (E)$ gives the local contribution to frequency selection. In the present discrete implementation, the vector ${\unicode{x1D63C}}\tilde {\boldsymbol q}$ is naturally partitioned into the four equation blocks associated with continuity, momentum and the two polymer-stress transport equations, thus one can write

(3.12) \begin{equation} {\unicode{x1D63C}}\tilde {\boldsymbol q} = \begin{bmatrix} ({\unicode{x1D63C}}\tilde {\boldsymbol q})_{\textit{cont}}\\[3pt] ({\unicode{x1D63C}}\tilde {\boldsymbol q})_{\textit{mom}}\\[3pt] ({\unicode{x1D63C}}\tilde {\boldsymbol q})_{\tau _{\textit{zz}}}\\[3pt] ({\unicode{x1D63C}}\tilde {\boldsymbol q})_{\tau _{\textit{rr}}} \end{bmatrix}. \end{equation}

Accordingly, the total endogeneity is decomposed as

(3.13) \begin{equation} E(z)=E_{\textit{cont}}(z)+E_{\textit{mom}}(z)+E_{\tau _{\textit{zz}}}(z)+E_{\tau _{\textit{rr}}}(z), \end{equation}

where

(3.14) \begin{equation} E_{\textit{cont}}=\tilde r^{\dagger }({\unicode{x1D63C}}\tilde {\boldsymbol q})_{\textit{cont}},\quad E_{\textit{mom}}=\tilde u^{\dagger }({\unicode{x1D63C}}\tilde {\boldsymbol q})_{\textit{mom}}, \quad E_{\tau _{\textit{zz}}}=\tilde \tau _{\textit{zz}}^{\dagger }({\unicode{x1D63C}}\tilde {\boldsymbol q})_{\tau _{\textit{zz}}},\quad E_{\tau _{\textit{rr}}}=\tilde \tau _{\textit{rr}}^{\dagger }({\unicode{x1D63C}}\tilde {\boldsymbol q})_{\tau _{\textit{rr}}}. \end{equation}

Together, wavemaker and endogeneity allow us to distinguish where the flow is most sensitive to localised perturbations (wavemaker) from where the eigenvalue is intrinsically produced by the dynamics of the unperturbed base state (endogeneity), thereby providing a sharp picture of the regions controlling global instability in stretched Newtonian and viscoelastic jets. In the present context, this distinction is especially useful because the region where the jet is most receptive need not coincide exactly with the region where the capillary, inertial and viscoelastic contributions combine to select the eigenvalue.

4. Selected case studies

This section presents a set of representative case studies designed to (i) benchmark the present base-flow and global-stability formulation against established Newtonian results and (ii) illustrate, in a systematic manner, how weak elasticity modifies both the onset and the structure of the global oscillatory instability. We focus on parameter combinations that lie close to marginality so that (a) the leading eigenpair is cleanly separated from the rest of the spectrum and (b) small changes in rheology produce measurable eigenvalue drift. The first case reproduces the nearly marginal Newtonian jetting configuration of Rubio-Rubio et al. (Reference Rubio-Rubio, Sevilla and Gordillo2013) and serves as a reference for our numerical implementation (base state, spectrum and eigenfunctions), after which we introduce viscoelasticity at fixed viscosity ratio to quantify the role of polymeric tension in the linear dynamics. The second case is a nearly marginal viscoelastic jet for which we validate the global-mode prediction by direct time marching of the nonlinear one-dimensional equations under a weak, transient inlet perturbation. Together, these two cases establish the reliability of the global framework and provide the building blocks for the mechanistic assessment developed below (wavemaker, endogeneity and their decomposition), as well as for the parametric trends in marginal conditions discussed subsequently.

Figure 1.

(a) Direct spectrum, where hollowed out squares mark the leading eigenvalues, (b) spatial evolution of the base polymeric stresses and (c) base jet radius, for $\textit{We} = 0.003$ , ${\textit{Bo}} = 1.81$ , $\varGamma = 5.83$ and selected rheology parameters (see legend). Dark blue crosses in (a) and (c) are Newtonian results from Rubio-Rubio et al. (Reference Rubio-Rubio, Sevilla and Gordillo2013), included for comparison.

Figure 1. Long description

The image contains three graphs. The first graph (a) is a scatter plot showing the direct spectrum with hollowed-out squares marking the leading eigenvalues. The x-axis represents the real part of the eigenvalues, and the y-axis represents the imaginary part. The second graph (b) consists of two line graphs showing the spatial evolution of the base polymeric stresses. The x-axis represents the spatial coordinate z, and the y-axis represents the stress components. The third graph (c) is a line graph showing the base jet radius as a function of the spatial coordinate z. The x-axis represents the spatial coordinate z, and the y-axis represents the jet radius. Dark blue crosses in graphs (a) and (c) are Newtonian results from Rubio-Rubio et al. (2013), included for comparison. The graphs compare different rheology parameters, with legends indicating Newtonian, beta = 0.5, De = 3, and beta = 0.5, De = 5.

4.2. The case $\textit{We} = 7\times 10^{-4}$ , ${\textit{Bo}} = 1.8$ , $\varGamma = 5.8$ , $\beta = 0.5$ , $\textit{De} = 5$

Figure 2.

Marginally stable jet for $\textit{We} = 7\times 10^{-4}$ , ${\textit{Bo}} = 1.8$ , $\varGamma = 5.8$ , $\beta = 0.5$ , $\textit{De} = 5$ . (a) Direct spectrum, for which the leading mode (marked with a square) is $\omega \simeq -1.211\times 10^{-6} \pm 0.08984\mathrm{i}$ . (b) Comparison between time marching and the leading linear eigenmode dynamics. (c) and (d) Direct and adjoint normalised eigenfunctions.

Figure 2. Long description

The image contains four graphs analyzing a marginally stable jet. The first graph (a) is a scatter plot showing the direct spectrum with the leading mode marked by a square. The x-axis represents the real part of the frequency, and the y-axis represents the imaginary part. The second graph (b) consists of two line graphs comparing time marching and the leading linear eigenmode dynamics. The top line graph shows the mean radius over time, while the bottom line graph shows the maximum axial stress over time. The third set of graphs (c) displays the real and imaginary parts of the direct normalized eigenfunctions for radial displacement, axial velocity, and axial stress as functions of the axial position. The fourth set of graphs (d) shows the real and imaginary parts of the adjoint normalized eigenfunctions for radial displacement, axial velocity, and axial stress as functions of the axial position. The real parts are represented by solid lines, and the imaginary parts by dashed lines. All values are approximated.

This case is chosen to be nearly marginal (weakly damped) in the viscoelastic regime, with a leading complex-conjugate pair $ \omega \simeq -1.211\times 10^{-6}\pm 0.08984\,\mathrm{i},$ so that the growth rate is essentially zero on the time scales of interest while a well-defined global oscillation frequency is selected. Figure 2(a) shows the direct spectrum, where the dominant pair lies closest to the imaginary axis and is well separated from the remaining (more strongly damped) eigenvalues, indicating that the long-time linear dynamics is governed by a single global Hopf-type mode rather than by a broadband convective response, as in the Newtonian case. To validate this prediction in the fully nonlinear time-dependent equations, we perform a direct time marching of (2.9)–(2.12) using the steady base state $\boldsymbol{q}_0(z)$ as the initial condition, and we trigger the response via a small, transient inlet perturbation of the velocity

(4.1) \begin{equation} u(0,t) = \textit{We}^{1/2}{\textit{Bo}}^{-1/4}\big [1+\delta \,\mathrm{e}^{-a(t-t^*)^2}\,\mathcal{H}(t-t^*)\big ], \end{equation}

where $\mathcal{H}$ is the Heaviside function and $\delta =0.01$ is sufficiently small for the ensuing evolution to remain in the linear regime and qualitatively independent of the forcing amplitude. The parameters $a=10$ and $t^*=10^{-3}$ control the width and onset of the notch and are selected to facilitate numerical start-up while providing a localised impulse-like excitation of the global dynamics. As shown in figure 2(b), both the mean radius and the maximum axial polymer stress $\tau _{\textit{zz}}$ extracted from the time marching rapidly lock onto the linear prediction. After a short transient, the signals oscillate at the eigenfrequency $\omega _i$ with an envelope consistent with $\exp (\omega _r t)$ (here essentially constant owing to $|\omega _r|\ll 1$ ), thereby confirming the global-mode interpretation. Finally, the spatial structure of the associated eigenfunctions (figure 2 c) further supports this picture. The direct mode exhibits an extended downstream support (consistent with a convective redistribution of perturbations along the stretched jet), whereas the adjoint mode is strongly localised near the inlet, indicating that the receptivity and mode selection are controlled predominantly by the near-nozzle region. This marked upstream localisation of the adjoint fields underscores that small perturbations imposed at, or near, the inlet are the most effective at exciting the observed oscillatory response.

Figure 3.

Wavemaker and endogeneity spatial distributions for three separate cases: $(\textit{We}, {\textit{Bo}}, \varGamma ) = (3\times 10^{-3},1.81, 5.83)$ (Newtonian), $(\textit{We}, {\textit{Bo}}, \varGamma , \beta , \textit{De}) = (1.85 \times 10^{-3}, 1.8, 5.8, 0.6, 3)$ and $(\textit{We}, {\textit{Bo}}, \varGamma , \beta , \textit{De}) = (7 \times 10^{-4}, 1.8, 5.8, 0.5, 5)$ . For endogeneity computations, solid lines indicate real part, and dashed lines indicate imaginary part.

Figure 3. Long description

The image contains three separate graphs showing wavemaker and endogeneity spatial distributions. Each graph represents a different case, with the top graph for a Newtonian case, the middle graph for a case with specific parameters, and the bottom graph for another set of parameters. The x-axis represents the spatial variable z, while the y-axis represents the normalized values of the wavemaker and endogeneity distributions. Solid lines indicate the real part of the endogeneity computations, and dashed lines indicate the imaginary part. The graphs show how these distributions vary with z for each case. All values are approximated.

5. Wavemaker and endogeneity: the origin of instability

Figure 3 compares wavemaker and endogeneity diagnostics for three nearly marginal configurations, namely, a Newtonian reference case (top panels) and two viscoelastic cases of increasing elasticity (middle and bottom panels). In all three cases the parameters are chosen close to the onset of global instability, so that the leading eigenpair is weakly damped (or nearly neutral) and the associated spatial diagnostics can be interpreted as the footprint of the global feedback loop. Throughout, solid and dashed lines denote the real and imaginary parts of the endogeneity, respectively; accordingly, $\textrm{Re} (E)$ quantifies the local contribution to the global growth/decay rate, while $\textrm{Im} (E)$ quantifies the local contribution to the global oscillation frequency. For the Newtonian case, the normalised wavemaker $S(z)$ is sharply localised near the inlet and decays rapidly over a few nozzle radii, while $S_{\textit{mom}}(z)$ is nearly indistinguishable from $S(z)$ , indicating that the feedback region is compact and that restricting the structural-sensitivity tensor to the $u$ -row captures essentially the same spatial support. Consistently, the endogeneity is dominated by the continuity/capillary pathway, i.e. $E_{\textit{cont}}$ is $O(10^{-1})$ locally and displays the expected sign structure in its real part (a positive lobe followed by a compensating negative lobe), so that the net integrated contribution remains small in the marginal limit, whereas $E_{\textit{mom}}$ is negligible in both real and imaginary parts. The imaginary component $\textrm{Im} (E_{\textit{cont}})$ is also confined to the upstream region, supporting the interpretation that the onset frequency is set primarily by the near-nozzle feedback core.

The viscoelastic cases retain the same qualitative global-oscillator signature but reveal a progressive and physically meaningful change in the spatial organisation of the sensitivity as elasticity is increased. For the intermediate case ( $\beta =0.6$ , $\textit{De}=3$ , middle panels), the wavemaker remains upstream peaked but its decay is already noticeably slower than in the Newtonian limit, producing a downstream tail that extends over a significantly larger fraction of the domain. This trend is further amplified in the more elastic case ( $\beta =0.5$ , $\textit{De}=5$ , bottom panels), for which $S(z)$ becomes markedly more distributed and remains non-negligible over an extended downstream interval. In other words, one of the most robust effects of introducing a viscoelastic dynamics in a gravitationally stretched jet is not merely a shift of the leading eigenvalue, but a pronounced spreading of the feedback-sensitive region. The instability is still selected upstream, yet the portion of the jet that participates in the global feedback loop grows with $\textit{De}$ . This behaviour is consistent with the fact that polymeric stresses are advected and relax over a finite time scale, so that elastic memory couples distant axial locations through the transport of tensile stress. At the same time, the separation between $S(z)$ and $S_{\textit{mom}}(z)$ becomes increasingly apparent as elasticity grows (most clearly in the bottom case), for which $S(z)$ exceeds $S_{\textit{mom}}(z)$ over much of the downstream tail, indicating that a momentum-row restriction underestimates the full structural sensitivity once viscoelastic degrees of freedom are active and confirming that the extra sensitivity is carried by the stress subspace rather than by the Newtonian momentum balance alone.

The endogeneity decomposition confirms this interpretation while highlighting that the dominant balance remains localised closer to the inlet than the wavemaker itself. In both viscoelastic cases, $E_{\textit{mom}}$ remains small, but an additional contribution associated with polymeric tension emerges through $E_{\tau _{\textit{zz}}}$ . This elastic contribution overlaps with the region where $E_{\textit{cont}}$ is active and provides a clear mechanism for eigenvalue selection. The onset results from a coupled capillary–elastic feedback pathway in which kinematic/capillary thinning is opposed by an additional tensile-stress response. The real parts of $E_{\textit{cont}}$ and $E_{\tau _{\textit{zz}}}$ exhibit comparable magnitudes in the upstream region and partially compensate one another, which is precisely the structure expected near marginality, where the global growth rate is the small residual of competing local contributions. Meanwhile, the imaginary parts (dashed curves) remain comparatively confined to the upstream region for all three cases, indicating that the oscillation frequency is still set predominantly by the near-inlet feedback core even though viscoelasticity distributes the sensitivity farther downstream. Taken together, figure 3 supports two conclusions. First, all three configurations behave as genuine global oscillators, i.e. both wavemaker and endogeneity decay well before the outlet, indicating that the selected mode is not an artefact of domain truncation. Second, viscoelasticity modifies the instability mechanism in a specific and important way: as $\textit{De}$ increases, the feedback-sensitive region becomes progressively more widespread while a distinct polymeric endogeneity contribution appears, showing that the marginal state is controlled by a spatially distributed capillary–elastic feedback loop rather than by a purely Newtonian, near-nozzle mechanism.

Figure 4.

Spatial metrics of the wavemaker over the $(\beta , \textit{De})$ parametric plane for $(\textit{We}, {\textit{Bo}}, \varGamma ) = (7 \times 10^{-4}, 1.8, 5.8)$ .

Figure 4. Long description

The image contains four graphs depicting spatial metrics of the wavemaker over the parametric plane. The top three graphs are line plots showing the variation of a parameter S with respect to z for different values of De (Deborah number) and beta (β). The first graph represents β equals 1, Newtonian, the second graph represents β equals 0.8, and the third graph represents β equals 0.4. Each graph includes multiple lines corresponding to different De values ranging from 0 to 6. The fourth graph is a heat map showing the variation of zS and Supp(S) with respect to De and β. The heat map uses a color scale to represent the values. The inset in the third graph shows a detailed view of the support of S. All values are approximated.

To support these claims, figure 4 quantifies the redistribution of the wavemaker with elasticity using both the full profiles $S(z)$ and the scalar metrics

(5.1) \begin{equation} z_S=\operatorname *{arg\,max}_z S(z), \qquad \mathrm{supp}(S)=\int _0^L \frac {S(z)}{\max _z S(z)}\,\mathrm{d}z, \end{equation}

for fixed $(\textit{We}, {\textit{Bo}}, \varGamma ) = (7 \times 10^{-4}, 1.8, 5.8)$ across the $(\beta , De)$ parametric plane (which do not necessarily correspond to marginal cases). The Newtonian reference case exhibits a sharply localised wavemaker confined to the immediate near-nozzle region. When viscoelasticity is introduced, the dominant effect is not an abrupt displacement of the peak away from the inlet, but a progressive broadening of the downstream tail of $S(z)$ as $\textit{De}$ increases, which is already visible for $\beta =0.8$ and becomes much more pronounced for $\beta =0.4$ . The heat maps confirm this trend quantitatively: the support length grows strongly as elasticity is increased, i.e. as $\textit{De}$ increases and $\beta$ decreases, whereas the peak location $z_S$ is comparatively less sensitive and shifts appreciably downstream only in the most elastic cases. Thus, viscoelasticity does not primarily relocate the instability core; rather, it preserves an upstream-selected feedback region while enlarging the axial portion of the jet over which tensile-stress memory participates in the global feedback loop. This provides a quantitative confirmation for the interpretation drawn from the wavemaker plots alone, namely that the main viscoelastic effect is a broadening of the feedback-sensitive region rather than a simple translation of its maximum.

Additionally, figure 5 quantifies this trend by reporting the fraction of total endogeneity activity carried by each equation, using the absolute budgets of $\textrm{Re} (E)$ and $\textrm{Im} (E)$ so as to measure the contribution of each equation into the endogeneity budget/eigenvalue production. In the Newtonian limit, both growth-rate and frequency selections are overwhelmingly dominated by $E_{\textit{cont}}$ , with only a minor contribution from $E_{\textit{mom}}$ . As $\textit{De}$ increases, however, a distinct axial polymer-stress pathway emerges through $E_{\tau _{\textit{zz}}}$ . Its relative weight grows systematically, especially in the growth-rate budget, while $E_{\tau _{\textit{rr}}}$ and $E_{\textit{mom}}$ remain comparatively small. Thus the viscoelastic modification of the instability is not a generic stress effect, but is carried primarily by the axial tensile response. At the same time, the continuity contribution remains dominant in the frequency budget, indicating that elasticity alters growth-rate selection more strongly than frequency selection.

Figure 5.

Endogeneity budget for $(\textit{We}, {\textit{Bo}}, \varGamma , \beta ) = (7 \times 10^{-4}, 1.8, 5.8, 0.5)$ and increasing $\textit{De}$ . The Newtonian case $\beta = 1$ is included as a baseline reference.

Figure 5. Long description

The bar graph compares the fraction of total activity for different components across various cases. The x-axis represents different cases, including Newtonian and De values ranging from 1 to 6. The y-axis represents the fraction of total activity. The graph is divided into two sections: Growth-rate-selection activity on the left and Frequency-selection activity on the right. Each bar is stacked and color-coded to represent different components: E tau r in purple, E tau zz in yellow, E mom in orange, and E cont in blue. The Newtonian case serves as a baseline reference. The data shows variations in the fraction of total activity for each component across the different cases. All values are approximated.

6. Influence of rheology on marginal conditions

Figure 6 generalises the marginal-curve analysis by varying the viscous level over three representative Kapitza numbers, $\varGamma =0.84$ , $5.83$ and $18.5$ , while keeping the viscosity ratio fixed at $\beta =0.8$ and sweeping $\textit{De}$ . Each row reports the critical Weber number $\textit{We}_c({\textit{Bo}})$ at the onset of the self-excited global oscillation and the corresponding critical frequency $\omega _{i,c}({\textit{Bo}})$ . For all three values of $\varGamma$ , the Newtonian reference curve decreases monotonically with $\textit{Bo}$ in both panels. Increasing the nozzle size, and therefore strengthening gravity relative to surface tension, shifts the transition to progressively smaller flow rates and to lower oscillation frequencies. This indicates that the global feedback loop becomes slower as the jet is more strongly stretched downstream. The middle row, $\varGamma =5.83$ , includes the Newtonian data of Rubio-Rubio et al. (Reference Rubio-Rubio, Sevilla and Gordillo2013) and recovers both the marginal branch and the onset frequency quantitatively, thereby extending the Newtonian benchmark of § 4.1 to the full marginal curve.

Figure 6.

Marginal curves $\textit{We}_c({\textit{Bo}})$ and corresponding frequencies $\omega _{i,c}({\textit{Bo}})$ for $\varGamma = 0.84, 5.83$ and $18.5$ . Dots represent individually computed marginal cases, joined by interpolating splines.

Figure 6. Long description

The image contains six line graphs arranged in a two-by-three grid. Each graph plots marginal curves and corresponding frequencies for different conditions. The x-axis represents the Bond number (Bo) on a logarithmic scale, while the y-axis represents different variables: the critical Weber number (We_c) on the left column and the critical frequency (ω_i,c) on the right column. The graphs compare Newtonian fluids and viscoelastic fluids with varying Deborah numbers (De). Dots represent individually computed marginal cases, joined by interpolating splines. The top row corresponds to a lower Γ value, the middle row to a moderate Γ value, and the bottom row to a higher Γ value. Each graph includes multiple lines representing different conditions: Newtonian fluid, and viscoelastic fluids with β = 0.8 and De values of 2, 4, and 6. The trends show how the critical Weber number and critical frequency change with the Bond number for different fluid properties.

Superimposed on this Newtonian baseline, viscoelasticity produces a systematic ordering with $\textit{De}$ across the whole $({\textit{Bo}},\varGamma )$ range considered. For fixed $({\textit{Bo}},\varGamma )$ , increasing $\textit{De}$ shifts both $\textit{We}_c$ and $\omega _{i,c}$ downward, so that $\textit{We}_c(\textit{De}=2)\gt \textit{We}_c(\textit{De}=4)\gt \textit{We}_c(\textit{De}=6)$ and $\omega _{i,c}(\textit{De}=2)\gt \omega _{i,c}(\textit{De}=4)\gt \omega _{i,c}(\textit{De}=6)$ , with the viscoelastic curves lying below the Newtonian one. In the present formulation this trend is interpreted in terms of an additional tensile-stress pathway. As $\textit{De}$ increases at fixed $\beta$ , polymeric stresses are transported over longer axial distances before relaxing. This strengthens the base-state tensile resistance and enhances the linear coupling by which stress perturbations feed back on the free surface through the axial-tension term. In effect, elastic tension competes with capillary-driven thinning and delays the onset of the global oscillatory instability, so that marginality is reached only when the through flow is reduced, meaning smaller $\textit{We}$ . The reduction in $\omega _{i,c}$ with $\textit{De}$ is consistent with the same picture. Marginality occurs at weaker advection and the polymer introduces a memory time scale that adds phase lag to the global loop. Both effects promote slower oscillations at onset. This interpretation is consistent with the wavemaker and endogeneity diagnostics of § 5, which show that elasticity broadens the feedback-sensitive region and introduces a non-negligible polymeric contribution to eigenvalue production in marginal configurations.

Varying $\varGamma$ modulates these trends in a physically informative way. Increasing $\varGamma$ shifts the overall Newtonian boundary downward in $\textit{We}_c$ and in $\omega _{i,c}$ , which reflects the stabilising role of viscous stresses in suppressing capillary-driven unsteadiness and permitting sustained jetting at lower flow rates. The incremental viscoelastic shift is noticeably weaker for the lowest viscous level, $\varGamma =0.84$ , where the curves nearly collapse and the dependence on $\textit{De}$ becomes mild. This near collapse is consistent with a regime in which onset is controlled by a rapid inertio-capillary balance localised near the inlet, leaving limited opportunity for polymeric stresses to accumulate and influence the marginal condition. It is also the regime where Rubio-Rubio et al. (Reference Rubio-Rubio, Sevilla and Gordillo2013) emphasise that additional upstream physics can become important for low-viscosity liquids issued through long injectors, in particular viscous relaxation of the exit velocity profile (see also Sevilla (Reference Sevilla2011)). By way of contrast, for $\varGamma =5.83$ and $18.5$ the viscoelastic shift is clearly resolved, especially at small-to-moderate $\textit{Bo}$ , and it weakens as $\textit{Bo}$ increases. This is consistent with the increasing dominance of gravity-driven stretching and advection at large $\textit{Bo}$ , which drives the marginal curves back towards the Newtonian scaling. Overall, figure 6 shows that moderate elasticity, $\textit{De}=O(1)$ , induces systematic shifts of both the critical flow rate and the onset frequency, and that these shifts depend on viscosity in a manner consistent with a coupled capillary–elastic feedback loop whose influence is strongest when the jet is not too strongly gravity dominated.