The polymer diffusive instability in highly concentrated polymeric fluids

The extrusion of polymer melts is known to be susceptible to `melt fracture' instabilities, which can deform the extrudate, or cause it to break entirely. Motivated by this, we consider the impact that the recently discovered polymer diffusive instability (PDI) can have on polymer melts and other concentrated polymeric fluids using the Oldroyd-B model with the effects of polymer stress diffusion included. Analytic progress can be made in the concentrated limit (when the solvent-to-total-viscosity ratio $\beta \rightarrow 0$), illustrating the boundary layer structure of PDI, and allowing the prediction of its eigenvalues for both plane Couette and channel flow. We draw connections between PDI and the polymer melt `sharkskin' instability, both of which are short wavelength instabilities localised to the extrudate surface. Inertia is shown to have a destabilising effect, reducing the smallest Weissenberg number ($W$) where PDI exists in a concentrated fluid from $W\sim 8$ in inertialess flows, to $W \sim 2$ when inertia is significant.


Introduction
The extrusion of plastic is a widespread industrial process in which plastics are shaped for countless applications.Plastic is melted and then forced through a slit (the 'die'), where it can then be shaped into sheets, tubes or moulded precisely.This molten plastic is a type of viscoelastic 'polymer melt', in which the fluid consists of long chain polymers, causing the fluid to feel an added force due to polymer stress.
It is known that the extrusion of polymer melts is an unstable process, as the extrudate can become distorted, or break entirely in a process called 'melt fracture' [1][2][3][4][5].This can place a limit on the speed of industrial plastic extrusion, slowing down the manufacturing process.Depending on the extrusion speed, different instabilities can be seen.The 'sharkskin' instability is characterised by the extruded plastic developing a regular short-wavelength distortion, and it is thought that its origins lie with the interaction of the flow with the outlet of the channel [6,7].The 'spurt-flow' instability alternates between the sharkskin and smooth textures, and the proposed mechanism for this is due to wall-slip [7][8][9].Even if these instabilities were suppressed, there is evidence that the pipe flow of polymer melts are intrinsically unstable due to a weakly non-linear subcritical bifurcation [3,4] that makes the system unstable to finite amplitude perturbations.This non-linear work was motivated by the apparent lack of a linear instability in polymer melts.
The recently discovered polymer diffusive instability (PDI) [10] is a linear instability that exists in an Oldroyd-B fluid and is localised at boundaries.It is elastic in origin, existing at vanishing Reynolds number, and requires the presence of polymer stress diffusion in the model.PDI provides a potential mechanism for the transition to elastic turbulence, a chaotic state characterised by a wide range of length scales in the absence of inertia, in planar geometries.The mechanism for elastic turbulence in curvilinear geometries is well understood -tension acts along curved streamlines which produces hoop stresses, generating a linear instability [11].However, in planar geometries there are no curved streamlines, and so this linear instability does not exist.Consequently, the process responsible for the transition to elastic turbulence in such geometries is currently under debate.Numerous mechanisms have been suggested, including linear instabilities [10,[12][13][14], finiteamplitude perturbations [15,16] and transient growth effects [17][18][19][20][21]. PDI represents one possible mechanism, and it has already been shown to trigger a transition from plane Couette flow to a sustained 3D chaotic state in a FENE-P fluid [10].
While previous work has investigated PDI in dilute fluids [10], we consider the concentrated limit of an Oldroyd-B fluid (where the solvent-to-total-viscosity ratio β → 0).This corresponds to highly concentrated polymeric solutions, as well as polymer melts.In §2 we examine the governing equations, and demonstrate that analytic progress can be made when β → 0 for both plane Couette flow and channel flow.We find the dispersion relation, identify the relevant scalings and detail the boundary layer structure present in PDI.In §3 we compare the analytic predictions of the growth rate to numerical results for plane Couette flow.We consider both inertialess and inertial flows as β → 0, and find that the growth rate (in units of the applied shear) diverges.In this limit, the smallest Weissenberg number, W, at which the system is unstable is W crit ∼ 8 for inertialess flows, and W crit ∼ 2 when inertia is significant.Numerically, we consider Weissenberg numbers of 1 ≲ W ≲ 1, 000, and Reynolds numbers 0 ≤ Re ≲ 100, 000 as well as solvent to total viscosity ratios 0.001 ≲ β.The order of magnitude of this last parameter is motivated by the use of the UCM fluid to model melt fracture [3,7,9,22], which assumes that β is sufficiently small that it can be neglected.In section §4 we check the robustness of PDI by changing the stress boundary conditions from vanishing diffusion to the Neumann condition of vanishing stress flux.Finally, there is a discussion in §5 in which we draw connections to the phenomena of melt fracture.

Formulation
We consider the flow of an Oldroyd-B fluid in 2D, with flow direction x and wall normal direction ŷ in a channel either driven by boundary motion (plane Couette flow) or an applied pressure gradient (channel flow).The effects of polymer stress diffusion are included, motivated by its presence in viscoelastic simulations to regularise hyperbolic constitutive relationships [23][24][25][26][27][28].

Governing Equations
The non-dimensional governing equations (e.g.[10]) are where T and C are the polymeric stress tensor and conformation tensor respectively.In the Oldroyd-B model they are simply related as follows In the above, we have scaled each quantity using the channel half-height h, a characteristic speed U 0 , and the total viscosity µ = µ p + µ s which sums the polymer and solvent viscosities (µ p and µ s respectively).In plane Couette flow, U 0 is taken to be the speed of the plates (u(x, y = ±h, t) = ±U 0 x) while in channel flow U 0 is defined in terms of the imposed pressure gradient along the channel (U 0 = −h 2 /µ dP/dx).There are 4 dimensionless parameters after scaling, where additionally we define ρ as the fluid density, λ as the relaxation timescale of the polymer, and δ as the polymer stress diffusivity.These dimensionless parameters are the Reynolds number Re, the Weissenberg number W , the solvent-to-total-viscosity ratio β and the polymer stress diffusion coefficient ε.

Boundary Conditions
All equations will be solved subject to no-slip and no-penetration conditions on the velocity at the boundaries, yielding four boundary conditions for our 2D system.

Channel flow
While there is evidence that wall-slip occurs for polymer melts [1,[7][8][9]29], it is not required for the existence of PDI, and so we do not include this in our model.We consider two types of boundary conditions on the polymer stress: The Type I condition corresponds to the polymer stress not diffusing at the walls.This is equivalent to the more typical, but less concise condition of applying equation ( 2) with ε = 0 at the wall [23,24], provided T has a continuous second derivative.The Type II Neumann condition instead corresponds to the stress flux vanishing at the walls.While the Type I condition is physically more relevant, we introduce the Type II condition as a way to check robustness to the choice of boundary conditions.We focus on Type I conditions, and show how the results change in §4 when type II conditions are used instead.

Base States
To find the base states we solve equations ( 1)-( 4) with u = U (y)x, T = T(y) and p = P (x), subject to the above boundary conditions.

Plane Couette flow
We impose dP/dx = 0 for plane Couette flow.The solutions for both Type I and Type II stress boundary conditions are This is independent of the stress boundary conditions as this solution has a constant stress.

Channel flow
A pressure gradient of dP/dx = −1 is imposed for channel flow.The laminar base state is altered by the presence of polymer diffusion.For Type I boundary conditions the exact solutions are We can see this base flow contains a boundary layer with lengthscale √ εW .For Type II conditions we find a narrower boundary layer of lengthscale √ εβW , with All of these laminar base flows were checked by comparing the analytic results to the numerical results obtained using the open-source python package Dedalus and its non-linear boundary value problem solver [30].Note that the presence of boundary layers due to polymer stress diffusion could provide a way to verify the validity of this model experimentally.Both boundary conditions for channel flow have boundary layers of lengthscale √ εW , which can be sizable.For example, a long chain polymer melt with ε = 0.001 and W = 100 has a boundary layer of length ∼ 0.3h or 15% of the channel width.See figure 1 for how the flow variables can be modified due to polymer stress diffusion.

Linearised Equations
We now linearise equations ( 1)-( 4) about the base states u = U (y)x and T = T(y) from §2.3.We introduce the polymeric stress perturbation τ and the streamfunction perturbation ψ, and assume τ , ψ ∝ e ik(x−ct) where c ∈ C is an eigenvalue to be determined.This results in where D := d dy .Choosing to use the streamfunction formulation immediately satisfies the incompressibility equation ( 3).We take the curl of (1) to eliminate the pressure, and then linearise it to give (5) which is the viscoelastic analogue to the Orr-Sommerfeld equation [31].Using ( 2) and ( 4), we then obtain the linearised constitutive equation for the stress perturbation as where we suppress the lower off diagonal element for clarity as all tensors are symmetric.The analysis for plane Couette flow and channel flow with vanishing stress diffusion at the boundary is very similar and so we keep all equations sufficiently general to cover both systems.

Asymptotics
We now rescale our variables so that we can identify which terms in ( 5) and ( 6) can be neglected in the polymer melt limit.To motivate our scaling, we observe that for inertialess flow c ∼ εW (1 − β)/β and k ∼ W (1 − β)/εβ.Then, noting PDI is localised at the boundaries [10], we introduce a boundary variable ŷ at the lower boundary, and rescale quantities as follows Using these rescalings, as well as the Taylor expansion of U around y = −1 equations ( 5) and ( 6) become and In the limit of β/W (1 − β) → 0 for finite εRe/β and W these become (11) This derivation makes use of the fact that for both plane Couette and channel flow, allowing us to treat the base stress as constant in the region where PDI is active.The same simplification also occurred for the base velocity U , allowing U − c to be treated as constant in the given limit.This is now a system with constant coefficients.Applying the operator (−i kĉ + k2 − D2 ) to equation (10), and then using equation (11) to rewrite each τ in terms of ψ gives (12) This admits solutions of the form ψ = ψ0 e i lŷ , where l satisfies the dispersion relation (13) which relates the complex eigenvalue ĉ to the spatial wavenumber pair ( k, l).We now consider the superposition of these solutions in order to satisfy the boundary conditions.

Boundary Conditions
For a given ĉ and k, equation ( 13) yields 6 possible l that can contribute to the eigenmode, yielding ψ = 6 n=1 ψ0,n e i ln ŷ . ( We are considering PDI localised to the lower boundary and so we can discard all l that do not decay as ŷ → ∞.This leaves only l with positive imaginary parts.One of these is an exact solution, l = i k, which appears as PDI is acting on the vorticity field, ( D2 − k2 )ψ = ∇2 ψ, and this component of the eigenmode has zero vorticity.We assume two more l exist with positive imaginary part, l+ , l− .We therefore have Using the no-penetration and no-slip conditions results in B ± = Af ± ( k, l+ , l− ) where f ± can be written analytically.We can then substitute the resultant ψ into equation ( 11) to obtain a 2nd order differential equations for each component of τ .
where the g i can be found exactly.Solving this differential equation introduces a new scale in the stress field, with decay rate ± k2 − i kĉ (the sign of this is set so that the resultant term decays in the far field).To simplify our notation, we will define the complex square root to pick the solution with a positive real part.Using the Type I stress boundary conditions then sets the coefficient of this introduced complementary function so that where the h i can be found analytically.We substitute ψ and τ into equation (10), and consider the resultant coefficients of each exponential on the left-hand side.The coefficient of exp − k ŷ trivially vanishes, while those of exp i l± ŷ just recover that the dispersion relation ( 13) is satisfied for l± .
Examining the coefficient of exp − k2 − i kĉ ŷ however gives a further condition.For the cases of plane Couette and channel flow with Type I stress boundary conditions, this extra condition is Solving equations ( 13) and ( 18) together gives the eigenvalue ĉpred .This was done by using the Newton-Raphson method to find zeros of L(ĉ) = L(ĉ, l+ (ĉ), l− (ĉ)), where l± (ĉ) are the solutions to equation ( 13) with positive imaginary part.In summary, the asymptotic limit pursued here has turned the original system with 4 parameters Re, β, ε and W into a constant-coefficient eigenvalue problem with only 2 parameters W and εRe/β, i.e.
In particular, the parameter ε now only appears in the parameter group εRe/β which measures the strength of inertia.

Results
We now compare our analytic predictions to numerical computations in the polymer melt limit.All numerical results are computed using the open-source python package Dedalus [30].As the results for plane Couette and channel flow are similar, we discuss plane Couette in the main text and place the key results for channel flow in Appendix A. We firstly consider creeping flows (Re = 0), and then investigate how this flow is altered due to the presence of inertia (Re > 0).

Inertialess flows: Re = 0
We begin by verifying the validity of the predicted scalings.To help with this we define the scaled and unscaled growth rates where R denotes taking the real part.Hence In addition, we define σ * and σ * to be the scaled and unscaled most unstable growth rates across all wavenumbers.Figure 2a shows the dispersion relations of σ against k as β varies, demonstrating that our analytic predictions match the numerics in the polymer melt limit.It also suggests that k = O(1) and σ = O(1) as β → 0, verifying our predicted scalings.Figure 2b shows that σ * tends to a constant, and that this limiting value of σ * is non-monotonic in W .The maximum growth rate σ * = W (1 − β)/β σ * therefore diverges at small β, as σ = O(1) as β → 0. We emphasise that the timescale on which this divergence occurs is the inverse of the applied shear rate h/U 0 so this divergence is not a feature of the choice of timescale.Physically, this means that the instability becomes stronger the more concentrated the polymer melt is and suggests that PDI is particularly of interest for concentrated polymeric fluids.
Figure 3 shows the neutral stability curves in (W , k) space when β or ε change.When Re = 0, the derived analytic results show that ĉ depends only on W and k in the small β/W (1 − β) limit (that is, (19) with εRe/β = 0), and hence the stability of the system is entirely shown by this plot.We show in figure 3a how the neutral curve changes as β is taken out of this limit.The analytic theory reproduces the β = 0.01 neutral curve almost exactly, and remains qualitatively accurate for viscosity ratios as large as β = 0.4.Similarly, in figure 3b we increase ε.The analytic theory is accurate for physically realistic values of polymer diffusion, with it remaining accurate for values as large as ε = 0.1.The reason the neutral curve stops being accurate for larger ε is due to the delocalisation of the eigenfunction near the boundary, when the analytic theory assumed localisation.The boundary layer thickness scales like εβ/W (1 − β) and so when ε is too large the eigenfunction interacts with the other wall.
The extrusion of plastics is typically modelled as the inertialess flow of a polymer melt [3,6,9], and so the above results are relevant to this system.We emphasise the result that the growth rate of PDI diverges in the limit of small β in units of U 0 /h, meaning that polymer melts are especially susceptible to this instability.We discuss the connection between PDI and the 'melt fracture' instabilities that are known to affect plastic extrusion in section 5 and argue there that the effects of PDI are potentially significant.

Flows with Inertia: Re > 0
We now introduce inertia and consider how this alters the instability.The analytic work demonstrated that Re only enters the equations in the parameter group εRe/β, and so this controls the strength of inertia.It will be useful to define to be the Reynolds number at which non-inertial behaviour transitions to inertial behaviour.The scalings and analysis already presented considered the limit of β/W (1−β) → 0 when Re/Re * is finite, however in this section we consider the effects of increasing Re/Re * with finite β/W (1 − β) ≪ 1.We use Re/Re * ∼ 1000 as the largest magnitude of this parameter.Bounding this parameter means it is finite, and hence provided the value of β/W (1 − β) is sufficiently small, the analytic predictions already presented remain valid.However, now that Re/Re * is no longer fixed, we incorporate it into our scalings to understand how k and c vary with it.
In figure 4a we consider the relationship between the most unstable wavenumber k * and Re/Re * .The most obvious feature of this plot is that k * is discontinuous, with it jumping to an inertial scaling of k * ∼ (Re/Re * ) 1/2 at around Re ∼ 50Re * .This transition point occurs at larger Re/Re * as W is increased.This discontinuity is due to the eigenmode having different behaviours at k ∼ (Re/Re * ) 0 and k ∼ (Re/Re * ) 1/2 .We examine this by considering the dispersion relation as we increase Re/Re * in figure 5.As Re/Re * increases, the instability in the non-inertial region of k ∼ (Re/Re * ) 0 is suppressed, while it is promoted in the inertial region of k ∼ (Re/Re * ) 1/2 .The discontinuity in k * appears when the inertial region becomes more unstable than the non-inertial region.
Figure 4b shows the dependence of the maximum growth rate σ * on Re/Re * , and it demonstrates that σ ∼ (Re/Re * ) 0 .In the limit of β → 0, for fixed Re, W and ε, the maximum growth rate curves collapse when plotted against Re/Re * = εRe/β showing that σ * = σ * (Re/Re * , W ) in this limit, as suggested by the analytic predictions.The magnitude of this growth rate is not sensitive to the magnitude of Re/Re * .A kink is seen in the maximum growth rate curves when the most unstable wavenumber transitions from non-inertial to inertial scalings.
Hence, for Re/Re * ≫ 1, the variables scale like Figure 6a shows how the neutral curve in the (W, k) plane changes when inertia is present.Increasing inertia destabilises the system, as more of the (W, k) plane becomes unstable.In particular, the smallest W at which the system is unstable for some wavenumber (denoted by W crit ) decreases as Re increases.As Re is increased, the neutral curve starts to bend upwards at Re ≈ 13Re * .This bend means that, for fixed W , there is a region of stable wavenumbers in between regions of unstable wavenumbers, as was seen for the Re = 2000 curve in figure 5.As Re increases past Re ≈ 13Re * , the unstable region of the neutral curve continues to expand, so that the stable region between the bends in the neutral curve gradually disappears (not shown).
In figure 6b we investigate how W crit changes with inertia, which is measured by the magnitude of εRe/β.There are two distinct regimes determining how low W crit can be.When Re ≪ Re * we get W crit ∼ 8, which is the same result obtained for β ≈ 0.5 at Re = 0 in [10].However, when Re ≫ Re * , this critical value drops to W crit ∼ 2. This behaviour continues for viscosity ratios as large as β = 0.2, as W crit still has a marked drop at Re ∼ Re * .
Lastly, we show neutral curves in the (W , β) plane in figure 7 where we consider the stability across all wavenumbers.This shows that inertia generically decreases W crit across the whole range of β.In the small β limit, this plot demonstrates again that W crit ∼ 8 when there is no inertia, but that including inertia lowers the critical Weissenberg number to W crit ∼ 2. The predicted neutral curves match the numerics as β → 0. In the inertialess case of Re = 0, the prediction for the inertialess case even appears to be accurate for viscosity ratios as high as β = 0.3.These neutral curves demonstrate that plane Couette flow is susceptible to PDI over a large range of parameters.This shows that the first peak occurs at non-inertial wavenumbers k ∼ (Re/Re * ) 0 and the growth in this region is suppressed as inertia increases.In contrast, the second peak is promoted by inertia, and exists at wavenumbers k ∼ (Re/Re * ) 1/2 .The plots also show that σ ∼ (Re/Re * ) 0 .

Prediction Validity
We now quantify the accuracy of the predictions based on the maximum growth rate σ * .Figure 8a shows how the error in this prediction varies with β for fixed W, Re and ε.It demonstrates that   W , k) plane with changing Re with ε = 0.0001 and β = 0.01.These set the transition from non-inertial to inertial behaviour to be at Re * = 100.Darkest to lightest shade shows numerics for Re = 0, 100, 200, 400, 800, 1200, 1300, 1400.Dotted lines denote the neutral curve predicted by the analytic theory.Note when inertia is considered the prediction is less accurate -we expect that if a smaller β were considered that the prediction would match the numerics more closely.b) The relationship between W crit (the smallest W at which the system goes unstable) and inertia.Each dashed curve tracks changing Re with constant ε = 0.0001 and (darkest to lightest) constant β = 0.001, 0.01, 0.1, 0.2.The transparent grey curve shows the analytic prediciton.Note each curve has constant Re, ε so as β → 0 we do not have εRe/β finite, which was required for the analytic work to be valid.This explains why the dotted prediction curves do not match the solid numeric curves exactly in the small β limit.These plots demonstrate that PDI is generic across β, and that inertia destabilises the system.
the error decays as β → 0, and it motivates defining β lim (Re) such that when β < β lim (Re), the predicted growth rate is within 5% of the actual value.β lim is independent of ε, as the exact rescaled governing equations ( 8) and ( 9) are independent of ε for plane Couette flow.To find β lim at a given Re, we consider the growth rate error curves for W ∈ {10, 12, 20, 50, 100, 200, 400}, and take β lim to be the largest β such that all errors on all curves are smaller than 5% for all β < β lim (Re).To illustrate this, figure 8a shows that β lim (Re = 0) = 0.090.Figure 8b shows how β lim varies with Re, showing that our predictions require smaller β to be valid as Re increases.The plot suggests that any β < 0.02 is sufficient for predictions to be accurate to within 5% for Re < 1000.

