Capillary Thinning of Viscoelastic Threads of Unentangled Polymer Solutions

In this paper, we theoretically consider the process of the capillary thinning of a polymer fluid thread bridging two large immobile droplets in the regime of highly stretched polymer chains. We first derive a new relation between the pressure p and the flow velocity v in unentangled polymer solutions, which is called the anti-Bernoulli law: it shows that p is higher where v is faster. Using this equation, it is shown that the flow field is asymptotically irrotational, in particular, in the thread/droplet transition zones (in the case, the negligible solvent viscosity and inertial effects). On this basis, we predict the free surface profile and the thread thinning law for the FENE-P model of polymer dynamics. The predictions are compared with recent theoretical results and some experimental data on capillary thinning.


Introduction
Long liquid cylinders tend to break in many droplets due to capillary surface forces. For a Newtonian liquid, this phenomenon, related to the well-known Plateau-Rayleigh instability, typically results in a localized pinching, leading to the formation of necks whose thicknesses rapidly decrease over time. In contrast, the thinning of polymer solution threads often result in a transient formation of nearly uniform (perfectly cylindrical) threads. The characteristic thinning time of such threads strongly increases with polymer molecular weight. This phenomenon was observed in numerous experiments [1][2][3][4], in particular those using the liquid filament rheometer (LFR) [5]. Most experiments show that the diameter of a polymer thread is thinning according to an exponential law [1,6,7].
Theoretically, such uniform exponential thinning of polymer filaments was deduced from rheological (force balance and constitutive) equations for the Oldroyd-B and FENE-P models [1,7,8]. There was, however, a problem with the early theories: the boundary condition at the filament ends was stated based on a plausible argument rather than a rigorous proof. ('We assume that the axial stress vanished because, in the LFR, the filament is attached to large stagnant drops on stationary end plates', as written in ref. [8].) Thus, Entov and Hinch [8] proposed that the axial polymer stress σ p in the thread of radius a 0 must be equal to the capillary pressure, σ p = γ/a 0 (here, γ is the surface tension) to obtain the total axial pressure equal to the atmospheric pressure well inside the large drops at the thread ends. A similar equation, but with the prefactor 1/2 in the rhs, was obtained earlier [9] based on the slender-body approximation, which, however, is not applicable in the decisive interfacial zone between the thread and an end-droplet. Later, Clasen et al. [7] derived a different condition: applying the slender-body approximation in a different way (see also [10]). Quite recently, the above equation was re-derived [11] based on the Onsager variational principle (in the regime of negligible viscous stress due to the solvent). Equation (1) was finally rigorously established [12] for the capillary collapse of a neo-Hookean elastic cylinder in the limit of small elastic modulus G. It was argued [12,13] that the neo-Hookean elastic model becomes z z m L f L t L t Figure 1. A thinning filament of radius a 0 and length L f connecting two semi-spherical droplets (a cross-section along the main axis is shown here); z m corresponds to the middle of the thread. L t L f is the size of the filament/droplet transition regions (shown in cyan).
The filament is typically thin, long, and nearly uniform during the polymer extension process [6]: L f a 0 (t) (2) where 2a 0 (t) is the thickness in the middle of the thread (z = z m ); a(z, t) a 0 (t) in the filament region. As the liquid is incompressible, the uniaxial extension rate in the thread is defined by the thinning rate (note that˙ r = ∂(ln a 0 )/∂t in Equation (3) is the (negative) extension rate in the radial direction, and that (by virtue of the incompressibility)˙ + 2˙ r = 0): = ∂v z /∂z −2∂(ln a 0 )/∂t (3) so the axial flow velocity is: v z ˙ (z − z m ) in the filament region of length L f . Outside this region, the liquid cross-section increases rapidly, leading (by virtue of the incompressibility) to a decrease of v z . Thus, the velocity maximum, v 0 = max|v z | ˙ L f /2, is reached near the filament ends.
Let us turn to the force balance in the transition zones next to the filament ends where the surface shape deviates from a perfect cylinder and its radius shows a significant zdependence. The transition region between the thread and the second droplet is depicted in Figure 2. (Note that the origin, z = 0, is located near the righthand end of the thread. We define the point z = 0 by the condition a(z = 0) = k e a 0 , where the factor k e can be chosen as convenient; here it was set k e ≈ 1.25.) As argued below, the axial size, L t , of this region (where the droplet surface significantly deviates from both cylindrical and spherical) is relatively small, L t ∼ a 0 L f . Note that a longitudinal cross-section including the cylindrical symmetry axis (z) is shown here. Thus, the upper curve corresponds to the cylindrical angle θ = 0, the lower curve to θ = 180 • , and the r-axis coincides with the Cartesian x 1 -axis. The inset clarifies the cylindrical coordinates used here.
To describe the flow and the free surface shape, we can benefit from two approximations. First, we neglect the inertial effects assuming that the 'inertial pressure' ∼ρv 2 is much lower than the capillary pressure γ/a (here, ρ is density and γ is the surface tension of the liquid): This is a natural condition since a uniform filament is not stable if inertial effects are significant [3]. Second, we assume that polymer stress (σ p ) dominates over the viscous stress of the solvent: σ p η s ∂v ∂z . Again, this condition is obvious well inside the thread (as, otherwise, the polymer effect would be rheologically negligible, and the liquid would be equivalent to the Newtonian solvent there); however, it must be assumed for the transition region. This assumption is justified if the capillary number C n (defined in Equation (71)) is small, as discussed in item 5 of Section 4 (see text around Equation (71)). It is also assumed that the liquid is incompressible: The master dynamical equation is then reduced simply to the force balance: where α, β are Cartesian components, σ αβ = σ αβ (x, t) is the total stress tensor at the point and Y ,β ≡ ∂Y/∂x β is a derivative of a variable Y along the coordinate x β . It involves two contributions due to the pressure field (p = p(x, t)) and the extra polymer stress: The polymer stress tensor is [15]: where F el = F el (R) is the elastic free energy of a polymer chain with end-to-end vector R. Angular brackets here mean the statistical average, c is the concentration of polymer segments (repeat units), and N is their number in one chain (a monodisperse polymer is considered in the present paper). For polymer liquids in marginal or theta-solvent conditions [16,17] (and for concentrated polymer solutions) a Gaussian chain model is applicable during the coil-stretch transition when the chains are far from their full extension, R L (here, L is the chain contour length), so the elastic energy is [15]: where T is temperature in energy units, R coil = N 1/2 b s = √ Ll K is the unperturbed coil size (b s is the polymer statistical segment, l K is its Kuhn segment). Hence, recalling Equation (8), we obtain: where we took into account that all the chains in a fluid element are strongly stretched along the same axis (defined by some unit vector m), so that the end-to-end vector R of most chains is nearly parallel to m, and therefore, R α R β const m α m β =R αRβ , wherē R =Rm, andR is the root-mean-square (rms) of the end-to-end distance. For the sake of simplicity, in what follows, we omit the 'bar' over R and R. (Thus, from now on, R is simply the rms end-to-end distance of polymer chains, which is widely used in polymer physics.) Note that we consider the regime of significant chain extension, R R coil , since, otherwise, σ p αβ would be comparable with the ideal gas pressure of polymer molecules, which is low for long chains. Note also thatR = R(x, t) in Equation (10) is the mean end-to-end vector which generally depends on the chain position x and time t. (More precisely, x is the position of the chain center-of-mass. In all the cases, we assume that the fields such as v(x, t) or R(x, t) do not change much on the polymer length-scale R coil . This condition is satisfied if R coil a 0 .) Obviously, well inside the thread, polymer chains must be stretched along its main axis (z): where R 0 isR (rms of polymer end-to-end distance) in the thread; thus, the only relevant component of the polymer elastic stress tensor is: (other components being subdominant). The model described above is virtually equivalent to the Oldroyd-B dumbbell model, as the polymer conformation is defined solely by the end-to-end vector R, and the polymer stress is quadratic in R. Using this model, it was established that the thinning of a polymer thread in the viscoelastic regime (R coil R L) is quasi-exponential [1,8]: where τ is the stress relaxation time (τ = τ p /2, where τ p is the longest relaxation time of a polymer coil), and t is the time passed since the beginning of coil stretching. Equation (12) (which is applicable in the transient regime during the polymer stretching stage) implies that the thinning rate˙ (cf. Equation (3)) is constant,˙ 2/(3τ). However, strictly speaking, the flow in the thread is not steady, as its radius a 0 = a 0 (t) is decreasing. Nevertheless, the flow in the entrance region (the 'corner flow') can be considered at a different angle: The characteristic thinning time is 1/˙ ∼ τ. Comparing it with the characteristic entering time This means that the entrance flow (at the interface between the thread and a large droplet) is quasi-stationary to a good accuracy, since the ratio a 0 /L f is typically small in experiments, a 0 /L f 0.01 [3,6,7,18]. Let us emphasize again that it is assumed here, following the previous theoretical studies [1,7,8,11,12], that the thread is perfectly cylindrical, the flow velocity in the thread is parallel to the main axis (z), and the chains are stretched along it. The force balance (cf. Equations (6) and (7)) in such a geometry demands that the pressure p is constant (independent of both radial and axial coordinates), p = γ/a 0 , and that the polymeric stress σ p ≡ σ p zz − σ p rr σ p zz is z-independent. In Section 3, we show that well inside the thread, the stress σ p is also constant in the radial direction (cf. text below Equation (56)). Note that, in contrast, σ p does depend on both z and r in the transition regions.
Since the entrance flow (see Figure 2) can be considered as steady, its time dependence can be neglected at the relevant time scale τ ent . The formal boundary conditions on the free liquid surface, defined by r = a(z), are: the surface must tend to a tube of radius a 0 at z/a 0 → −∞ (here and below, z/a 0 → −∞ means that z < 0 and L f |z| a 0 ), and to a sphere of large radius (which can be approximated by a flat plane) at r/a 0 → +∞. (Note that we use cylindrical coordinates z = x 3 , r = x 2 1 + x 2 2 , as shown in Figure 2.) In the limit of no inertia, the equation of motion is reduced to the local force balance, which can be written as (see Equations (6), (7), and (10)). (Most equations below are approximate, just like Equation (10). However, the '=' sign is used there for simplicity, except where it is important to emphasize their approximate nature.) The pressure at the free surface is defined by its curvature C: where a(z) → a 0 at z/a 0 → −∞, a ,zz ≡ ∂ 2 a ∂z 2 → 0 at z/a 0 → +∞ (16) and p is the excess pressure on the top of atmospheric pressure. The flow is therefore defined by an interplay of capillary and viscoelastic effects.
The basic Equations (5) and (14) involve two unknown fields: polymer extension, R = R(x), and velocity v = v(x). However, in the entrance zone, these fields can be easily related to each other. The idea is that the polymer chains must be extended along the stationary streamlines, and their fragments (blobs) must simply follow the stream (i.e., a blob at point x moves with the velocity v(x)). In fact, the relative motions of the blobs due to elastic forces can be neglected as being too slow since the relevant polymer relaxation time (τ p ) is much longer than τ ent (see Equation (13)). As a result, the chain extension vector R stays nearly proportional to the flow velocity v for kinematic reasons: For a steady flow, the above equation follows from the kinematic equations dv/dt ≡ ∂v/∂t + v · ∇v = v · ∇v, dR/dt = R · ∇v (cf. Equation (80) with 1/τ p → 0), which ensure that if R = Bv at some point (which is the case for z < 0, a 0 |z| L f ), the latter equation must stay valid (with B = const ) along the whole streamline in the transition region where the rate 1/τ p can be neglected. Here, d/dt means the full rate of change 'following the fluid' for a material element.
Of note, Equation (17) remains valid in the extended transition region, including the adjacent end part of the thread at z < 0, |z| < ∆, where a 0 ∆ L f . The condition ∆ L f ensures that the transition time (along the whole transition region including the ∆-part) is much shorter than the polymer relaxation time τ p , so that the polymer relaxation process can be neglected in the ∆-part. On the other hand, the condition a 0 ∆ is sufficient to state that at −z ∆ the thread is nearly perfectly cylindrical, so at z = −∆ its radius is a a 0 , the rms chain extension is R R 0 and the flow velocity is v v 0 also being uniform in the cross-section. Then, applying Equation (17) at z = −∆ we obtain B R 0 /v 0 and therefore B is the same for all the streamlines. Equations (10) and (17) obviously agree with the convected Jeffreys (or Oldroyd-B) models with infinite relaxation time [14] (cf. item 11 of the Discussion).
The force balance Equation (14) then becomes: and we took into account the flow incompressibility from Equation (5). Equations (18) and (5), together with the boundary conditions in Equation (15), an obvious condition v n = 0 at the free surface (where v n is the flow velocity component normal to the surface), and v → v 0 at z/a 0 → −∞, p → 0 at z/a 0 → +∞ (20) completely define the flow, the pressure field, and the shape of the free surface in the entrance zone. (Recall that z/a 0 → −∞ physically means that z < 0, a 0 |z| L f ; in a similar fashion, z/a 0 → +∞ means z a 0 . Note also that we consider the limit of a very large droplet, so the excess pressure deeply inside it must tend to 0.) Obviously, Equations (15) and (16) also imply that: There is only one essential non-dimensional parameter in the problem: the ratio X ≡ σ p a 0 /γ of polymer stress σ p to the capillary pressure well inside the filament, p 0 = γ/a 0 . It was set to X = 1 in ref. [8], while X = 2 was derived in ref. [12]. (Note that the quantity X f = 1+X 2 was denoted X in ref. [5]; X f was defined there as the ratio of the net tensile force in the thread to 2πγa 0 . As follows from Table I of ref. [5], the parameter X = 2X f − 1 is as important as X f itself for Newtonian liquids considered there. This parameter (X) is even more important for polymer liquids as it highlights a relation between the main quantities of interest, the polymer stress σ p , and the capillary pressure γ/a.) Below, we show that X can be easily obtained based on the equations considered above. To this end, let us first analyze the pressure p variation along a streamline. Equation (18) gives (with v ≡ |v| and the coordinate u equal to the curvilinear distance along the streamline): (22) leading to ∂p/∂u = Ã /2 ∂v 2 /∂u; hence, This equation is applicable in the entrance zone, where the velocity is changing fast along streamlines. It connects the local pressure and velocity just like the Bernoulli's principle does for an ideal fluid. We shall refer to Equation (23) as the anti-Bernoulli law (due to the 'minus' sign leading to a higher pressure for a faster flow). The const in Equation (23) can be easily found for the flow of Figure 2: Recall that each streamline is coming from the uniform thread, where both p and v are constant over a cross-section. Therefore, in the whole entrance zone (|z| L f ). Note that Equations (24) and (23) are valid if the Newtonian solvent friction is negligible as assumed below Equation (4) and as discussed in item 5 of the Discussion section (cf. Equation (71)). The rhs in the above equation equals 0 since both p and v must vanish deeply in the droplet region at z a 0 (cf. Equation (20)). Thus, Equations (20), (21), and (24) imply that and we obtain σ p in the thread using Equations (11) and (19): Hence, the parameter X introduced below Equation (21) is: The total thread tension force T (defining the total momentum flux through its cross-section) is therefore: The two terms in the middle of the above equation represent the surface and the internal axial stress contributions to the total tension T . Equations (27) and (28) agree with the results of ref. [12]. It is worthwhile to note that Equations (17)-(24) are valid for a virtually arbitrary flow geometry (involving a long channel as a part of it). Using Equations (18) and (24), we obtain: Let us now take into account that, typically, the experimental setup is such that both the free surface and the flow are axisymmetric. This symmetry is already implied in Equation (15). For the flow velocity, it means that its axial and radial components (v z and v r ) do not depend on the polar angle θ, while the polar velocity component v θ must vanish (here, we use cylindrical coordinates (z, r, θ), cf. Figure 2): The factor in brackets in Equation (29) is related to the vorticity where αβγ is the Levi-Civita symbol. Therefore, Equation (29) is equivalent to: For an axisymmetric flow, ω is always parallel to the azimuthal direction (ω r = ω z = 0); hence, ω and v are orthogonal. Equation (32) then gives ωv ≡ 0, so: everywhere except at the stagnation points (where v = 0), which are not expected in the entrance zone. Therefore, the flow is irrotational: where ϕ = ϕ(r, z) is a potential field which must satisfy the Laplace equation (in view of Equation (5)): Furthermore, on using Equations (24) and (25), the boundary condition, Equation (15), transforms to: while another condition at the free surface (v n = 0) simply says that the boundary curve r = a(z) must be a streamline with r → a 0 at z/a 0 → −∞ (obviously, the axisymmetric problem is essentially two-dimensional, so we can set θ = 0 and r = x 1 without loss of generality). Further statements, which can facilitate numerical solution of the above equations for the velocity field v(x), can be proven in an elementary way. At large z a 0 , the flow is asymptotically similar to the electrostatic field of an electric charge: v Qx where x is the position vector, r = |x| is the distance from the origin, and Q is obtained based on the incoming volume flow rate, which is equal to πa 2 0 v 0 . Equation (37) implies that at r a 0 , the pressure decreases as: The pressure therefore becomes negligible in the droplet far from the entrance point, so the droplet shape there must be close to the static shape defined solely by the Laplace pressure. It leads to the condition of constant mean curvature, C → 0, in the case of large droplet, so [4]: where a ∼ a 0 . The above equation is consistent with Equation (3).18 of ref. [12].
Next, let us note that the structures of Equations (35) and (36) are such that a simple rescaling r → r/a 0 , z → z/a 0 , and v → v/v 0 makes the problem free of any parameters (in other words, the solution is self-similar). It is therefore enough to solve the above equations for a 0 = 1, v 0 = 1 and Q = 1/2. The unique rescaled flow field and the free surface curve a(z) can be found numerically in two steps: (1) by setting a trial a(z) and then solving the Laplace Equation (35) in the region z > −z max , r < r max , and r < a(z) (see Figure 2) with Neumann boundary conditions: ∇ϕ · n = 0 at the free surface (r = a(z)), ∇ϕ · n = 1 in the filament cross-section at z = −z max , and ∇ϕ · n = const at the hemisphere r = r max inside the droplet (here n is unit vector normal to the surface); (2) by adjusting a(z) iteratively so as to satisfy the anti-Bernoulli equation . The resultant thread shape and streamlines in the entrance zone (obtained for z max = 6 and r max = 30 taken to obtain a well-converged solution to an accuracy of ∼0.5%) are shown in Figure 3. In the thread region, z → −∞, the free surface tends to a cylindrical shape in an exponential fashion: where k ≈ 1.393. (The constant k > 0 is the lowest root of the characteristic equation 2kJ 0 (k) = (1 + k 2 )J 1 (k), which comes from the main correction to the potential field ϕ at −z 1: ϕ − ϕ 0 const J 0 (kr)e −kz , where ϕ 0 = z corresponds to a perfectly uniform flow. Here J 0 , J 1 are Bessel functions.) The asymptotic behavior in the opposite regime, z 1, can be deduced from the force balance, which says that the total momentum flux in the cylindrical part (the thread tension T = 3πγ, cf. Equation (28)) must be equal to the momentum flux, T + , across a hemisphere r = const 1, z > 0. The contribution of both polymer stress and pressure to T + can be neglected since they rapidly tend to zero at large r: σ p ∼ p ∝ 1/r 4 (cf. Equation (38)); hence, pr 2 → 0 at r → ∞. Therefore, T + at r 1 is defined solely by the surface tension γ: The equation T + = 3πγ leads to This asymptotic behavior was predicted [12] using a similar argument. However, the point that the polymer stress can be neglected in the droplet region was not proven there. Of note, the numerical solution for a(z) shown in Figure 3 is in harmony with the asymptotic laws, Equations (40) and (42).

