Technical Note (draft) on the Conformation- and Concentration-Dependent Drag (C2D2) model for unentangled solutions of flexible polymers and numerical algorithms for rheological properties in steady and unsteady flows.

Technical Note (draft) on the Conformationand1 Concentration-Dependent Drag (C2D2) model for unentangled2 solutions of flexible polymers and numerical algorithms for3 rheological properties in steady and unsteady flows.4 Ranganathan Prabhakar5 (Dated: December 22, 2025)6 1 Abstract Status and purpose. This document is a technical report intended to support understanding, use, and development of the accompanying codebase. It is a working draft: notation, derivations, algorithmic details, and default numerical settings may evolve as the software matures. Portions of the model description and primary equations may later be adapted into one or more journal manuscripts; those manuscripts will cite a versioned, archived copy of this report and the corresponding code release. 7 I. INTRODUCTION8 The rheology of polymer solutions reflects the interplay of intramolecular and intermolec-9 ular hydrodynamic interactions (HI), finite extensibility (FE), and concentration effects.10 While excluded-volume and entanglement interactions are also important, HI are present11 in all regimes and therefore provide a useful starting point for theory. Classical treatments12 such as the Zimm model, which averages HI at equilibrium, established the importance of HI13 in dilute solutions and explained universal scaling in linear viscoelastic properties. However,14 the equilibrium-averaging assumption breaks down in flows with Wi &1, where single-chain15 Brownian dynamics (BD) simulations with HI show strong conformation dependence.16 To address this, several closure approximations have been developed. Consistent-17 averaging [1; 2], diagonalization approaches [3], and Gaussian approximations [4] extend18 beyond Zimm and successfully capture systematic trends in BD simulations of single chains.19 When combined with finite extensibility via the FENE-P closure, these models reproduce20 many qualitative features observed in unsteady and steady rheology. However, such closures21 exist only for the dilute, single-chain limit. This creates a fundamental disconnect: the very22 observation of strong viscoelastic effects in a complex flow already implies that polymer23 stresses are comparable to viscous stresses, and that molecules are sufficiently stretched24 to interact hydrodynamically with one another. Thus, experiments necessarily probe con-25 ditions beyond the range of validity of dilute single-chain models, while existing models26 cannot predict the regimes that experiments explore. At present, there are no tractable27 multimode models that capture the key physics of unentangled solutions at concentrations28 c/c?=O(1), where hydrodynamic screening and spectrum evolution become essential.29 2 FIG. 1. Concentration–stretch state space for the C2D2 model, showing regimes dominated by tension blobs (dilute and weakly-screened semidilute solutions), and correlation blobs (stronglyscreened semidilute solutions): Along the horizontal axis at stretch-ratio E/E0= 1 are equilibrium coil states with varying overlap-ratios (c/c∗); along the the vertical axis at c/c∗= 0 are isolated, partially-stretched chains. (a) The Zimm state of equilibrium ideal, isotropic, isolated random coils. (b) Isotropic equilibrium coils at critical overlap. (c) Semidilute polymer solution with interpenetrating isotropic equilibrium coils where chains are random walks of correlation blobs with intramolecular HI fully screened between those blobs. (d) Concentrated polymer solution at equilibrium with fully-screened intramolecular HI and fully Rouse chains. (e) Partially-stretched chain at infinite dilution consisting of a nearly-linear array (blob-pole) of tension blobs with weak Batechelor intramolecular HI. (f) Tension-blob poles just beginning to weakly interact hydrodynamically. (g) Tension-blob poles with weakly-screened intramolecular HI; the blue circles represent the Mackaplow-Shaqfeh screening length (h) Correlation blobs begin to dominate as core blobs. (i) Partially-stretched fully-Rouse chains in a concentrated polymer solution. (j) Fully-stretched isolated chain with weak Batechelor intramolecular HI. (k) Fully-stretched fully-Rouse chain in a concentrated polymer solution. 3 A natural framework for describing this crossover is the blob picture [5]. At concentra-30 tions above overlap (c>c?), coils can be described as strings of correlation blobs, within31 which intramolecular HI are retained but screened beyond the blob size by neighbouring32 chains. Under flow, stretched chains are described as arrays of tension blobs, below which33 thermal fluctuations dominate and above which anisotropy becomes evident. In both cases,34 HI act within blobs but are weakened between them, leading to a crossover from Zimm-35 like to Rouse-like dynamics as blobs shrink with concentration or stretch. The concentra-36 tion–stretch state space is sketched in Fig. 1. It highlights the dilute regime (tension blobs37 only), the intermediate semidilute regime (weakly screened HI), and the strongly screened38 regime dominated by correlation blobs. A tractable multimode constitutive model must39 incorporate this physics in order to describe polymer solutions relevant to experiments and40 CFD simulations.41 Recently, we proposed the single-mode conformationand concentration-dependent drag42 (C2D2) model, which incorporates HI screening through blob arguments and showed quali-43 tative agreement with experimental trends in uniaxial extension and capillary-breakup flows.44 The present work extends this framework in two directions. First, we formulate a multi-45 mode C2D2 model that connects microscale parameters (e.g. number of Kuhn segments,46 c/c?) to macroscale rheology, while allowing the relaxation spectrum to evolve continuously47 from Zimm-like to Rouse-like as chains stretch or overlap. Second, we introduce a variant of48 the FENE-PM model, denoted F-PME, that simplifies analysis and enables semi-analytical49 solutions for steady-state properties. The classical FENE–P dumbbell and its multimode50 version have long been used to describe the rheology of dilute polymer solutions. However,51 as emphasized by Wedgewood et al. [6] in their development of the FENE–PM chain model,52 the FENE–P description is limited to small number of modes since the calculations scale as53 the square of the number of modes, Nm. The FENE–PM chain of Wedgewood et al. reduces54 this computational cost by mode-specific terms by a mean value that allows a normal-modes55 transformation and diagonalization. This in turn leads to O(Nm) computation and also en-56 ables analytical treatment in steady shear and elongational flows. The FENE-PME variant57 proposed here uses the same reasoning to achieve O(Nm) but uses the chain extension Eper58 mode instead of the mean value in the FENE-PM model.59 Our aim is to make available a fast, tractable constitutive model that retains the essential60 physics of HI and finite extensibility, and that is suitable for CFD simulations of complex61 4 flows of polymer solutions. The focus of this paper is therefore on model formulation,62 analysis, and computational strategies, with sample results presented for linear viscoelastic63 spectra, steady, and unsteady uniaxial extensional flows. We obtain results for steady-state64 rheological properties in uniaxial extensional flows using the pseudo-arclength continuation65 (PAC) method. This is particularly useful for the coil-stretch hysteresis regime in extensional66 flows. Results for unsteady extensional are obtained using an exponential time-differencing67 (ETD) scheme. More detailed results for extensional and shear flows, including comparison68 against Brownian Dynamics simulations or experiments will be taken up in subsequent69 works. We also provide open-source code alongside this work so that interested researchers70 can carry out their own systematic validation.71 The remainder of this paper is organized as follows. Section II introduces the multimode72 C2D2 formulation, the FENE-PME variant, and the blob-based arguments for the dynamic73 relaxation time spectrum. Section III reports predictions for linear viscoelastic spectra,74 steady and unsteady extensional rheology. Finally, Section ?? summarizes the main findings75 and outlines directions for future comparison with experiments and applications in CFD.76 The Supplementary Information (SI) provides further details of the blob-model and the77 analysis underlying numerical computations of steady-state and unsteady properties.78 II. MODEL EQUATIONS79 A. Modal tensors and evolution equations80 We consider an ideal (phantom) chain at a given concentration, represented as a sequence81 of Nscoarse-grained subchains or entropic springs. We shall use the subscript ‘0’ to refer to82 the quiescent (thermodynamic) equilibrium state at a given polymer concentration, whereas83 the subscript ∞indicates the fully-stretched state. The instantaneous mean end-to-end84 distance is E, with E0its equilibrium value and E∞the contour length. We take E0to85 represent the coil size at equilibrium, corresponding to a volume E3 0. Thus, chains begin86 to overlap once their number density exceeds the critical overlap concentration c∗=E−3 0.87 Under Θ-conditions, the chain has Nk=E2 ∞/E2 0Kuhn segments of length bk=E2 0/E∞. A88 typical coarse-grained subchain contains Nk/NsKuhn segments; at equilibrium its mean-89 squared size is Q2 0=E2 0/Ns=b2 kNk/Ns, and its contour length is Q∞=E∞/Ns=bkNk/Ns.90 5 Any subchain possesses a spectrum of dynamical responses: the slowest relaxations are those that are coupled to motions of the chain as a whole, whereas the fastest are dominated by motions internal to the subchain itself. Motivated by this physical picture and by the structure of model equations obtained with standard closure approximations for bead–spring models incorporating HI, we assume that the relevant dynamics can be captured by a discrete set of Nmmodes and, for definiteness, take Nm=Ns. We introduce tensors Mp(p= 1, . . . , Nm) that track the modal contributions of the conformations of the subchains. These modal tensors are such that (i) at the quiescent (thermodynamic) equilibrium state, Mp,0= (Q2 0/3) Ifor all p, and (ii) the chain conformation tensor describing the average shape of a polymer chain is A= Nm X p=1 SpMp.(1) The trace of this tensor is interpreted as the mean-squared end-to-end distance: E2= Nm X p=1 Sptr Mp.(2) The end–to–end weight of mode pis determined by the orthogonal transform Π that maps spring–connector coordinates to normal modes: Sp=Nm X i=1 Πip2, p = 1, . . . , Nm.(3) For linear chains with free ends, Sp= 0 for even pand is nonzero only for odd p, with PNm p=1 Sp=Nm. Using tr Mp, 0=Q2 0then gives E2 0=NmQ2 0. As discussed later, intramolecular HI leads to Zimm–like spectra for isolated coils which change toward Rouse–like spectra when HI is weakened by stretch or intermolecular screening. In principle, the orthogonal matrix could vary across this crossover; in practice, Zimm and Rouse choices yield very similar Sp∼p−2for the slow odd modes and differ only near p∼O(Nm). For simplicity, we therefore adopt the Rouse weights throughout: for odd p, Sp=2 Nm+ 1 cot2ϑp, ϑp=pπ 2(Nm+ 1),(4) and Sp= 0 for even p.91 The modal tensors evolve under flow, entropic elasticity, finite extensibility (FE), and a HI-controlled relaxation spectrum. Their dynamics are governed by dMp dt=κ·Mp+Mp·κ>−1 λpfpMp−Q2 0 3I, p = 1, . . . , Nm,(5) 6 where κis the transpose of the velocity-gradient, Iis the identity, and λpdepends on chain92 size and concentration (specified below). The FE function fpis defined so that fp= 1 at93 equilibrium; hence for κ=0the steady state is Mp,0= (Q2 0/3) I. As chains approach full94 extension, fp→ ∞. We consider three FE closures.95 (i) FENE–P: In this classical closure, each mode stiffens according to its own trace: fp=1−Q2 0 Q2 ∞1−tr Mp Q2 ∞−1 .(6) (ii) FENE–PM: In this closure introduced by Wedgewood et al. [6], all modes share a common FE factor based on the mode-averaged conformation M: fp=f=1−Q2 0 Q2 ∞1−tr M Q2 ∞−1 ,M=1 Nm Nm X p=1 Mp.(7) (iii) FENE–PME: We introduce a new variant wherein the common FE factor is based on the mode-averaged mean-squared end-to-end distance (Eq. (2)): fp=f=1−Q2 0 Q2 ∞1−E2/Nm Q2 ∞−1 =E2 ∞−E2 0 E2 ∞−E2,(8) where, in this variant, we define Q2 ∞=E2 ∞/Nm. This choice simplifies analysis and96 enables efficient numerical methods while retaining all the required qualitative features97 of the predictions of either the FENE-P or FENE-PM models.98 The connection to macroscopic rheology is made by relating the modal tensors to the polymeric contribution to the total stress. We adopt a mean-field picture in which each chain evolves in an effective medium created by all other chains. In this setting, the polymeric stress (subscript ‘p’ denotes the polymer contribution, not the mode index) is given by the Kramers expression τp=−c H Nm X p=1fpMp−Q2 0 3I,(9) where H= 3 kbT/Q2 0. This ensures τp=0at equilibrium since fp= 1 and Mp,0= (Q2 0/3) I.99 Interchain interactions enter only through the concentration-dependence of the relaxation100 times. If this concentration dependence is ignored, then the polymeric stress contribution is101 linear in concentration. In contrast, a concentration-dependent relaxation spectrum makes102 rheological properties nonlinear functions of c. We next turn to the question of how chain103 stretching and polymer concentration determine the relaxation time spectrum.104 7 B. Blobs and the relaxation spectrum105 While the actual chain has Nkmodes down to the Kuhn scale, coarse-graining each chain into Nssubchains resolves only the Nm=Nsslowest modes; the remaining fast modes are assumed relaxed and contribute only to equilibrium properties. We denote the largest relaxation time (mode 1) by λ≡λ1, dropping the subscript when unambiguous. Its equilibrium value at a given concentration is λ0. In the dilute limit c→0, λ0→λz(the Zimm time). In the C2D2 framework, λdepends on the instantaneous chain conformation, so λ6=λ0. We quantify the stretch effect (at fixed c) through Λ ≡λ/λ0. To describe the full spectrum we also define the relative times ˆ λp≡λp/λ (with ˆ λ1≡1), so that λp=λ0Λˆ λp, p = 1, . . . , Nm.(10) Below, we first quantify how the chain-averaged friction coefficient ζdepends on the stretch106 ratio E/E0and the overlap ratio c/c∗, and how this determines λ0(c/c∗) and Λ(E/E0, c/c∗).107 We then characterise how the relative spectrum {ˆ λp}varies with stretch and overlap.108 1. Chain friction and the longest relaxation time109 The dynamic relaxation spectrum is governed by hydrodynamic interactions (HI) be-110 tween chain segments, which in turn depend on the instantaneous stretch E/E0and the111 (equilibrium) coil-overlap c/c∗. As outlined in the Introduction, C2D2 uses blob concepts to112 compute the average chain friction at any instant. The tension blob is the scale below which113 local structure is an isotropic random walk; above it, stretch-induced anisotropy is apparent.114 The correlation blob is the scale within which segments of a chain interact hydrodynami-115 cally with each other unhindered by segments from other chains. Together, E/E0and c/c∗ 116 determine the size ξband number Nbof core blobs—the scale at which Zimm dynamics117 still prevails locally within ideal, isolated, isotropic random-walk domains containing Nk/Nb 118 segments (Fig. 1; SI). In dilute or weak-screening regimes, tension blobs are the core blobs;119 in strong-screening regimes, correlation blobs are the core blobs.120 The chain-averaged friction coefficient is written as ζ=ζbNb G,(11) 8 where Nbζbis the Rouse drag of Nbcore blobs of size ξb, and Gis a Batchelor-type cor-121 rection that accounts for weak residual intramolecular HI between nearly aligned blobs (a122 “blob pole”) [7; 8]. We have G= 1 at equilibrium and in the strong-screening regime; for123 partially stretched chains in dilute or weak-screening regimes, G=G(E/E0, c/c∗) with an124 interpolation summarised in the SI.125 Within an isotropic core blob, the strength of unscreened intramolecular HI is set by the Kuhn-scale parameter h∗ k≡ak √π bk ,(12) where the segment hydrodynamic radius akrelates to the Stokes friction of a single Kuhn segment, ζk= 6π ηsak, in a solvent of viscosity ηs. Given h∗ k, the Kirkwood–Riseman (KR) procedure (SI) yields the mean friction of any ideal, isolated, isotropic random-walk segment containing nKuhn steps via a draining function α(h∗ k, n): ζb=α(h∗ k, Nk/Nb)ηsξb, ξb=bkrNk Nb .(13) Thus αquantifies how the drag of a random-walk coil differs from the Stokes drag of an126 impenetrable sphere of the same size. In the non-draining limit of infinitely large random127 walks, α→αnd ≈5.8 (independent of h∗ k); finite coils with h∗ k.0.25 are partially draining,128 with α < αnd but increasing toward αnd as ngrows at fixed h∗ k.129 At equilibrium, the chain friction at any given c/c∗follows from (11) with G= 1. For c/c∗≤1 the whole isotropic coil is a single blob; for c/c∗>1 the coil is a Rouse chain of Nc correlation blobs (each of friction ζc). Hence ζ0=       ζz, c/c∗≤1, ζcNc, c/c∗>1, (14) with the Zimm whole-coil drag (dilute, equilibrium) ζz=α(h∗ k, Nk)ηsE0,(15) and, for c/c∗>1, ζcgiven by (13) with Nb=Nc= (c/c∗)2and ξb=ξc=bkpNk/Nc.130 We now evaluate the equilibrium longest relaxation time λ0and the ratio Λ needed in (10). In the dilute limit, the longest time is the Zimm time. In the non-draining limit, Zimm 9 [7]: P1(N) = 0.91 −0.02 √N+0.10 N,(39) P2(N) = 1.6−1.1 √N−0.5 N.(40) Thus, for a given h∗ kand N,α(h∗ k, N) can be evaluated directly.189 As N→ ∞,αapproaches the universal non-draining limit αnd ≃5.8, independent of h∗ k.190 When N= 1, α(h∗ k,1) = 6π3/2h∗ k, recovering the drag of a single Kuhn segment. Figure 3191 shows the smooth interpolation between these limits.192 α 1 10 1/ N 0.01 0.1 1 0.5 h* K = 0.02 0.03 0.2 FIG. 3. Draining function α(h∗ k, N) as a function of 1/N for different h∗ k. Symbols: Kirkwood–Riseman estimates. Curves: polynomial fits. Circle: universal non-draining value αnd ≈5.8. SI-2. AVERAGE FRICTION COEFFICIENT AND INTERPOLATION SCHEME193 FOR G194 Given the stretch ratio E/E0≥1 and the concentration ratio c/c∗, the number and size of the tension or correlation blobs are obtained as Nt= (E/E0)2, ξt/E0= (E/E0)−1(41) Nc= (c/c∗)2, ξc/E0= (c/c∗)−1,if c/c∗≥E/E0,(42) 16 where Ntand ξtare the number and isotropic size of tension blobs, and Ncand ξcare the195 corresponding quantities for correlation blobs.196 There are three principal regions in this state space. The dilute regime is defined by the197 regime where c/c∗<(E/E0)−3<1 where only t-blobs are present. In the intermediate198 semidilute regime where (E/E0)−3< c/c∗< E/E0, t-blobs are the core blobs and hydro-199 dynamic screening is weak. C-blobs are the core blobs in the semidilute regime with strong200 hydrodynamic screening; this occurs when c/c∗> E/E0.201 In the weak-screening regime, the hydrodynamic screening length is not ξcbut the value χgiven by Mackaplow and Shaqfeh [16] for a suspension of aligned rods of length Eand diameter ξtsuch that χ ξt =1 φt ln 1 φt1/2 ,(43) where the volume fraction of t-blob poles, φt=c E ξ2 t=c/c∗ E/E0 .(44) The expected scaling results for the drag coefficient are:202 •Dilute coiled chains (c/c?1, E/E0= 1): ζ=ζz.203 •Dilute stretched chains (c/c?(E/E0)−3): ζ≃ζtNt/ln(2Nt).204 •Weakly screened semidilute ((E/E0)−3c/c?E/E0): ζ≃ζtNt/ln(2χ/ξt).205 •Strongly screened semidilute (c/c?E/E0): ζ=ζcNc.206 •Concentrated (c/c?≥√Nk): ζ=ζR.207 To interpolate smoothly between these scaling results, once the regime and the relevant core blob is identified, the effective chain friction in the C2D2 model is expressed as ζ=ζbNb G, ζb=α(h∗ k, Nk/Nb)ηsξb, ξb=bkpNk/Nb,(45) with Nb=Ntor Nc, and ζb=ζtor ζcdepending on the regime, where ζtand ζcdenote208 the blob frictions in the tensionand correlation-blob regimes, respectively. The draining209 function αhas been discussed above in SI-1.210 17 The factor Gis the Batchelor-type correction interpolating across dilute, weakly screened, and strongly screened regimes. At equilibrium (when E/E0= 1) and in the strong-screening regime (E/E0> c/c? > 1), G= 1. In the dilute and in the weak =-screening regimes where where tension blobs are the core blobs, G=G1 G2 ,(46) where G1= 1 + ln A+ (B−1)J A+B,(47) G2= 1 + ln B (1/A) + B,(48) A=R/ξt, B = 1 + χ/ξt,(49) J= 1 + (1/2)(A−1) + (2/e)(A−1)2 A.(50) This form ensures recovery of the correct limits at the weak/strong screening boundary211 and avoids discontinuities in G. The interpolation scheme above unfortunately produces212 shallow unphysical local minima when E/E0is close to 1 or near the weak/strong screening213 boundary. To remove these, we could calculate first ζwith the interpolation formula above,214 and then choose max (ζ, ζ0) as the final value of the friction coefficient. This approach will215 however eliminate the possibility of the friction coefficient decreasing below ζ0as chains216 are stretched, which is physically possible for small values of h∗ k. Therefore, the unphysical217 minima can be left in place while being careful about interpreting any consequential feature218 of those minima in predictions.219 SI-3. PSEUDO-ARCLENGTH CONTINUATION FOR MULTI-MODE C2D2220 MODELS AT STEADY STATE221 The dimensionless modal tensors evolve according to222 ˙ Mp= Wi (eκ·Mp+Mp·eκ >)−θp(Z)Mp+σp(Z)I,(51) where Zis a vector of a set of auxiliary scalar variables and θp(Z) = fp(Z) λp(Z)=fp(Z) Λ(Z)ˆ λp(Z), σp(Z) = 1 3λp(Z)=1 3 Λ(Z)ˆ λp(Z).(52) 18 The set of Na+ 1 auxiliary variables, Z={X1, X2, . . . , XNa, Y }, where Y=E2satisfy a223 system of Na+ 1 closures,224 Fj(Z,{T(Mp)})=0, j = 1, . . . , Na+ 1.(53) While the closure functions may more generally involve other invariants of Mp, for the225 present we shall assume they only depend on the trace functions of the modal tensors. We226 require that the last closure equation is227 FNa+1 = Nm X p=1 SpT(Mp)−Y= 0 .(54) The three FENE variants are distinguished by their auxiliary closure equations.228 a. FENE–PME There are no auxiliary variables other than Y=E2, that is, Na= 0.229 F1= Nm X p=1 SpT(Mp)−Y= 0. b. FENE–PM The auxiliary variables are Z= (X, Y ) which satisfy: F1=1 Nm Nm X p=1 T(Mp)−X= 0,(55) F2= Nm X p=1 SpT(Mp)−Y= 0.(56) c. FENE–P The auxiliary variables Z= (X1, . . . , XNm, Y ) satisfy Fp=T(Mp)−Xp= 0, p = 1, . . . , Nm,(57) FNm+1 = Nm X p=1 SpT(Mp)−Y= 0.(58) Wedgewood et al. [6] showed that the FENE–PM chain equations can be algebraically230 manipulated to yield analytical or semi–analytical steady–state results in simple shear and231 elongational flows. Here we take a similar algebraic approach to the multi-mode C2D2232 framework above, organizing the analysis around admissibility and pole singularities that233 govern the steady–state manifold in incompressible extensional flows and planar mixed flows.234 This structural viewpoint applies uniformly across the entire multi-mode C2D2, rather than235 to a single model. The central idea is that, instead of Wi, the steady–state solution is236 parameterized by Y≡E2; Wi is obtained as a function of E2instead. A pseudo–arclength237 19 continuation (PAC) method is used to obtain the steady–state solution parameterized by Y238 (with the Newton’s method being used as a fallback). Before implementing this numerical239 scheme, the steady states in extensional flows and planar mixed flows are analyzed to firstly240 show that, for any given Y, a unique Wi exists that ensures that the components of the241 modal conformation tensors are non–negative. This analysis provides admissibility criteria242 that need to be checked as the steady–state solutions are obtained with the PAC method.243 The analysis also provides a fast approximation for Wi and the rest of the variables under244 certain conditions.245 At steady state, the p–th conformation mode satisfies246 0= WieκMp+MpeκT−θp(Z)Mp+σp(Z)I, p = 1, . . . , Nm,(59) and the auxiliary variables Z≡(X1, . . . , XNa, Y ) are obtained as the solution of a system of closure equations for the steady state: Fj(Z) = 0, j = 1, . . . , Na,(60a) FNa+1(Z,Wi) ≡ Nm X p=1 SpTp(Z,Wi) −Y= 0, Sp>0, Nm X p=1 Sp= 1.(60b) Here Tpare the modal trace functions obtained after rearranging the steady–state equa-247 tions to eliminate the modal components and express them explicitly in terms of Zand248 Wi. The specific forms of these closure systems for FENE–PME, FENE–PM and FENE–P249 were introduced above. At steady state, they reduce respectively to one, two, and Nm+ 1250 coupled algebraic equations. We analyze their common structure in what follows, focusing251 on admissibility and pole singularities.252 1. Pseudo–arclength continuation253 The steady auxiliary–variable closure for a homogeneous flow defines a nonlinear manifold in an augmented state space. Let Z∈RNa+1 collect the auxiliary closure variables (e.g. Z= (X, Y ) for FENE–PM and FENE–P, and Z= (Y) for FENE–PME), and let ηdenote a scalar continuation coordinate. Depending on the flow and admissibility structure, ηmay be taken as Wi or as a reparameterisation of Wi (Sec. 3). The steady closure can be written abstractly as F(Z, η) = 0,F:RNa+2 →RNa+1,(61) 20 and hence defines a one–dimensional solution manifold in (Z, η)–space.254 Pseudo–arclength continuation (PAC) advances along this manifold using a predictor– corrector scheme. Let ˜ z= (Z, η)∈RNa+2. Given a converged point ˜ zkand an estimate of the local tangent ˆ tk, a predictor is formed as ˜ zk+1 pred =˜ zk+ ∆skˆ tk,(62) where ∆skis the continuation step size. The corrector then solves the augmented system ˜ G(˜ z) =  F(˜ z) g(˜ z) =0,(63) where gis a scalar constraint that fixes the arclength phase. A standard choice is the hyperplane condition g(˜ z) = ˆ tk T ˜ z−˜ zk+1 pred= 0.(64) Any smooth, transverse scalar constraint yields a valid PAC formulation. In practice, when255 Zis high–dimensional (e.g. FENE–P with Z∈RNm+1), it is advantageous to impose the256 arclength constraint in a low–dimensional physically interpretable subspace, while enforcing257 the full closure F=0; we return to this implementation choice in Sec. 3.258 The tangent ˆ tkmay be obtained from the previous converged points (secant approx-259 imation), or by solving the bordered linear system associated with (63). The nonlinear260 corrector (63) is solved using Newton iterations with globalization (Sec. 3), and ∆skis261 adjusted adaptively to maintain robustness and efficiency.262 2. Pole structure and admissibility263 Consider a vorticity–free incompressible extensional flow with Wi >0. In the principal264 frame of the flow,265 eκ= diag(α1, α2, α3), α1+α2+α3= 0, αmax ≡max iαi>0, Mpis diagonal at steady state. Defining θp>0 and σp>0 as before (Eqn. (52)),266 θp−2Wi αiMp, ii =σp, i = 1,2,3. Hence,267 Mp, ii =σp θp−2Wi αi , Tp(Wi,Z) = σp 3 X i=1 1 θp−2Wi αi .(65) 21 a. Poles and admissibility. If Zis held fixed and Wi is varied, simple poles occur when268 θp−2 Wi αi= 0 i.e. when Wi = θp 2αi . Physical admissibility requires Mp, ii >0 for all i. Since σp>0 and θp>0, this implies269 θp>2 Wi αmax for each mode p. (66) Therefore a global admissible interval for Wi is270 0<Wi <θmin(Z) 2αmax , θmin(Z) := min 1≤p≤Nm θp(Z).(67) For any state on the state-state manifold, the condition above gives a global upper bound271 on Wi for which the multimode system has admissible solutions.272 b. Uniqueness of physically admissible Wi.Given an admissible Z, we shall show in273 subsequent paragraphs the following:274 1. For physically realistic closures in which the finite–extensibility factors are nondecreas-275 ing with Yalong the closure manifold (i.e. ∂fp/∂Y ≥0), the weighted sum in the last276 auxiliary closure function277 FNa+1(Wi; Z) = Nm X p=1 SpTp(Wi; Z)−Y is negative at Wi = 0, i.e. FNa+1(0; Z)<0 whenever Y > Y0.278 2. The trace function Tp(Wi,Z) is strictly increasing in Wi on the admissible interval (67).279 The second point implies that FNa+1(Wi; Z) is strictly increasing in Wi on (67). Moreover,280 as Wi ↑θmin(Z)/(2αmax) within this interval, at least one denominator in (65) vanishes,281 hence PpSpTp(Wi,Z)→+∞and FNa+1 →+∞. By the intermediate value theorem, there282 is exactly one zero crossing of FNa+1 on the admissible interval. Therefore, for each fixed283 admissible Z, there exists a unique physically admissible Wi.284 c. Sign of FNa+1 at Wi = 0.From (65), Tp(0,Z)=3σp/θp= 1/fp(Z), hence285 FNa+1(0,Z) = Nm X p=1 Sp 1 fp(Z)−Y. Let fp,0:= fp(Z0) denote equilibrium values, so Tp(0,Z0)=1/fp,0and therefore286 FNa+1(0,Z0) = Nm X p=1 Sp 1 fp,0−Y0= 0. 22 Subtracting this equilibrium identity from FNa+1(0,Z) yields287 FNa+1(0,Z) = Nm X p=1 Sp1 fp(Z)−1 fp,0−(Y−Y0). Along the steady–state manifold out of the quiescent equilibrium in extensional flows, Y >288 Y0. Further, for physically realistic models, as chains stretch we have fp(Z)≥fp,0for all p289 (equivalently ∂fp/∂Y ≥0). Hence each difference 1/fp(Z)−1/fp,0≤0, so FNa+1(0,Z)≤290 Y0−Y < 0.291 d. Monotonicity of Tpin Wi at fixed Z.For any mode p, abbreviate Tp(Z,Wi) ≡T,292 θ≡θp(Z)>0, σ≡σp(Z)>0. Differentiating (65) with respect to Wi (at fixed Z) gives293 ∂T ∂Wi =σ 3 X i=1 2αi θ−2 Wi αi2= 2σ H(θ, Wi; α), H(θ, Wi; α) := 3 X i=1 αi (θ−2Wi αi)2. We now show H > 0 on the admissible interval θ > 2Wi αmax for both uniaxial and biaxial294 extensional flows. Planar extensional flows are a degenerate case of either kind of extensional295 flow (obtained with a= 0 or 1 in the expressions below).296 Uniaxial extension: By appropriately defining Wi, the shape tensor of any uniaxial exten-297 sional flow can be parameterized as α1=αmax = 1, α2=−a,α3=a−1, a∈[0,1]. Then,298 defining299 h(x) := x (θ+ 2Wi x)2, for any a,300 H(θ, Wi, a) = 1 (θ−2Wi)2−h(a)−h(1 −a), In the admissible regime, θ > 2Wi αmax = 2Wi, and h(a) is increasing and concave on301 a∈[0,1]. Then, by Jensen’s inequality302 h(a) + h(1 −a) 2≤h1 2=1 (θ+ Wi)2. Therefore303 H(θ, Wi, a)≥1 (θ−2Wi)2−1 (θ+ Wi)2=3 Wi (2θ−Wi) (θ−2Wi)2(θ+ Wi)2>0. Biaxial extension: Write α1=−1, α2=a,α3=αmax = 1 −a, with a∈[0,1/2] so304 αmax ∈[1/2,1]. Redefine305 h(y) := y (θ−2Wiy)2, 23 so that306 H(θ, Wi, a) = −1 (θ+ 2Wi)2+h(a) + h(1 −a). On the admissible interval θ > 2Wi αmax = 2Wi(1 −a) we have, for all y∈[0, αmax],307 θ−2Wiy≥θ−2Wiαmax >0, so the derivatives308 h0(y) = θ+ 2Wiy (θ−2Wiy)3>0, h00(y) = 8 Wi(θ+ Wiy) (θ−2Wiy)4>0 show his increasing and convex on [0, αmax]. Since a, 1−a∈[0, αmax] and the midpoint is309 always 1 2, Jensen’s inequality yields310 h(a) + h(1 −a)≥2h1 2=1 (θ−Wi)2. Therefore, for all θ > 2Wi αmax,311 H(θ, Wi, a)≥1 (θ−Wi)2−1 (θ+ 2Wi)2=3 Wi (2θ+ Wi) (θ−Wi)2(θ+ 2Wi)2>0. Hence H > 0 on the full admissible interval θ > 2Wi αmax (not merely for θ > 2Wi), and312 thus ∂WiTp>0 there in biaxial extension.313 a. Planar mixed flows314 The transpose of the velocity gradient of a planar mixed flow can be written as315 κ= Wi bκ(ω),bκ(r) =      1r0 −r−1 0 0 0 0     , where the flow–type parameter ris the ratio of the vorticity to the strain rate, Wi =316 ˙γ. This representation makes the analysis below convenient. Other planar mixed–flow317 parameterizations can be mapped to the form above by a change of the flow-type variable318 and rescaling the strength.319 For mode p, with θp=fp/λpand σp= 1/(3λp) as before, the steady state solution for320 the modal tensors is of the form321 Mp=     ApCp0 CpBp0 0 0 Dp     . 24 Solving the steady-state balance (59) yields the closed–form components Ap=σpθ2 p+ 2θpWi + 4r2Wi2 θpθ2 p−4(1 −r2) Wi2,(68a) Bp=σpθ2 p−2θpWi + 4r2Wi2 θpθ2 p−4(1 −r2) Wi2,(68b) Cp=−4σprWi2 θpθ2 p−4(1 −r2) Wi2,(68c) Dp=σp θp .(68d) a. Poles and admissibility. From (68), all in–plane components share the common322 denominator323 θpθ2 p−4(1 −r2) Wi2. Since θp>0, the only potential singularity that arises as Wi is varied at fixed Zis when324 θ2 p−4(1 −r2) Wi2= 0 or when Wi = θp 2√1−r2(0 ≤r < 1).(69) Thus:325 •Strain–dominated mixed flows (0 ≤r < 1): admissibility requires326 0<Wi <θmin(Z) 2√1−r2, θmin(Z) := min 1≤p≤Nm θp(Z),(70) which guarantees finiteness of all in–plane components for every mode p.327 •Rotation–dominated cases (r≥1, including simple shear at r= 1): the factor in (69)328 is strictly positive for all Wi, so no pole occurs from the kinematics.329 To be physically valid, the modal tensors must be positive definite, even if the off–diagonal Cpis negative. First, Dp=σp/θp>0 since σp>0 and θp>0. Next, when (69) holds with strict inequality (i.e. in the admissible set), the common denominator is positive; the numerator of Apis a sum of nonnegative terms, and hence Ap>0. To understand the sign of Bp, we first check the in–plane 2 ×2 block by evaluating its determinant: ApBp−C2 p=σ2 p θ2 p (θ2 p+ 2θpWi + 4r2Wi2)(θ2 p−2θpWi + 4r2Wi2)−(4rWi2)2 θ2 p−4(1 −r2)Wi22 =σ2 p θ2 p (θ2 p+ 4r2Wi2)θ2 p−4(1 −r2)Wi2 θ2 p−4(1 −r2)Wi22 =σ2 p θ2 p θ2 p+ 4r2Wi2 θ2 p−4(1 −r2)Wi2>0. 25 where Znis obtained from {Mp(εn)}. With these constants, the strain–form evolution (??)442 reduces, over the step, to the linear inhomogeneous matrix ODE443 dMp dε =eκMp+Mpeκ >−θn p WinMp+σn p WinI.(80) a. Exact ETD1 update (frozen coefficients). Define the 3 ×3 step matrix444 An p:= eκ−θn p 2 WinI, so that (80) is dMp dε =An pMp+Mp(An p)>+σn p WinI.The exact solution over the strain step is445 Mn+1 p=ehAn pMn peh(An p) >+Zh 0 esAn pσn p WinIes(An p) >ds. (81) Let Xn pbe the (unique) solution of the 3 ×3 continuous Lyapunov equation446 An pXn p+Xn p(An p)>=−σn p WinI.(82) Then the integral in (81) has the identity447 Zh 0 esAQesA >ds =ehAXehA >−X for the pair (A,Q) satisfying AX +XA>=−Q, which yields the compact update448 Mn+1 p=ehAn pMn p+Xn peh(An p) >−Xn p.(83) This is our ETD1 scheme. It integrates exactly the frozen kinematic part and the stiff linear449 relaxation, while adding the isotropic fluctuation through the closed–form Lyapunov term.450 b. Evaluation of the Lyapunov solve. In the ETD1 update the auxiliary matrix Xn pis451 defined as the symmetric solution of the continuous Lyapunov equation452 An pXn p+Xn p(An p)>=−σn p WinI,An p=eκ−θn p 2WinI. Because An pis 3 ×3, this system admits closed forms in canonical flows:453 In the principal frame eκ= diag(α1, α2, α3), so An p= diag(λ1, λ2, λ3) with λi=αi−454 θn p/(2Win). The Lyapunov equation decouples componentwise, yielding455 Xii =−(σn p/Win) 2λi , Xij = 0, i 6=j. 32 c. Evaluation of the matrix exponentials. In (83) the update involves congruence with456 ehAn pand eh(An p) >. Since An p=eκ−(θn p/2Win)Iis 3 ×3, its exponential admits closed forms457 in canonical flows:458 In the principal frame,459 ehAn p= diageh(α1−θn p/2Win), eh(α2−θn p/2Win), eh(α3−θn p/2Win). Thus in the extensional flows, ETD1/ETD2 reduce to congruence operations with either460 explicit diagonal exponentials (extension) or simple 2 ×2 formulas. This is also easily461 extended to planar mixed flows.462 d. Step size h.At each stage the linear operator has the form463 An=eκ−θn p 2 WinI, so its eigenvalues are those of eκshifted by −θn p/(2Win). For uniaxial/biaxial extension,464 these are real.465 In ETD, ehAnis evaluated exactly, so strongly damped directions that normally cause466 stiffness place no stability restriction on h. Practical limitations on harise from (i) accuracy467 when tracking rapidly decaying modes, and (ii) the nonlinear variation of θp(Z), σp(Z) with468 the closure variables. Near admissibility boundaries (“pole–hugging” states), these coeffi-469 cients vary sharply with Y, so adaptivity must shrink hto resolve the nonlinear changes.470 [1] H. C. ¨ Ottinger. Generalized Zimm model for dilute polymer solutions under theta conditions.471 J. Chem. Phys., 86:3731–3749, 1987.472 [2] H. C. ¨ Ottinger. A model of dilute polymer solutions with hydrodynamic interaction and finite473 extensibility. I. basic equations and series expansions. J. Non-Newtonian Fluid Mech., 26:474 207—-246, 1987.475 [3] A. J. Kishbaugh and A. J. McHugh. A discussion of shear-thickening in bead-spring models.476 J Non-Newton. Fluid Mechanics, 34:181–206, 1990.477 [4] R. Prabhakar and J. R. Prakash. Gaussian approximation for finitely extensible bead-spring478 chains with hydrodynamic interaction. J. Rheol., 50:561–593, 2006.479 [5] M. Rubinstein and R. H. Colby. Polymer physics. Oxford University Press, London, UK,480 2003.481 33 [6] L. E. Wedgewood, D. N. Ostrov, and Bird R. B. A finitely extensible bead-spring chain model482 for dilute polymer solutions. J. Non-Newtonian Fluid Mech., 40:119–139, 1991.483 [7] R. Prabhakar, S. Gadkari, T. Gopesh, and M. J. Shaw. Influence of stretching induced self-484 concentration and self-dilution on coil-stretch hysteresis and capillary thinning of unentangled485 polymer solutions. J. Rheol., 60:345–366, 2016.486 [8] R. Prabhakar, C. Sasmal, D. A. Nguyen, T. Sridhar, and J. R. Prakash. Effect of stretching-487 induced changes in hydrodynamic screening on coil-stretch hysteresis of unentangled polymer488 solutions. Phys. Rev. Fluids, 2(1):011301, January 2017.489 [9] M. Doi and S. F. Edwards. The Theory of Polymer Dynamics. Oxford University Press, 1986.490 [10] G. G. Fuller and L. G. Leal. The effects of conformation-dependent friction and internal491 viscosity on the dynamics of the nonlinear dumbbell model for a dilute polymer solution. J.492 Non-Newtonian Fluid Mech., 8:271–310, 1981.493 [11] M N Hounkonnou, C Pierleoni, and J-P Ryckaert. Liquid chlorine in shear and elongational494 flows: A nonequilibrium molecular dynamics study. J. Chem. Phys., 97(12):9335–9344, De-495 cember 1992.496 [12] P. Debye and A.M. Bueche. Intrinsic viscosity, diffusion and sedimentation rate of polymers497 in solution. J. Chem. Phys., 16:573–579, 1948.498 [13] R. G Larson. The Structure and Rheology of Complex Fluids. Oxford University Press, 1999.499 [14] J G Kirkwood and J Riseman. The intrinsic viscosities and diffusion constants of flexible500 macromolecules in solution. J. Chem. Phys., 1948.501 [15] H. Yamakawa. Modern Theory of Polymer Solutions. Harper and Row, New York, 1971.502 [16] M. B. Mackaplow and E. S. G. Shaqfeh. A numerical study of the rheological properties of503 suspensions of rigid, non-Brownian fibres. J. Fluid Mech., 329:155–186, 1996.504 34