Robustness to Boundary Conditions
So far we have considered the results of plane Couette flow with vanishing diffusion at the boundary (∇ 2 T| y=±1 = 0).We consider here how changing the boundary condition to vanishing stress flux (DT| y=±1 = 0) changes PDI, and we find that the instability is sensitive to this change.The laminar base flow does not change with the boundary condition in the case of plane Couette flow whereas a new boundary layer is introduced in the base profile for channel flow (see figure 1).We plot in figure 9 the neutral stability curve in the (W, β) plane for Couette flow with this Neumann boundary condition, and we see that PDI is suppressed in the small β limit.There is now a minimum β at which PDI exists, with PDI penetrating to lower β when Re is increased, showing more evidence that inertia destabilises the system.This means that in the concentrated limit of a polymeric fluid PDI does not exist for Neumann conditions on the polymer stress.There is also a maximum W at which the instability exists.PDI is therefore sensitive to a change in boundary conditions.However, the diffusion-free boundary conditions seem more likely to be physically relevant than the Neumann condition presented here but this, of course, is not certain.

Connection to Melt Fracture
Our work indicates that the polymer diffusive instability (PDI) is relevant for polymer melts where β is very small, due to the divergence of the growth rate of PDI as β → 0. PDI could therefore play a role in the melt fracture of extruded plastics.The 'sharkskin' instability is a short wavelength distortion seen on the surface of extruded plastics that is associated with melt fracture [29,32] and we now explore whether this could in fact be PDI.We begin by comparing the observed wavelengths and those of the PDI. Figure 13a in [32] shows the sharkskin wavelength found experimentally in a channel with half-width h is in the range, λ sharkskin ∈ (h/10, h).Using the non-dimensionalisation and notation introduced earlier, this corresponds to a wavenumber of k ∈ (2π, 20π).We will consider what value the polymer stress diffusion constant ε should take for inertialess PDI to be unstable for such k.Taking W ∼ 10, β ∼ 0.1, then, from figure 3, k ∼ 0.04 at the onset of instability.Since then ε ∈ (4×10 −5 , 4×10 −3 ), which is a reasonable value for the polymer diffusion to take [10].Hence the wavelength at which PDI is unstable in a plastic extrudate is comparable to the wavelength of the observed sharkskin instability.
For PDI to be able to affect the extruding process, the growth rate must be large enough that any perturbation at the start of the process is sufficiently amplified by the time the plastic has exited the die.PDI is active at a lengthscale εβ/W (1 − β) away from the boundary, and so the deforming material travels down the tube with velocity U ∼ εβ/W (1 − β).Hence the instability grows over timescale T ∼ L W (1 − β)/εβ where L is the length of the die in units of half-channel height.The growth rate of PDI (from figure 2b) is σ ∼ 0.001W (1 − β)/β, and so the amplification of a perturbation over time T is exp Taking the parameter values, ε ∼ 10 −3 , W ∼ 5, β ∼ 0.1 and L ∼ 10, yields a huge amplification of 3 × 10 41 .However, while clearly strong enough, PDI does not reproduce the relationship between wavelength and shear rate that is seen for the sharkskin instability.Specifically, the wavenumber of sharkskin is seen to decrease with shear (see figure 13a in [32] again), while the most unstable wavenumber associated with PDI increases with shear (see figure 10).