FENE-P Model
The theory described above was developed for the Oldroyd-B model. Below, it is generalized to account for nonlinear finite extensibility effects. In the general case, the polymer chain tension ∂F el ∂R depends on the end-to-end distance R in a nonlinear fashion (cf. Equation (9)): where s = R/L is the stretching degree (whose maximum value is 1), and R 2 The factor κ FE (s) must tend to 1 at low s, when the Gaussian model is valid (recall that we consider polymers in marginal or theta solvents). By contrast, κ FE must generally strongly increase for 1 − s 1 to provide a finite extensibility (FE): κ FE → ∞ at s → 1. The whole function κ FE (s) is, however, not universal: it depends on the polymer flexibility mechanism. One of the most popular FE models was introduced by Warner [19]. It serves as a part of the FENE-P dumbbell model [14,20,21] with Using Equations (43) and (44), we obtain: where N K = L/l K = L 2 / Nb 2 s is the number of Kuhn segments in a polymer chain. The elastic force, Equation (43), defines the polymer stress tensor, cf. Equation (8): where G = cT/N is the shear elastic modulus defining the polymer elastic response at short times (t τ). (Note that internal relaxation modes are ignored in the FENE models.)

Equations of Motion in the Entrance Zone
We now turn to dynamical equations applicable in the transition zone where the flow is non-uniform but is virtually steady (cf. Equation (13) and the text below it). In this regime, as argued above Equation (17), the polymer extension (the end-to-end vector) R is proportional to the flow velocity v, cf. Equation (17). Note that polymer contraction due to elastic force can be neglected in the entrance region even if the chains are stretched close to full extension (s close to 1) because s rapidly decreases hydrodynamically (affinely with the fluid element) down to s ∼ 0.5 upon entering the transition region due to a fast deceleration of the flow there. Indeed, the rate of affine contraction of chains in the transition region is defined by the flow velocity gradient, |∂v z /∂z| ∼ v 0 /a 0 ∼ (L f /a 0 )τ −1 ; it is much faster than the rate (∼1/τ) for the non-Gaussian elastic relaxation, since the factor L f /a 0 (the ratio of the thread length to its radius) is very large, as stated in Equation (2) (and as observed experimentally). Therefore, considering the entrance flow for the FENE-P model, we can treat the flow as quasi-stationary setting τ → ∞ and using Equation (17) just like for the Oldroyd-B model.
The factor B in Equation (17) is constant along a streamline, but strictly speaking, it may vary among the streamlines. In Section 2, we provided a plausible argument showing that B is constant in the whole entrance zone. This property is proved below on more general grounds. Remarkably, even with a variable B equation, R = Bv still ensures that the field R(x) is solenoidal: since v and ∇B are necessarily orthogonal. Using Equations (6), (7), (46), and (47) we obtain: Recall that we neglect inertial effects and the stress contribution due to the solvent viscosity for the reasons considered in Section 2. Next, taking into account that F el can be considered as a function of R 2 , one obtains: where The above equations lead to: Considering the pressure gradient along a streamline, using Equation (51), we obtain: where the function F = F (R) must satisfy the equation: and u is a curvilinear coordinate (the contour length) along a streamline. From Equation (53), we obtain: so that Equation (52) leads to: along each streamline. Taking into account that F (0) = 0 and that both p and the chain extension R must vanish in the droplet far away from the thread, we find that const = 0 in Equation (55), so: Equation (56) is valid in the whole entrance region including the thread-end region (z < 0, a |z| L f ), where p p 0 and v z v 0 are nearly uniform both axially and radially (in a cross-section). Using Equation (56), we find that the same must be true for the chain extension R and the polymer stress σ p related to it. Recalling that both σ p and R = R 0 are z-independent within the thread (see text below Equation (13)), we conclude that these quantities do not change in the radial direction in the whole thread. This means, in particular, that the factor B = R 0 /v 0 in equation R = Bv (cf. Equation (17)) is the same for all streamlines. Equation (56) can be therefore written in the form: which generalizes the anti-Bernoulli law, Equation (24). For the FENE model, Equations (45) and (54) give: where R = sL andF In the regime s 1, both F (R) and the elastic energy density F el (R) = c N F el (R) are nearly quadratic in R and are equal: Therefore, in this case, the anti-Bernoulli law also says that the local pressure is equal to the elastic energy density. However, in the opposite regime of nearly fully stretched chains, 1 − s 1, F (R) becomes much larger than F el (R): Based on the anti-Bernoulli law, we have shown in Section 2 that the flow must be irrotational for the Oldroyd-B model: ω ≡ ∇ × v = 0 (cf. Equation (33)). This property is not preserved in the more general case considered in this section: it is replaced by a more general relation derived below (cf. Equation (63)).
For symmetry reasons, ω has only one nonzero component (= ω) in the azimuthal direction perpendicular to the rz-plane. Using local coordinates with axis u tangential to a streamline (at a given point O) and axis u ⊥ perpendicular to it (cf. Figure 4), we obtain: where v ⊥ is the velocity component along axis u ⊥ . Noting that ∂v ⊥ /∂u is proportional to the local curvature, C, of the streamline: and using Equations (54) and (56), we find: Comparing the above equation with Equation (51) and recalling Equation (17), we finally obtain: Hence, in the general case, ∂v/∂u ⊥ = ∂v ⊥ /∂u and ω = 0. Note, however, that Equations (62) and (63) still ensure that ω = 0 for a straight streamline, C = 0. This is true, in particular, both well inside the droplet and in the thread. Therefore, Equations (37) and (38) remain valid at r a 0 .