Physical or artificial?
We now discuss the possibility that PDI may be an artificial instability.This is an issue because the lengthscale of the boundary layer obtained for PDI approaches that of the polymer lengthscale, suggesting that the continuum assumption is not obviously valid [10].We consider the largest PDI lengthscale, which comes from the inertialess case (the most unstable wavenumber grows with inertia in figure 4a) and so this case is most likely to provide the required separation of scales for the continuum assumption to be reasonable.
Figure 4a suggests that the most unstable wavenumber for inertialess PDI occurs at k * ≈ 0.04, and hence dimensionally we have a most unstable wavelength of 150 εβ/W (1 − β)h.As in [10], the Stokes-Einstein relation is used to consider the ratio of the PDI boundary layer lengthscale to the polymer gyration radius R for parameters relevant to low density polyethylene (LDPE), which is known to show signs of the sharkskin instability [1].At temperatures of T ∼ 500K, LDPE has µ s ∼ 100P a • s, R ∼ 10 nm [33] and its relaxation spectra consists of a range of relaxation times as large as λ ∼ 100 s and as small as λ ∼ 10 −4 s [34].For W ∼ 10 and β ∼ 0.1 the PDI lengthscale is therefore where k B is the Boltzmann constant.This shows that the lengthscale of PDI is sensitive to the choice of λ.Large λ makes the lengthscale of PDI much larger than R, meaning the continuum assumption is valid as there is the required separation of scales.However, small λ makes PDI act on a lengthscale smaller than R, invalidating the continuum assumption.Due to this, it is unclear how physical PDI is for polymer melt systems.

Conclusion
The polymer diffusive instability is analytically tractable in the limit of β/W (1−β) → 0 for finite εRe/β.The eigenvalues of PDI can be predicted when boundary conditions of vanishing diffusion are used.We find β lim (Re) such that when β < β lim (Re), a prediction of the maximum growth rate is within 5% of its true value, and we see that β < 0.02 yields accurate predictions for Re < 1000.
In the limit of β/W (1 − β) → 0, the smallest W at which an instability is seen (W crit ) depends on the strength of inertia, measured by εRe/β.Inertia is shown to destabilise the system, with it allowing the system to become unstable at lower W than in the inertialess case.When εRe/β < 0.1, we obtain W crit ∼ 8, while when εRe/β > 5 we obtain W crit ∼ 2. Similar behaviour is even seen outside of the concentrated limit, when β ∼ 0.2.PDI is sensitive to a change in boundary conditions however, as changing the diffusion-free boundary conditions to vanishing stress flux boundary conditions changes where in parameter-space the instability is operative.In particular, this change suppresses PDI in the small β limit.
PDI seems particularly relevant to the case of polymer melts since the growth rate of PDI is so large.This indicates it could be relevant to real extrusion devices where the most unstable wavenumber is similar to the wavenumber of the observed 'sharkskin' instability for typical parameters.However, in some parameter ranges, the wavelength over which PDI operates approaches the lengthscale of the polymers themselves which invalidates the continuum viscoelastic model used (here Oldroyd-B but also FENE-P as discussed in [10]).We are therefore left to contemplate whether a 'better' continuum model would suppress PDI in these troublesome parts of parameter space where PDI may not be physical.Hopefully experiments can help resolve this as well as studies using more sophisticated continuum models.We hope to report on the latter (at least) in the near future.