Flow in the Thread and the Thinning Law
Turning to the flow in the cylindrical part of the thread, we first recall that it is uniform in a cross-section, so both the pressure p and the polymer stress σ p αβ do not depend either on the axial coordinate z or r in the thread. The pressure p = p 0 in the thread is defined by Equation (56), where R = R 0 is the chain extension: On the other hand, the polymer stress σ p = σ p zz − σ p rr is (cf. Equation (46)): The s-dependence of the ratio X = σ p /p 0 is shown in Figure 5. It is clear that σ p /p 0 is not a constant: for s 1, the above equations lead to σ p /p 0 = 2, while for 1 − s 1, the result is σ p /p 0 = 1. The kinetics of the chain extension for the FENE-P model can be described with evolution equation: Here,˙ is the extension rate in the axial direction, which is related to the thread thinning rate, see Equation (3), and R = R 0 . The second term in the rhs of Equation (66) is proportional to the elastic force [20]. (Note that a thermal noise term is not present in Equation (66) It leads (using Equation (44)) to the following ODE for s = s(t): Equation (68) is valid for the FENE-P model. It can be integrated to give (with substitution y ≡ s 2 ): − ln y + y + 4 1 + y 2y + (1 − y) ln(1 − y) The above equation defines the time dependence of the stretching degree, s(t), which is shown in Figure 6. The time dependence of the ratio σ p /p 0 is indicated there as well.
Obviously, this ratio strongly decreases near the breakup point (the const in Equation (69) was chosen to obtain the breakup at t = 0). The resultant thinning law, a 0 (t) = γ/F (R 0 ) = 2γ 3GN K 1 F (s) , is drawn in Figure 7. The straight dashed line highlights the fact that a 0 (t) follows the classical exponential law, a 0 (t) ∝ exp − t 3τ , [1], as long as s = R 0 /L 0.5.

Discussion
1. In the present paper, we investigated the dynamics of a viscoelastic liquid bridge (a cylindrical ligament) connecting two large droplets (see Figure 1). Its thinning is governed by capillary and viscoelastic forces and is coupled to the flow both in the ligament and in the transition zone between the thread and a droplet (see Figure 2). As a fluid, we consider an unentangled polymer solution whose rheology can be described by either the Oldroyd-B or the FENE-P dumbbell models. The developed theory is focused on the regime of strongly stretched polymers whose conformational tensor R α R β is dominated by a single eigenvalue (here, R is the end-to-end vector of a polymer chain). The models we used assume that the stress tensor contribution of a polymer chain is defined by the vector R and therefore is proportional to R α R β . The prefactor in such a relation is either constant (cf. Equation (10) for Oldroyd-B model) or depends on the degree of chain stretching s = R/L (for the FENE-P model). Both models ignore hydrodynamic interactions (HDI) in the system (see point 9 below).
As one of the main results, we have shown that for the Oldroyd-B model, the capillarydriven flow is nearly irrotational and satisfies the Laplace equation. In addition, we have established a general relation between the pressure and velocity fields (valid for a broad class of polymer flows) which is referred to as the anti-Bernoulli law (cf. Equation (23); an equivalent static conservation law was derived for a neo-Hookean elastic solid [12]). Using these results, we obtained a self-similar solution (involving a single length-scale, the radius a 0 of the thread) for the shape of the free surface and for the pressure and velocity fields. In cylindrical coordinates (r, z), the free surface is given by r = a(z) = a 0 f (z/a 0 ). The numerically obtained free surface profile a(z) for a 0 = 1 and the corresponding flow structure are shown in Figure 3. The calculated a(z) is compared in Figure 8 with the similarity solution for the bead-string structure of a neo-Hookean elastic body obtained in ref. [12]. (Note that l e = 2 1/3 a 0 was taken as a unit length in ref. [12], so the data have been rescaled accordingly.) The obvious good agreement supports the conclusion of refs. [12,13] that (as conjectured earlier in refs. [22,23]) the neo-Hookean model asymptotically reproduces the geometry of a polymer liquid bridge formed in the course of its capillary thinning (according to the Oldroyd-B constitutive equation) in the regime of sufficiently long polymer relaxation time τ. Of note, it was recently shown that the universal interfacial shapes obtained in ref. [12] (and, therefore, the theoretical profile of Figures 3 and 8) are in agreement with high-resolution experimental data on aqueous solutions of high-molecular-weight polymers (PEO) and biopolymers (hyaluronic acid) [24,25]. 2. The conservation (anti-Bernoulli) law defining pressure in terms of the flow field (cf. Equation (24)) was obtained here based on the steady-flow dynamical equations (Equations (6) and (17)). Of note, in the form of Equation (56), p = F (R), this law is akin to the conservation law established for neo-Hookean solids [12]. The main difference is that Equation (56) is more general: it also accounts for non-linear elastic effects related to the finite extensibility of polymer chains.
3. To quantitatively establish the thinning kinetics, a 0 = a 0 (t), one needs a 'closure' relation between the pressure p 0 and the axial polymer stress σ p in the thread. This relation cannot be obtained by considering a cylindrical thread as such. Long ago, Entov and Hinch [8] assumed that the total axial stress in the thread (σ p − p 0 ) must be equal to that in the droplet, i.e., to zero (if atmospheric pressure is subtracted from p 0 ): p 0 = σ p . It is found here that the rigorous anti-Bernoulli law, Equation (24), demands a different relation, p 0 = σ p /2, in agreement with more recent predictions for the Oldroyd-B model [7,11,12]. It is remarkable that a generalization of the anti-Bernoulli law for polymer chains with finite extensibility, Equation (57), leads to a variable ratio X = σ p /p 0 which depends on the degree of polymer extension s = R/L. While for s 1 we obtain σ p /p 0 = 2, this ratio decreases with s (see Figure 5), reaching its minimum, σ p /p 0 = 1, in the regime of completely stretched chains, s → 1. Incidentally, it is the latter ratio that was assumed in ref. [8].
It may seem that the model of a nearly uniform thread cannot be used for s close to 1 since in this regime of highly extended chains the liquid behaves nearly as a Newtonian fluid with renormalized viscosity [4,[26][27][28]] η * η s 1 + π 18k H c N L 3 , where k H 0.5 ln(1/φ) is a hydrodynamic factor [29] and φ is volume fraction of polymer. Thus, the classical Plateau-Rayleigh instability (undulations of the thread shape) must emerge in this regime (at s close to 1, s > s c ∼ 0.5), so that the uniform thread approximation seems to fail. The point, however, is that: (i) Before the instability (at s < s c ), the undulations are weak, and their amplitude (in radial direction) due to thermal fluctuations isã 0 ∼ T/γ. (Here we make a plausible assumption that the amplitude of thermal fluctuations of the radius in the dynamically stable thread are comparable to that for a statically stable liquid cylinder of the same diameter.) Therefore,ã 0 a 0 , since, typically, T/γ 1 nm, while a 0 1 nm. (ii) The growth rate of their amplitude (coming from the Rayleigh theory [30]) is |d(lnã)/dt| ∼ 1 6 γ aη * , so it is comparable with the thread thinning rate, [4,5]˙ /2 ∼ 1 6 γ aη * . Therefore, as the thread thins by a factor of k, the undulations grow by the same factor, so the thread radius a becomes comparable with the undulation amplitudeã for k = k * ∼ √ a 0 /ã 0 1. Considering a = a(s) as a function of s, we find, according to the theory of Section 3.3 that a(s)/a(0.5) ∼ 1 − s for s 0.5, that is, 1/k ∼ 1 − s. The above argument therefore shows that the thread undulations remain small (ã a) if 1 − s 1/k * . Hence, the theory of Section 3.3 is actually applicable in the regime where the chain extension degree s is close to the maximum (=1): this is true in the region 1 1 − s 1/k * since k * 1. In this regime, the ratio X indeed becomes very close to 1, as stated in the previous paragraph.
4. The developed theory was applied to obtain the thread thinning dynamics for the FENE-P model. While this problem was treated already in ref. [8], it is treated here for the first time in a rigorous way. The numerically obtained thinning law, a 0 (t) (see Figure 7) includes both the elasto-capillary (exponential thinning) and the terminal stage (quasi-linear thinning). It is based on the asymptotically exact dependence of the ratio X = σ p /p 0 on the degree of chain stretching established in Section 3 (see also point 3 above).
The predicted thinning law is universal (up to arbitrary horizontal and vertical shifts) when plotted in semilog scale vs. t/τ. The theory is compared with experimental data on semidilute, unentangled polymer solutions (cf. Figure 5 of ref. [18]) in Figure 9, where t/τ = 0 for the theoretical curves is set to be the putative breakup point (a 0 = 0). One can observe a good agreement down to the crossover regime; however, at later times (on the right to a red circle), the experimental curve deviates from the prediction. We attribute this discrepancy to a new regime most probably associated with the blistering (pearling) instability [2,31], which is normally accompanied by the polymer/solvent phase separation and formation of secondary solvent droplets (which can be both flow-and capillary-induced [28,29,[32][33][34]). Interestingly, according to Figure 9, this transition occurs where the chains are not yet fully stretched, at s ≈ 0.5, corresponding to X = σ p /p 0 ≈ 1.75. 5. Considering the thinning of a cylindrical polymer thread, we assumed that the pressure p there is constant in the radial direction (i.e., it is homogeneous in a cross-section). This assumption was adopted in the previous theoretical studies on capillary thinning [1,7,8,11,12]. In the case of negligible inertia, it simply follows from the force balance in the cylindrical geometry since the flow velocity v is parallel to the main axis, and so is the elastic force due to stretched polymer chains. In addition, however, the theory developed in Sections 2 and 3 implies that the axial polymer stress σ p is also uniform (independent of r) in the thread. For the Oldroyd-B model, this condition is justified in the limit of vanishing solvent viscosity as follows: First, it is straightforward to show that the argument (cf. Section 2) leading to the anti-Bernoulli law, Equation (24), is generally valid in the form: (where σ p = σ p uu is the polymer stress component along the streamline). Obviously, in the thread, σ p σ p and (as argued above) p = const . Therefore, Equation (70) implies that σ p is also constant (r, z-independent) there. The conjecture that for negligible solvent viscosity the stress profile in the thread is uniform [12] is thus rigorously proved using a dynamical approach. (Of note, it was also established [12] that the radial distribution of the axial stress is uniform for a soft neo-Hookean elastic solid. Hence, here, again, we confirm the correspondence between viscoelastic and elastic models [13].) In Section 3.2, we show that this statement is also valid for the models with finite chain extensibility.  Figure 5 of ref. [18]). Red circles indicate the onset of significant deviations between theoretical and experimental curves. The theoretical breakup points are indicated with red arrows for each red curve. For the two highest concentrations, we also indicate the reduced time t/τ (with t = 0 corresponding to the theoretical breakup). Thin red dashed lines indicate the exponential asymptotic behavior at short times (according to Equation (12)).
Numerical studies of the Oldroyd-B model with a finite Newtonian solvent viscosity η s [12] reveal a weak radial dependence of σ p , which was attributed to an effect of η s . We tend to agree with this interpretation. The effect of η s (which was neglected in the present study) is unimportant as long as η s ∂v ∂z σ p , which is true everywhere (including the entrance region) if a capillary number is small: The above condition was adopted in the previous sections. Note that for any η s , this number (C n ) can be made however small, provided that the polymer time τ is sufficiently long. For a finite η s , the polymer stress σ p in the thread may vary along the radius. We expect that (at least for a small C n ) a maximal relative variation of σ p in the radial direction should be defined by C n . It may also be worth mentioning that an increase in the solvent's viscosity does not necessarily lead to a higher C n since the ratio η s /τ depends mainly on the polymer molecular weight rather than on η s (note that τ is nearly proportional to η s if polymer volume fraction is low [15,35]). 6. It is well-known [7,11,12,36] that the elasto-capillary thread thinning (according to an exponential law, Equation (12)) is preceded by a fast initial process of the thread formation, which is primarily governed by capillary and inertial forces. In a typical scenario, a cylinder of initial radius R 0 and length > 2πR 0 undergoes the Plateau instability. A thread of radius a 0i R 0 is formed as a result of this process [7,11,12,36]. The characteristic time of this initial stage is inertial in nature, t i ∼ ρR 3 0 /γ 1/2 . In the inertial regime, t i also defines the period of oscillations of the emerging semi-droplets (connected by the thread). Their damping rate is proportional to the solvent viscosity η s defining the dissipation rate, so the damping time is t damp ∼ t 2 i γ/(η s R 0 ) ∼ ρR 2 0 /η s . (Note that t damp t i as long as R 2 0 η 2 s /(ργ), which literally means that initially the system falls in the inertial regime.) This time must not exceed τ (since, otherwise, the thread thinning would be strongly perturbed by droplet oscillations). Hence, we have to demand that This condition is compatible with Equation (71) if which simply means that t i τ. Note that the condition (73) also ensures that inertial effects are negligible in the thread: ρv 2 0 γ/a 0 , as was assumed in Section 2 (cf. Equation (4)). Thus, Equations (71) and (72) define the range of allowed η s : (Here, we assumed that the filament length is comparable with R 0 , L f ∼ 3R 0 , which was the case in the previous numerical studies [7,12,23].) Obviously the conditions of Equation (74) are always satisfied for a sufficiently long polymer relaxation time τ. 7. The fluid dynamics considered here are dominated by the longitudinal polymer stress σ p (along a streamline). This property is quite natural in the regime of strongly stretched polymer chains, when the polymer chain end-to-end distance R is much larger than the equilibrium coil size R coil , R R coil . In fact, it was implicitly assumed that polymer chains are strongly stretched during an initial fast inertio-capillary stage of polymer thread development [36], i.e., before the elasto-capillary regime of the thread thinning considered in this paper. Such initial process leads to a uniaxial polymer stress field, cf. Equation (10), both in the thread (of radius a 0 ) and in the droplet/thread transition region of size ∼ a 0 . However, the applicability of this approximation in the droplet at large distances from the transition region, r a 0 , may be questioned. Indeed, R ∝ v ∝ 1/r 2 (cf. Equations (17) and (37)) strongly decreases there with the distance r, so the chains become less stretched in the radial direction (along the streamline), and, more importantly, they become expanded in the lateral directions: the lateral size R ⊥ grows as R ⊥ ∼ R coil r/a 0 . It leads to an emergence of the lateral stress σ ⊥ ∼ G(r/a 0 ) 2 . The stress growth stops as soon as the travel time t v (r) for a fluid element to reach the distance r from the origin becomes comparable with the stress relaxation time τ. Using Equation (37) with v 0 =˙ L f /2 = L f /(3τ), we obtain t v (r) ∼ r 3 /(v 0 a 2 0 ), so that t v ∼ τ is reached at r = r * ∼ L 1/3 f a 2/3 0 .
The maximum lateral stress is, therefore, σ ⊥ ∼ G L f /a 0 2/3 . This stress should lead to an additional curvature ∆C ∼ σ ⊥ /γ of the free surface, which can be neglected if ∆Cr 1, leading to the condition: Here, E c is nearly equivalent to the elasto-capillary number introduced in refs. [7,12] (with L f ∼ 3R 0 ). It is noteworthy that the same condition, E c 1, also ensures that the chains are strongly stretched during the initial inertio-capillary process: R R coil , a 0 L f (cf. refs. [7,12]).
8. To sum up the previous two points, there are three main conditions of validity of the theory: Equations (71), (74), and (75). It is noteworthy that Equations (71) and (75) are similar and are actually equivalent in the dilute/semidilute transition regime since the polymer contribution η p to the total viscosity is η p = Gτ (for the Oldroyd-B model) and η p ∼ η s at c ∼ c * , where c * is the coil overlap concentration. Furthermore, in the semidilute regime (c > c * ), Equation (71) follows from Equation (75).
9. In the present study, we totally neglected the hydrodynamic interactions (HDI), which are not incorporated either in the Oldroyd-B or FENE-P models. Such a simplification, which was also adopted in the previous theoretical works on the subject [1,7,8,11,12,37], may be appropriate in the semidilute or concentrated solution regime (c > c * ). The effects related to topological constraints (entanglements between polymer chains) [38] are disregarded here as well, so the region of applicability should be formally restricted to the unentangled polymer regime. It is worth noting, however, that both models employed here seem to also work well for dilute polymer solutions, at least as far as the capillary thinning kinetics are concerned [8,28,32,39]. Moreover, the main effect of the HDI is to modify the polymer relaxation time, which is completely irrelevant in the entrance zone with high deformation rate. Therefore, the approach to describe the flow in such zones (cf. Sections 2 and 3.2) also remains valid for dilute polymer solutions. 10. In this paper, we ignored any effects of polymer/solvent separation assuming that polymer concentration c = const . It is well-known that an inhomogeneous polymer distribution in a solution may be generated by a flow, for example, due to the stressconcentration coupling (SCC) effect [32,[40][41][42]. This effect is not relevant for the cylindrical thread where the stress is uniform, so the concentration is uniform as well. Moreover, while in principle, the SCC effect is present in the entrance zone, it is very weak there: in fact, we argued (cf. Section 2) that polymer chains are extended along the streamlines which are curved in the entrance zone. This curvature generally leads to a net elastic force perpendicular to the streamline, f ⊥ , acting on an internal part of the chain, f ⊥ ∼ f el R/a, where R is the chain length, 1/a is the typical streamline curvature, and f el = ∂F el /∂R is the elastic tension of the chain (cf. Equation (9)). The force, f ⊥ , leads to the perpendicular velocity, v ⊥ , of the chain relative to the solvent, v ⊥ ∼ (R/τ)(R/a). The typical time the chain spends in the entrance zone is t ∼ a/v, so its typical lateral displacement is Taking also into account that R a L f , we find that u ⊥ /a 1, so the SCC effect is indeed negligible and c const in the system we considered.
11. The polymer stress tensor defined in Equation (8), with the conditions specified before Equation (10), can be written as: Its rate-of-change prescribed by the Oldroyd-B model can be expressed in terms of the upper convective derivative: as (here, d/dt ≡ ∂/∂t + v · ∇ is the full rate-of-change referred to a moving physical element of the fluid). The above equation is equivalent to the standard Oldroyd-B equation (see, e.g., Equation (2.2) in ref. [12]). (Note that σ p αβ here is the full polymer stress, while the polymer stress σ (p) in the standard Oldroyd-B model does not include the isotropic equilibrium polymer stress contribution δ αβ cT/N, so σ p αβ = σ (p)αβ + δ αβ cT/N.) To further simplify the equations, we took into account that, in the thread (and, in particular, near its mid-point), all the chains are strongly stretched in a similar way along the filament axis (z-axis). Therefore, the polymer stress tensor there must be dominated by a single eigenvalue corresponding to an eigenvector R parallel to this axis, as argued around Equation (10): The prefactor in the above equation is such that it makes sure that |R| is equal to the rootmean-square (rms) end-to-end distance of polymer chains in a fluid element. Equations (77) and (78) (76) and (79) together are equivalent to Equation (10)), it will be also valid at all the points downstream (as long as the polymer chains remain strongly stretched in a single direction), with vector R changing according to where τ p = 2τ. To derive Equation (80), we neglected the second term in the rhs of Equation (78): this term is small because |σ p | cT/N since the polymer (the dumbbell) is strongly stretched.
In the case of a steady flow, it is enough to know that Equation (10) is applicable in a cross-section of the thread: then, it must also be valid in the whole volume downstream (as long as the eigenvalue of the tensor σ p associated with the eigenvector R remains strongly dominating). Note also that (again for a steady flow) the full rate of change of velocity is: Interestingly, Equations (80) and (81) show that if R = Bv at some point at t = 0, the two vectors remain parallel at all points downstream (along the streamline which includes the initial point): R = B(t)v with B(t) = B(0) exp(−t/τ p ). The validity of this statement can be demonstrated simply by substitution of the solution, R = B(t)v, in Equation (80) to make sure that this equation is satisfied using Equation (81).

Summary
We studied the thinning dynamics of a liquid bridge containing long flexible polymer chains in the viscoelastic (elasto-capillary) regime where the thread diameter decreases according to the classical exponential law and the chains are highly stretched with respect to their equilibrium coil size, so that the polymer stress tensor can be approximated by a vector dyad, σ p αβ ∝ R α R β . It is shown that the liquid flow in the transition zone between the thread and a large end-droplet can be considered as quasi-stationary and that the time spent by a polymer chain in this zone is much shorter than the polymer relaxation time τ. Moreover, it is rigorously demonstrated based on the full Oldroyd-B equations (cf. item 11 of the Discussion) that if the polymer stress tensor in a material element is a vector dyad initially, it will remain a dyad provided that the polymer relaxation time is very long, τ → ∞. This statement is valid for any flow, whether it is irrotational or not, steady or not, and for any initial orientation of R. If, in addition, the flow is steady and initially the polymer end-to-end vector R is parallel to the velocity v, then R will stay parallel to v and proportional to it at all points down the streamline. The validity of Equation (17) in the transition zone is justified this way.
Using the Oldroyd-B model in the transition zone, we established a general relation between the pressure and the flow velocity for the case of negligible solvent viscosity η s . The relation says that the excess pressure is proportional to the square of velocity. Using this relation (termed the anti-Bernoulli law) and the axial symmetry of the flow, we also found that the flow must be irrotational. These results allow us to obtain the flow field and the free surface shape in the transition zone using an obvious electrostatic analogy (cf. Figure 3). The obtained surface profile is in good agreement with asymptotically exact Equations (42) and (40) and with recent numerical results for a similar model [12] obtained with a different theoretical approach (cf. Figure 8).
Using the developed theory, we show that the ratio of the polymer normal stress difference σ p to the capillary pressure p 0 in the thread is X = σ p /p 0 = 2 in the elastocapillary regime if the polymer relaxation time τ is sufficiently long (more precisely, when both capillary numbers, C n and E c , are small, cf. Equations (71) and (75)), which is in agreement with recent theoretical studies [11,12]. Furthermore, it is also proven that σ p is uniform in a cross-section of the thread, σ p = const , in the limit η s → 0 or for a long polymer relaxation time, τ → ∞.
The proposed theory (including the anti-Bernoulli law, cf. Equation (57)) is also generalized to account for the case of polymer chains with finite extensibility (the FENE-P dumbbell model). It shows (cf. Section 3) that the thread thinning turns significantly faster than the classical exponential law, i.e., Equation (12), if the degree of chain stretching s = R/L exceeds s ∼ 0.4. In this regime, the maximum flow velocity v 0 considerably increases with respect to the classical prediction v 0 = L f /(3τ), while the ratio X = σ p /p 0 decreases down to X 1.
The developed theory paves the way to consider other flow regimes (such as a flow out of a droplet into a filament) and other rheological effects, for example, the effects of solvent viscosity η s or of an increasing filament length L f for the capillary thinning dynamics. We expect that the exponential thinning law, Equation (12), should be significantly modified if η s /τ γ/L f (i.e., C n 1). In this regime, the solvent viscosity η s must be very important for the transition regions near the thread ends. Such effect was considered to some extent based on the Onsager variational principle [11,37]. However, the effect of η s on the flow field in the transition regions was not analyzed in these studies. This effect could be a subject of a separate publication.