Figure 1 :
Figure 1: The channel flow base state showing a) Txx b) Txy and c) U with boundary conditions of no diffusion (blue) and no stress-flux (orange), plotted alongside the base flow obtained when ε = 0 (black dotted).Parameters are W = 200, ε = 0.001, β = 0.05 and Re = 0.This shows how the presence of polymer stress diffusion can modify the base flow in channel flow.

Figure 6 :
Figure 6: a) Neutral stability curves in the (W , k) plane with changing Re with ε = 0.0001 and β = 0.01.These set the transition from non-inertial to inertial behaviour to be at Re * = 100.Darkest to lightest shade shows numerics for Re = 0, 100, 200, 400, 800, 1200, 1300, 1400.Dotted lines denote the neutral curve predicted by the analytic theory.Note when inertia is considered the prediction is less accurate -we expect that if a smaller β were considered that the prediction would match the numerics more closely.b) The relationship between W crit (the smallest W at which the system goes unstable) and inertia.Each dashed curve tracks changing Re with constant ε = 0.0001 and (darkest to lightest) constant β = 0.001, 0.01, 0.1, 0.2.The transparent grey curve shows the analytic prediciton.

Figure 7 :
Figure 7: The neutral curve in the (W , β) plane for plane Couette flow with diffusion-free boundary conditions at ε = 0.0001 and Re = 0, 500, 2000.The right figure zooms in on the small β region of the left figure, and uses a log scale.Solid lines correspond to numeric solutions, dotted lines correspond to the predictions based on the analytic results.Note each curve has constant Re, ε so as β → 0 we do not have εRe/β finite, which was required for the analytic work to be valid.This explains why the dotted prediction curves do not match the solid numeric curves exactly in the small β limit.These plots demonstrate that PDI is generic across β, and that inertia destabilises the system.

Figure 8 :
Figure 8: a) The relative error in the predicted maximum growth rate σ * pred as β varies for ε = 0, Re = 0 and (darkest to lightest) W = 10, 12, 20, 100, 200, 400.For β < β lim (Re = 0) = 0.090 (vertical dashed line), all predictions are within 5% of the true value (horizontal dashed lines).b) β lim as a function of Re.These plots verify that our predictions are asymptotically accurate, and demonstrate the size of β at which the predictions become valid.