1 Introduction
An illustration of the problem studied in this paper. In (a), a sequence of experimental images showing the pinch-off of a drop of polymer solution (polyacrylamide) between two solid plates, forming a uniform thread of radius 
$h_{thr}(t)$ (Deblais, Velikov & Bonn Reference Deblais, Velikov and Bonn2018b). In (b,c), our numerical simulation of the Oldroyd-B model for flexible polymer solutions, with dimensionless parameters (defined in § 2.3), chosen as in Turkoz et al. (Reference Turkoz, Lopez-Herrera, Eggers, Arnold and Deike2018): 
$De=60$, 
$Oh=3.16$, 
$S=\unicode[STIX]{x1D702}_{s}/\unicode[STIX]{x1D702}_{0}=0.25$. The elasto-capillary number is 
$E_{c}=\unicode[STIX]{x1D702}_{p}R_{0}/(\unicode[STIX]{x1D706}\unicode[STIX]{x1D6FE})=0.0395$. Times are made dimensionless with the capillary time scale 
$\unicode[STIX]{x1D70F}$, and lengths are given in units of the initial bridge radius 
$R_{0}$. In (b), 
$h_{thr}(t)$ is shown to decrease exponentially in time, with decay exponent 
$1/3\unicode[STIX]{x1D706}$, where 
$\unicode[STIX]{x1D706}$ is the polymer time scale (blue line). For comparison, we also include a simulation with infinite relaxation time, keeping 
$E_{c}$ constant. In (c), a close up of the corner region between the thread and the drop, showing the axial stress 
$\unicode[STIX]{x1D70E}_{zz}$. The stress is very high inside the thread, and decays very rapidly in the transition region.
When a drop of water falls from a tap, it will break into two or more pieces under the action of surface tension (Basaran Reference Basaran2002; McKinley Reference McKinley2005; Eggers & Villermaux Reference Eggers and Villermaux2008). Similarly, holding a drop between two solid plates which are then separated to make this liquid bridge unstable, breakup is observed at some finite time 
$t_{0}$. Near breakup time, drop profiles and the velocity field are characterised by similarity solutions, which describe the evolution toward smaller and smaller scales.
If even a small amount of long, flexible polymer is added, the singularity is inhibited (Bazilevskii et al. Reference Bazilevskii, Voronkov, Entov and Rozhkov1981; Chang, Demekhin & Kalaidin Reference Chang, Demekhin and Kalaidin1999; Amarouchene et al. Reference Amarouchene, Bonn, Meunier and Kellay2001; Clasen et al. Reference Clasen, Eggers, Fontelos, Li and McKinley2006; Morrison & Harlen Reference Morrison and Harlen2010; Eggers & Fontelos Reference Eggers and Fontelos2015), and a thread of almost uniform thickness is formed, as shown in figure 1. This is because polymers are stretched in the extensional flow near breakup, but resist the stretching, producing an axial stress, sometimes quantified as an extensional viscosity. This slows down the pinching, and results in a uniform thread thickness. In the framework of the Oldroyd-B model for polymer solutions, the characteristic relaxation time 
$\unicode[STIX]{x1D706}$ of the polymer selects a constant extension rate 
$\dot{\unicode[STIX]{x1D716}}$ inside the thread, which leads to an exponential thinning
$$\begin{eqnarray}h_{thr}=h_{0}\exp \left(-\frac{t}{3\unicode[STIX]{x1D706}}\right),\end{eqnarray}$$ of the thread radius (Bazilevskii, Entov & Rozhkov Reference Bazilevskii, Entov, Rozhkov and Oliver1990; Clasen et al. Reference Clasen, Eggers, Fontelos, Li and McKinley2006; Eggers & Villermaux Reference Eggers and Villermaux2008). This is well confirmed by experiments (Bazilevskii et al. Reference Bazilevskii, Entov, Rozhkov and Oliver1990; Anna & McKinley Reference Anna and McKinley2001; Amarouchene et al. Reference Amarouchene, Bonn, Meunier and Kellay2001; Clasen et al. Reference Clasen, Eggers, Fontelos, Li and McKinley2006) in both high and low viscosity solvents, as well as in full numerical simulations of the Oldroyd-B equations (Etienne, Hinch & Li Reference Etienne, Hinch and Li2006; Bhat et al. Reference Bhat, Appathurai, Harris, Pasquali, McKinley and Basaran2010; Morrison & Harlen Reference Morrison and Harlen2010; Turkoz et al. Reference Turkoz, Lopez-Herrera, Eggers, Arnold and Deike2018), which have become a standard model for the description of solutions of long, flexible polymers (Bird, Armstrong & Hassager Reference Bird, Armstrong and Hassager1987; Larson Reference Larson1999; Morozov & Spagnolie Reference Morozov, Spagnolie and Spagnolie2015), including capillary flows (Eggers & Villermaux Reference Eggers and Villermaux2008). The stress contains contributions from both the polymers themselves, and a Newtonian contribution coming from the solvent. Its main simplifying features are that polymers are modelled with a single characteristic relaxation time of the polymer 
$\unicode[STIX]{x1D706}$, and that each polymer strand is assumed infinitely extensible, so that the extensional viscosity can become arbitrarily large. For that reason the Oldroyd-B model is the same as the so-called Hookean dumbbell model (Bird et al. Reference Bird, Armstrong and Hassager1987), where polymers are modelled by Hookean springs, with spherical beads attached at either end. We will measure the relaxation time relative to the capillary time 
$\unicode[STIX]{x1D70F}=\sqrt{\unicode[STIX]{x1D70C}R_{0}^{3}/\unicode[STIX]{x1D6FE}}$, where 
$\unicode[STIX]{x1D70C}$ is the density, 
$\unicode[STIX]{x1D6FE}$ the surface tension and 
$R_{0}$ the initial bridge radius; 
$De=\unicode[STIX]{x1D706}/\unicode[STIX]{x1D70F}$ is called the Deborah number.
In figure 1(b,c) we show the results of our own simulations, to be discussed in more detail below. After a short initial collapse of the liquid bridge, elastic stresses become important and delay pinch-off. In the case of a finite relaxation time (
$De=60$), the thread thickness decays exponentially in accordance with (1.1) (blue curve). We were able to follow this decay over more than one decade (from 
$h_{thr}/R_{0}\approx 0.26$ at the beginning of the exponential regime down to 
$h_{thr}/R_{0}\approx 0.014$), which takes place over 500 units relative to the capillary time scale. As we will see below, this means that the thread, and the transition region connecting it to the drop is well converged toward the asymptotic state we aim to investigate.
As the red curve in figure 1(b) we also show the thread radius for a simulation in which we have taken the limit of infinite relaxation time. As a result, elastic stresses do not relax, and a stationary thread thickness is reached at a point where elastic and surface tension forces balance. Below, we will show that the time-dependent problem of an exponentially thinning thread is closely related to the stationary elastic problem. Both are described by the same similarity solution, respectively, in the limit of long times, and of vanishing elasto-capillary number (given that contributions from the solvent are small).
The exponential law (1.1) was proposed in Entov (Reference Entov, Giesekus and Hibberd1988) and Bazilevskii et al. (Reference Bazilevskii, Entov, Rozhkov and Oliver1990). To calculate the prefactor 
$h_{0}$, one must take into account the initial conditions, in order to know the total stress that has built up inside the thread. However, even without following the history of the collapse, we will show on the basis of the full three-dimensional axisymmetric viscoelastic flow equations, but assuming that effects of the solvent viscosity are negligible, that
$$\begin{eqnarray}\unicode[STIX]{x1D70E}_{zz}=\frac{2\unicode[STIX]{x1D6FE}}{h_{thr}},\end{eqnarray}$$ where 
$\unicode[STIX]{x1D6FE}$ is the surface tension. Thus, by measuring the thread thickness, and also recording the extension rate 
$\dot{\unicode[STIX]{x1D716}}=-2{\dot{h}}_{thr}/h_{thr}$ of the flow inside the thread, one can infer the extensional viscosity 
$\unicode[STIX]{x1D702}_{E}=\unicode[STIX]{x1D70E}_{zz}/\dot{\unicode[STIX]{x1D716}}$, and thus use the device as a rheometer (Bazilevskii et al. Reference Bazilevskii, Entov, Rozhkov and Oliver1990; Amarouchene et al. Reference Amarouchene, Bonn, Meunier and Kellay2001).
Renardy (Reference Renardy1995) used lubrication theory to investigate the exponential thread thinning within the Oldroyd-B model. The lubrication description is valid as long as the interface slope is small, so the flow inside the thread can be described within a one-dimensional set of equations. In addition, Renardy (Reference Renardy1995) made pioneering use of Lagrangian variables, introduced into the one-dimensional description in Markovich & Renardy (Reference Markovich and Renardy1985). The thread was then matched to the neighbouring drop, where the slope of the interface increases exponentially, as we will show below. As a result, lubrication theory breaks down, and Renardy (Reference Renardy1995) obtained the incorrect result 
$\unicode[STIX]{x1D70E}_{zz}=\unicode[STIX]{x1D6FE}/(2h_{thr})$, which differs from (1.2) by a factor of 1/4.
Entov & Hinch (Reference Entov and Hinch1997) revisited the problem, but also allowed for a linear superposition of relaxation times instead of a single time 
$\unicode[STIX]{x1D706}$. By considering a force balance within the thread alone, and making an assumption about the tension inside the thread, Entov & Hinch (Reference Entov and Hinch1997) concluded that 
$\unicode[STIX]{x1D70E}_{zz}=\unicode[STIX]{x1D6FE}/h_{thr}$, which again differs from the correct result (1.2). Indeed, it was pointed out by Clasen et al. (Reference Clasen, Eggers, Fontelos, Li and McKinley2006) that a proper description of the thread cannot result from a force balance within the uniform thread alone (along which pressure gradients vanish), since the driving force of the motion results from the capillary pressure difference between the thread and the drops which connect to it (cf. sequence in figure 1a). To calculate the resistive force, one has to consider the flow of liquid emptying from the thread into the drops (of which in figure 1 there is one on either end of the thread). Thus, the task is to understand the transition region connecting the thread and the drop, which has the form of a similarity solution (Clasen et al. Reference Clasen, Eggers, Fontelos, Li and McKinley2006; Turkoz et al. Reference Turkoz, Lopez-Herrera, Eggers, Arnold and Deike2018): if both the axial and radial coordinates are rescaled by the thread thickness 
$h_{thr}$, the drop profile in the corner region connecting the thread and the drop falls onto a universal curve.
To deal with the transition toward the drop, Clasen et al. (Reference Clasen, Eggers, Fontelos, Li and McKinley2006) used an improved version of the lubrication approximation, obtained by keeping the full curvature term (Eggers Reference Eggers1993; Eggers & Dupont Reference Eggers and Dupont1994), so that surface tension contributions are accounted for correctly inside the drop. By solving the similarity equations in the limit of vanishing solvent viscosity contribution, it was shown that (1.1) holds within the improved lubrication approximation. It was also shown that this similarity solution is identical to the solution found previously by Entov & Yarin (Reference Entov and Yarin1984) for the collapse of a soft solid under surface tension. In fact, it is known that, at least in the one-dimensional approximation (Fontelos Reference Fontelos2003), in the limit of infinitely long relaxation times 
$\unicode[STIX]{x1D706}\rightarrow \infty$ the viscoelastic equations turn into the equations of large-deformation neo-Hookean elasticity.
However, even the improved lubrication theory cannot represent the flow inside the corner exactly if the slope is no longer small (see figure 1). It is also seen that the stress distribution has a significant radial dependence (cf. figure 1c), which cannot be captured by a one-dimensional model. It was therefore thought (Clasen et al. Reference Clasen, Eggers, Fontelos, Li and McKinley2006) that (1.2) could only be an approximation to the true relation between 
$\unicode[STIX]{x1D70E}_{zz}$ and 
$1/h_{thr}$, with a constant of proportionality different from 2. Here, we show that (1.2) can in fact be derived from the full three-dimensional equations without computing the corner similarity solution explicitly, making use of a new conserved quantity along the thread, also satisfied by the improved lubrication equation. As a consequence, the result of Clasen et al. (Reference Clasen, Eggers, Fontelos, Li and McKinley2006) is correct by serendipity.
Recently, Zhou & Doi (Reference Zhou and Doi2018) presented a ‘0-dimensional’ approximation of the thinning thread, representing the problem by a finite number of judiciously chosen scalar variables. The formation of threads, separated by drops, is modelled as a phase separation problem (Xuan & Biggins Reference Xuan and Biggins2017). Using a variational principle, Zhou & Doi (Reference Zhou and Doi2018) were able to reproduce the analytical result (1.2) of Clasen et al. (Reference Clasen, Eggers, Fontelos, Li and McKinley2006).
The present paper aims to close the gap between previous calculations based on various approximations, by instead solving the full three-dimensional axisymmetric Oldroyd-B equations, both numerically and theoretically. Since this model assumes a linear relationship between the stress and the extension of polymers, the exponential thinning regime described by (1.1) will continue forever, which allows us to focus on this critical part of the pinching dynamics. If finite extensibility of the polymer is taken into account (Fontelos & Li Reference Fontelos and Li2004; Renardy Reference Renardy2004; Wagner et al. Reference Wagner, Amarouchene, Bonn and Eggers2005), as in the so-called FENE-P model (Bird et al. Reference Bird, Armstrong and Hassager1987), pinching proceeds in a localised fashion, similar to the Newtonian case (Renardy Reference Renardy1995, Reference Renardy2002b; Fontelos & Li Reference Fontelos and Li2004). Other models that have been studied with view to including a more realistic, nonlinear response of the stress are the Giesekus model (Renardy Reference Renardy2001; Renardy & Losh Reference Renardy and Losh2002), the Phan-Thien–Tanner model (Renardy Reference Renardy2002a) or the generalised Newtonian fluid (Renardy Reference Renardy2002a); an overview of these calculations can be found in Renardy (Reference Renardy2004).
However, this mode of breakup does not appear to be a realistic description of experiment, as instead (Oliveira & McKinley Reference Oliveira and McKinley2005; Sattler, Wagner & Eggers Reference Sattler, Wagner and Eggers2008; Sattler et al. Reference Sattler, Gier, Eggers and Wagner2012) the polymer thread undergoes a spatially periodic instability, in the course of which the highly stressed thread partially relaxes, to form a series of small droplets. This has been called the blistering instability (Oliveira & McKinley Reference Oliveira and McKinley2005; Sattler et al. Reference Sattler, Wagner and Eggers2008, Reference Sattler, Gier, Eggers and Wagner2012), whose theoretical description (Eggers Reference Eggers2014) has recently been confirmed experimentally (Deblais et al. Reference Deblais, Velikov and Bonn2018b).
In the next section, we present the Oldroyd-B equations of motion. We use our own code, based on mapping the physical domain onto a simple rectangular domain (Herrada & Montanero Reference Herrada and Montanero2016; Deblais et al. Reference Deblais, Herrada, Hauner, Velikov, van Roon, Kellay, Eggers and Bonn2018a), together with a logarithmic transformation of the polymeric stress tensor (Fattal & Kupferman Reference Fattal and Kupferman2004; Coronado et al. Reference Coronado, Arora, Behr and Pasquali2007), in order to study the self-similar structure of the polymeric pinching problem in detail. We find that, contrary to previous assumptions, owing to the stress contribution from the solvent, the stress distribution inside the thread has a non-constant radial distribution. A similar non-trivial radial dependence of the stress has been discussed previously in Yao & McKinley (Reference Yao and McKinley1998). This suggests the similarity solution in general is non-universal, and dependent on boundary conditions. However, the contribution from the solvent is often small for realistic parameter values, as measured by a capillary number 
$\overline{v}_{0}$ based on the solvent viscosity, which we identify below. As a result, the main focus of our paper is on calculating the universal similarity solution in the limit of vanishing solvent contribution; the non-universal part comes from the solvent contribution.
In § 3 we present the similarity equations describing the corner region between the polymer thread and the drop including solvent contributions. Since inertia is found to be subdominant asymptotically, the tension in each cross-section is conserved. To solve the similarity equations for vanishing solvent parameter 
$\overline{v}_{0}$, in § 4 we consider a related problem in nonlinear elasticity: the collapse of a cylinder of elastic material under surface tension. As the elastic bridge deforms, elastic stresses build up and eventually balance surface tension to establish a new stationary equilibrium state. For small elastic modulus 
$G$, a thread of almost uniform thickness is formed, whose radius goes to zero in the limit of vanishing 
$G$. We show that this limit is described by a similarity solution, which is identical to the similarity solution for the time-dependent collapse of a viscoelastic fluid bridge.
Using Lagrangian coordinates with respect to the reference state of a cylindrical elastic bridge, we derive a simplified set of similarity equations for the elastic problem in the limit 
$G\rightarrow 0$. The elastic similarity equations are shown to obey a conservation law involving the pressure and the elastic energy, which is closely related to Eshelby’s elastic energy–momentum tensor (Eshelby Reference Eshelby1975). Physically, it represents invariance of the equations under a relabelling of Lagrangian coordinates.
With the help of this conservation law, we are able to derive (1.2), without having to solve the full three-dimensional similarity equations. However, the spatial structure of the similarity solution is different from what the earlier one-dimensional description predicts. We then solve the similarity equations numerically to yield the self-similar interface profile and the stress distribution in the interior. This then solves the problem both of the time-dependent pinching of a viscoelastic thread, and the collapse of an elastic thread.
In a final discussion, we compare numerical results to our similarity theory, and outline future directions of research.
2 Equations of motion and numerical simulations
The Oldroyd-B model can be described by the following set of equations (Bird et al. Reference Bird, Armstrong and Hassager1987; Morozov & Spagnolie Reference Morozov, Spagnolie and Spagnolie2015):
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}\left(\frac{\unicode[STIX]{x2202}\boldsymbol{v}}{\unicode[STIX]{x2202}t}+(\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{v}\right)=-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}_{p}+\unicode[STIX]{x1D702}_{s}\unicode[STIX]{x0394}\boldsymbol{v}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D748}_{p}}{\unicode[STIX]{x2202}t}+(\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\unicode[STIX]{x1D748}_{p}=(\unicode[STIX]{x1D735}\boldsymbol{v})^{T}\boldsymbol{\cdot }\unicode[STIX]{x1D748}_{p}+\unicode[STIX]{x1D748}_{p}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\boldsymbol{v})-\frac{\unicode[STIX]{x1D748}_{p}}{\unicode[STIX]{x1D706}}+\frac{\unicode[STIX]{x1D702}_{p}}{\unicode[STIX]{x1D706}}\dot{\unicode[STIX]{x1D738}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}=0, & \displaystyle\end{eqnarray}$$ where 
$\boldsymbol{v}$ is the velocity field, 
$\dot{\unicode[STIX]{x1D738}}=(\unicode[STIX]{x1D735}\boldsymbol{v})+(\unicode[STIX]{x1D735}\boldsymbol{v})^{T}$ the rate-of-deformation tensor and 
$\unicode[STIX]{x1D748}_{p}$ is the extra polymeric stress. Constants characterising the model are the fluid density 
$\unicode[STIX]{x1D70C}$, the polymeric viscosity 
$\unicode[STIX]{x1D702}_{p}$, the solvent viscosity 
$\unicode[STIX]{x1D702}_{s}$ and the polymer relaxation time 
$\unicode[STIX]{x1D706}$. As seen on the right of the momentum equation (2.1), the total stress on a fluid element is the sum of two contributions, coming from the solvent and from the polymer.
The Oldroyd-B model can be derived as the continuum version of a solution of non-interacting model polymers, each of which consists of two beads, experiencing Stokes’ drag, and connected by a Hookean spring. Thus on one hand, springs become stretched when beads move apart as described by the flow. On the other hand, the resulting tension in the spring contributes to the polymeric stress, described by the equation of motion (2.2); (2.3) enshrines incompressibility of the flow.
In the limit of small shear rates (so that polymers are hardly stretched at all), (2.1)–(2.3) describe a Newtonian fluid of total dynamic viscosity 
$\unicode[STIX]{x1D702}_{0}=\unicode[STIX]{x1D702}_{s}+\unicode[STIX]{x1D702}_{p}$. Polymeric stress relaxes at a rate 
$\unicode[STIX]{x1D706}$, as seen from the second to last term on the right of (2.2), while stretching by the flow is described by the first two terms on the right. In particular, in the case of an extensional flow with constant extension rate 
$\dot{\unicode[STIX]{x1D716}}$, the polymeric stress will grow exponentially if 
$\dot{\unicode[STIX]{x1D716}}\unicode[STIX]{x1D706}>1/2$. This follows from a comparison of the stretching and relaxation terms. We will see below that in fact 
$\dot{\unicode[STIX]{x1D716}}\unicode[STIX]{x1D706}=2/3$ inside the thread.
From a continuum perspective, the first four terms of (2.2) are known as the ‘upper convected derivative’
$$\begin{eqnarray}\overset{\triangledown }{\unicode[STIX]{x1D748}}_{p}=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D748}_{p}}{\unicode[STIX]{x2202}t}+(\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\unicode[STIX]{x1D748}_{p}-(\unicode[STIX]{x1D735}\boldsymbol{v})^{T}\unicode[STIX]{x1D748}_{p}-\unicode[STIX]{x1D748}_{p}(\unicode[STIX]{x1D735}\boldsymbol{v}).\end{eqnarray}$$Its form can be derived purely from the requirement that the polymeric stress ought to transform consistently as a covariant tensor (Morozov & Spagnolie Reference Morozov, Spagnolie and Spagnolie2015).
On account of the axisymmetry of the problem, the velocity field can be written 
$\boldsymbol{v}=v_{r}\boldsymbol{e}_{r}+v_{z}\boldsymbol{e}_{z}$ in cylindrical coordinates, and the polymeric stress tensor is 
$\unicode[STIX]{x1D748}_{p}=\unicode[STIX]{x1D70E}_{rr}\boldsymbol{e}_{r}\otimes \boldsymbol{e}_{r}+\unicode[STIX]{x1D70E}_{rz}\boldsymbol{e}_{r}\otimes \boldsymbol{e}_{z}+\unicode[STIX]{x1D70E}_{zz}\boldsymbol{e}_{z}\otimes \boldsymbol{e}_{z}+\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D703}\unicode[STIX]{x1D703}}\boldsymbol{e}_{\unicode[STIX]{x1D703}}\otimes \boldsymbol{e}_{\unicode[STIX]{x1D703}}$. The stress boundary condition at the free surface is
$$\begin{eqnarray}\boldsymbol{n}\boldsymbol{\cdot }(\unicode[STIX]{x1D748}_{p}+\unicode[STIX]{x1D702}_{s}\dot{\unicode[STIX]{x1D738}})=(p-\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D705})\boldsymbol{n},\end{eqnarray}$$where
$$\begin{eqnarray}\unicode[STIX]{x1D705}=\frac{1}{h(1+h_{z}^{2})^{1/2}}-\frac{h_{zz}}{(1+h_{z}^{2})^{3/2}},\quad \boldsymbol{n}=\frac{\boldsymbol{e}_{r}-\boldsymbol{e}_{z}h_{z}}{(1+h_{z}^{2})^{1/2}},\end{eqnarray}$$ are (twice) the mean curvature and the surface normal, respectively. If 
$h(z,t)$ is the thread profile, the kinematic boundary condition becomes
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}t}+v_{z}(z,h)\frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}z}=v_{r}(z,h).\end{eqnarray}$$Powerful free-surface codes for flows of Newtonian and viscoelastic fluids have been under development since the 1980s (Kistler & Scriven Reference Kistler and Scriven1984; Szady et al. Reference Szady, Salamon, Liu, Bornside, Armstrong and Brown1995; Bhat, Basaran & Pasquali Reference Bhat, Basaran and Pasquali2008; Mitsoulis & Tsamopoulos Reference Mitsoulis and Tsamopoulos2017), making ingenious use of the conformation tensor (Pasquali & Scriven Reference Pasquali and Scriven2002). Simulations of viscoelastic liquid bridges have been performed since the 1990s (Yao, McKinley & Debbaut Reference Yao, McKinley and Debbaut1998; Yao & McKinley Reference Yao and McKinley1998; Yao, Spiegelberg & McKinley Reference Yao, Spiegelberg and McKinley2000) to investigate filament stretching devices, and more recently to study capillary-driven viscoelastic flows (Etienne et al. Reference Etienne, Hinch and Li2006; Bhat et al. Reference Bhat, Appathurai, Harris, Pasquali, McKinley and Basaran2010; Morrison & Harlen Reference Morrison and Harlen2010; Turkoz et al. Reference Turkoz, Lopez-Herrera, Eggers, Arnold and Deike2018), the subject of the present paper.
However, we are not aware of a simulation which looked in detail into the transition region connecting the viscoelastic thread, where stresses are very large, and the drop, where stresses have relaxed, as illustrated by the colour map in figure 1(c). In order to study the asymptotics of the viscoelastic thread requires a regime where the thread is both much thinner than the drop (to ensure scale separation for a similarity analysis to be valid), and where elastic forces are strong compared to capillary forces, resulting in a elastic thread of uniform thickness. This means we aim to trace the exponential collapse of the thread for as long as possible (for more than a decade of perfectly exponential decay in the simulation analysed in this paper), at a Deborah number as large as possible. In our simulation (see figure 1) 
$De=60$, which means that exponential decay is observed over 500 dimensionless time units. In Turkoz et al. (Reference Turkoz, Lopez-Herrera, Eggers, Arnold and Deike2018), a similar benchmark was reached, but there remain problems with this numerical method, which prevent us from obtaining smooth stress profiles.
The numerical problem of solving (2.1)–(2.3) with boundary conditions (2.5) and (2.7) under conditions of breakup is highly demanding, for a number of reasons. First, the Oldroyd-B equation is well known to exhibit numerical instabilities (for reasons which are not well understood) if 
$\dot{\unicode[STIX]{x1D716}}\unicode[STIX]{x1D706}$ is of order unity (Fattal & Kupferman Reference Fattal and Kupferman2004). This is sometimes known as the high Weissenberg number problem (HWNP). Second, it is crucial to correctly describe the balance between surface tension forces and the force exerted by the polymer, both of which go to zero in the limit of vanishing thread thickness (the tension in the thread goes to zero, since it is the stress multiplied by the cross-sectional area of the thread).
2.1 Logarithmic transformation
In order to avoid the HWNP as 
$\dot{\unicode[STIX]{x1D716}}\unicode[STIX]{x1D706}$ is of order unity, a partial (two-dimensional) matrix-logarithm transformation of the polymeric stress tensor 
$\unicode[STIX]{x1D748}_{p}$ is applied in this work (Fattal & Kupferman Reference Fattal and Kupferman2004; Turkoz et al. Reference Turkoz, Lopez-Herrera, Eggers, Arnold and Deike2018). To this end we write the polymeric stress in terms of the conformation tensor 
$\unicode[STIX]{x1D63C}$ according to
$$\begin{eqnarray}\unicode[STIX]{x1D70E}_{zz}=\unicode[STIX]{x1D702}_{p}(\unicode[STIX]{x1D608}_{zz}-1)/\unicode[STIX]{x1D706},\quad \unicode[STIX]{x1D70E}_{rr}=\unicode[STIX]{x1D702}_{p}(\unicode[STIX]{x1D608}_{rr}-1)/\unicode[STIX]{x1D706},\quad \unicode[STIX]{x1D70E}_{rz}=\unicode[STIX]{x1D702}_{p}\unicode[STIX]{x1D608}_{rz}/\unicode[STIX]{x1D706};\end{eqnarray}$$ the equation of motion for 
$\unicode[STIX]{x1D63C}$ is
$$\begin{eqnarray}\overset{\triangledown }{\unicode[STIX]{x1D63C}}=-\frac{1}{\unicode[STIX]{x1D706}}(\unicode[STIX]{x1D63C}-\unicode[STIX]{x1D644}).\end{eqnarray}$$ Now, 
$\unicode[STIX]{x1D63C}$ is written as the logarithm of a symmetric, positive–definite matrix 
$\unicode[STIX]{x1D74D}$: 
$\unicode[STIX]{x1D63C}=\ln \unicode[STIX]{x1D74D}$. Equation (2.2), which allows us to advance the stresses 
$\unicode[STIX]{x1D70E}_{rr}$, 
$\unicode[STIX]{x1D70E}_{zz}$ and 
$\unicode[STIX]{x1D70E}_{rz}$ in time, is replaced with three equations for the elements of the 
$2\times 2$ matrix 
$\unicode[STIX]{x1D74D}$. The first step is to decompose this matrix as
$$\begin{eqnarray}\unicode[STIX]{x1D74D}=\left(\begin{array}{@{}cc@{}}\unicode[STIX]{x1D713}_{zz} & \unicode[STIX]{x1D713}_{rz}\\ \unicode[STIX]{x1D713}_{rz} & \unicode[STIX]{x1D713}_{rr}\end{array}\right)=\boldsymbol{R}\left(\begin{array}{@{}cc@{}}\ln (\unicode[STIX]{x1D6EC}_{1}) & 0\\ 0 & \ln (\unicode[STIX]{x1D6EC}_{2})\end{array}\right)\boldsymbol{R}^{T},\end{eqnarray}$$
$\ln (\unicode[STIX]{x1D6EC}_{1})$ and 
$\ln (\unicode[STIX]{x1D6EC}_{2})$ being the eigenvalues of 
$\unicode[STIX]{x1D74D}$, and the columns of 
$\boldsymbol{R}$ containing their normalised eigenvectors: 
$\boldsymbol{R}\boldsymbol{R}^{T}=\unicode[STIX]{x1D644}$. The second step is to construct the conformation tensor 
$\unicode[STIX]{x1D63C}$,
$$\begin{eqnarray}\unicode[STIX]{x1D63C}=\left(\begin{array}{@{}cc@{}}\unicode[STIX]{x1D608}_{zz} & \unicode[STIX]{x1D608}_{rz}\\ \unicode[STIX]{x1D608}_{rz} & \unicode[STIX]{x1D608}_{rr}\end{array}\right)=\boldsymbol{R}\left(\begin{array}{@{}cc@{}}\unicode[STIX]{x1D6EC}_{1} & 0\\ 0 & \unicode[STIX]{x1D6EC}_{2}\end{array}\right)\boldsymbol{R}^{T}.\end{eqnarray}$$ Now, the stresses 
$\unicode[STIX]{x1D70E}_{zz}$, 
$\unicode[STIX]{x1D70E}_{rz}$ and 
$\unicode[STIX]{x1D70E}_{rr}$ can be expressed as function of the elements of 
$\unicode[STIX]{x1D74D}$ with the help of (2.8).
Finally, the equation of motion for 
$\unicode[STIX]{x1D74D}$ is
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D74D}}{\unicode[STIX]{x2202}t}+(\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\unicode[STIX]{x1D74D}-(\unicode[STIX]{x1D734}\unicode[STIX]{x1D74D}-\unicode[STIX]{x1D74D}\unicode[STIX]{x1D734})-2\unicode[STIX]{x1D63D}=\frac{1}{\unicode[STIX]{x1D706}}[\unicode[STIX]{x1D63C}^{-1}-\unicode[STIX]{x1D644}],\end{eqnarray}$$ with matrices 
$\unicode[STIX]{x1D734}$, 
$\unicode[STIX]{x1D63D}$ given by
$$\begin{eqnarray}\unicode[STIX]{x1D734}=\boldsymbol{ R}\left(\begin{array}{@{}cc@{}}0 & \unicode[STIX]{x1D714}\\ -\unicode[STIX]{x1D714} & 0\end{array}\right)\boldsymbol{R}^{T},\quad \unicode[STIX]{x1D63D}=\boldsymbol{R}\left(\begin{array}{@{}cc@{}}\unicode[STIX]{x1D614}_{11} & 0\\ 0 & \unicode[STIX]{x1D614}_{22}\end{array}\right)\boldsymbol{R}^{T}.\end{eqnarray}$$Here
$$\begin{eqnarray}\unicode[STIX]{x1D648}=\left(\begin{array}{@{}cc@{}}\unicode[STIX]{x1D614}_{11} & \unicode[STIX]{x1D614}_{12}\\ \unicode[STIX]{x1D614}_{21} & \unicode[STIX]{x1D614}_{22}\end{array}\right)=\boldsymbol{R}^{T}(\unicode[STIX]{x1D735}\boldsymbol{v})_{2d}^{T}\boldsymbol{R},\quad (\unicode[STIX]{x1D735}\boldsymbol{v})_{2d}=\left(\begin{array}{@{}cc@{}}{\displaystyle \frac{\unicode[STIX]{x2202}v_{z}}{\unicode[STIX]{x2202}z}} & {\displaystyle \frac{\unicode[STIX]{x2202}v_{r}}{\unicode[STIX]{x2202}z}}\\[12.0pt] {\displaystyle \frac{\unicode[STIX]{x2202}v_{z}}{\unicode[STIX]{x2202}r}} & {\displaystyle \frac{\unicode[STIX]{x2202}v_{r}}{\unicode[STIX]{x2202}r}}\end{array}\right),\end{eqnarray}$$ and 
$\unicode[STIX]{x1D714}=[\unicode[STIX]{x1D6EC}_{2}\unicode[STIX]{x1D614}_{12}+\unicode[STIX]{x1D614}_{21}\unicode[STIX]{x1D6EC}_{1}]/[\unicode[STIX]{x1D6EC}_{2}-\unicode[STIX]{x1D6EC}_{1}]$. In the case that 
$\unicode[STIX]{x1D6EC}_{1}=\unicode[STIX]{x1D6EC}_{2}$, matrices 
$\unicode[STIX]{x1D734}$, 
$\unicode[STIX]{x1D63D}$ are replaced by
$$\begin{eqnarray}\unicode[STIX]{x1D734}=\mathbf{0},\quad \unicode[STIX]{x1D63D}={\textstyle \frac{1}{2}}[(\unicode[STIX]{x1D735}\boldsymbol{v})_{2d}^{T}+(\unicode[STIX]{x1D735}\boldsymbol{v})_{2d}].\end{eqnarray}$$ The primary quantity in this numerical scheme is 
$\unicode[STIX]{x1D74D}$, for which we solve (2.12). All other quantities, such as 
$\unicode[STIX]{x1D63C}$, are then derived from it.
2.2 Mapping technique
The numerical technique used in this study is a variation of that developed in Herrada & Montanero (Reference Herrada and Montanero2016). The spatial physical domain occupied by the fluid is mapped onto a rectangular domain 
$0\leqslant \unicode[STIX]{x1D702}\leqslant 1$, 
$0\leqslant \unicode[STIX]{x1D709}\leqslant L$ by means of the coordinate transformation 
$\unicode[STIX]{x1D702}=r/h(z,t)$ and 
$\unicode[STIX]{x1D709}=z$. It is worth mentioning that the mapping method can also be applied to situations where the free surface has overhangs (Ponce-Torres et al. Reference Ponce-Torres, Montanero, Herrada, Vega and Vega2017; Deblais et al. Reference Deblais, Herrada, Hauner, Velikov, van Roon, Kellay, Eggers and Bonn2018a). Now each variable (
$v_{z}$, 
$v_{r}$, 
$p$, 
$\unicode[STIX]{x1D713}_{rr}$, 
$\unicode[STIX]{x1D713}_{rz}$, 
$\unicode[STIX]{x1D713}_{zz}$, 
$\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D703}\unicode[STIX]{x1D703}}$ and 
$h$) and all its spatial and temporal derivatives, which appear in the transformed equations, are written as a single symbolic vector. For example, seven symbolic elements are associated with the axial velocity, 
$v_{z}$: 
$v_{z}$, 
$\unicode[STIX]{x2202}v_{z}/\unicode[STIX]{x2202}\unicode[STIX]{x1D702}$, 
$\unicode[STIX]{x2202}v_{z}/\unicode[STIX]{x2202}\unicode[STIX]{x1D709}$, 
$\unicode[STIX]{x2202}^{2}v_{z}/\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D702}$, 
$\unicode[STIX]{x2202}^{2}v_{z}/\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D709}$, 
$\unicode[STIX]{x2202}^{2}v_{z}/\unicode[STIX]{x2202}\unicode[STIX]{x1D702}\unicode[STIX]{x1D709}$ and 
$\unicode[STIX]{x2202}v_{z}/\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}$, 
$\unicode[STIX]{x1D70F}=t$ being the transformation for time. The next step is to use a symbolic toolbox to calculate the analytical Jacobians of all the equations with respect to the symbolic vector. Using these analytical Jacobians, we generate functions which can be evaluated later point by point once the transformed domain is discretised in space and time. In this work, we used the MATLAB tool matlabFunction to convert the symbolic Jacobians and equations in MATLAB functions.
Next, we carry out the temporal and spatial discretisation of the transformed domains. The spatial domain is discretised using a set of 
$n_{\unicode[STIX]{x1D709}}=850$ equally spaced collocation points in the axial (
$\unicode[STIX]{x1D709}$) direction, while a set of 
$n_{\unicode[STIX]{x1D702}}=20$ of equally spaced collocation points are used in the radial (
$\unicode[STIX]{x1D702}$) direction. Fourth-order finite differences are employed to compute the collocation matrices associated with the discretised points. In the time domain, second-order backward finite differences are used for the discretisation.
The final step is to solve at each time step the nonlinear system of discretised equations using a Newton procedure, where the numerical Jacobian is constructed using the spatial collocation matrices and the saved analytical functions. Since the method is fully implicit, a relatively large fixed time step, 
$\text{d}t/\unicode[STIX]{x1D70F}=0.1$, could be used in the simulation. Since 
$\unicode[STIX]{x1D706}$ sets a fixed time scale on which the thread evolves, the simulation remains resolved as the thread thins. We made sure that the results remained virtually the same when the step size was halved.
2.3 Numerical results and parameters
As a representative example, we simulate a liquid bridge with periodic boundary conditions, using the parameters of Turkoz et al. (Reference Turkoz, Lopez-Herrera, Eggers, Arnold and Deike2018). The initial condition is a cylinder of radius 
$R_{0}$, with a sinusoidal perturbation of small amplitude 
$\unicode[STIX]{x1D716}$ added to it,
$$\begin{eqnarray}\frac{h(z,0)}{R_{0}}=1-\unicode[STIX]{x1D716}\cos \left(\frac{z}{2R_{0}}\right),\quad -2\unicode[STIX]{x03C0}\leqslant \frac{z}{R_{0}}\leqslant 2\unicode[STIX]{x03C0},\end{eqnarray}$$ and 
$\unicode[STIX]{x1D716}=0.05$. Only half of the domain was simulated, using symmetry.
We define the dimensionless parameters of the simulation as in Clasen et al. (Reference Clasen, Eggers, Fontelos, Li and McKinley2006). First, the Deborah number 
$De=\unicode[STIX]{x1D706}/\unicode[STIX]{x1D70F}$ is the dimensionless relaxation time, compared to the capillary time scale 
$\unicode[STIX]{x1D70F}=\sqrt{\unicode[STIX]{x1D70C}R_{0}^{3}/\unicode[STIX]{x1D6FE}}$. The Ohnesorge number 
$Oh=\unicode[STIX]{x1D702}_{0}/\sqrt{R_{0}\unicode[STIX]{x1D70C}\unicode[STIX]{x1D6FE}}$ measures the relative importance of viscous and inertial effects, while 
$S=\unicode[STIX]{x1D702}_{s}/\unicode[STIX]{x1D702}_{0}$ is the solvent viscosity ratio. The ratio 
$G=\unicode[STIX]{x1D702}_{p}/\unicode[STIX]{x1D706}$ defines an elastic constant; it measures the effective shear modulus of the viscoelastic fluid on time scales shorter than the relaxation time 
$\unicode[STIX]{x1D706}$. Its dimensionless version
$$\begin{eqnarray}E_{c}=GR_{0}/\unicode[STIX]{x1D6FE}=\unicode[STIX]{x1D702}_{p}R_{0}/(\unicode[STIX]{x1D706}\unicode[STIX]{x1D6FE}),\end{eqnarray}$$is known as the elasto-capillary number.
In most of the plots shown below, we choose the same parameter values as those in Turkoz et al. (Reference Turkoz, Lopez-Herrera, Eggers, Arnold and Deike2018): 
$De=60$, 
$Oh=3.16$, 
$S=\unicode[STIX]{x1D702}_{s}/\unicode[STIX]{x1D702}_{0}=0.25$. The elasto-capillary number is 
$E_{c}=0.0395$. Our mapping technique allows for the description of all fields with high accuracy, avoiding the singular behaviour of the stress near the free surface reported previously (Turkoz et al. Reference Turkoz, Lopez-Herrera, Eggers, Arnold and Deike2018). The use of the logarithmic transformation described in § 2.1, together with a fully implicit formulation, allows for superior stability. As a result, we were able to follow the dynamical evolution over dimensionless times in excess of 
$t/\unicode[STIX]{x1D70F}=500$, while the thread radius is two orders of magnitude smaller than the original radius, and the exponential decay is observed over more than a decade. In Etienne et al. (Reference Etienne, Hinch and Li2006), 
$De$ is similar to ours, but exponential decay is followed over only 1 dimensionless time unit. In Bhat et al. (Reference Bhat, Appathurai, Harris, Pasquali, McKinley and Basaran2010), on the other hand, exponential behaviour is followed for a decade, but the Deborah number is only 0.1 in that case, and so the evolution is only followed for 0.6 time units.
3 Similarity description
3.1 Solution inside the thread
As shown in Clasen et al. (Reference Clasen, Eggers, Fontelos, Li and McKinley2006), and illustrated in figure 1, there are three different regions involved in this problem: first, a thread of uniform thickness, inside which polymers are highly stretched by an extensional flow of constant extension rate 
$\dot{\unicode[STIX]{x1D716}}$,
$$\begin{eqnarray}v_{r}(r,z,t)=-\dot{\unicode[STIX]{x1D716}}r/2,\quad v_{z}(r,z,t)=\dot{\unicode[STIX]{x1D716}}z.\end{eqnarray}$$As a result, using the kinematic boundary condition at the interface, the thread radius behaves as
$$\begin{eqnarray}h_{thr}=\overline{h}_{0}R_{0}e^{-\dot{\unicode[STIX]{x1D716}}t/2}\equiv \overline{h}_{0}\ell ,\quad \ell =R_{0}e^{-\dot{\unicode[STIX]{x1D716}}t/2},\end{eqnarray}$$ where 
$\overline{h}_{0}$ is a dimensionless constant. Second, fluid is injected from the thread into a drop, in which elastic stresses have relaxed almost completely, and the force balance is dominated by capillarity. Third, connecting the thread and the drop is a corner region, the typical size of which is set by the thread thickness, for which we aim to find a similarity description. A typical velocity scale is set by the velocity 
$v_{0}=\dot{\unicode[STIX]{x1D716}}L/2$ at the exit of the thread, based on the linear axial velocity field 
$v_{z}$ (cf. (3.1)), and a vanishing velocity at the middle of the thread (on account of symmetry), of total length 
$L$.
We begin by finding the solution inside the thread, where the radius is assumed uniform, shrinking at an exponential rate, consistent with the incompressible extensional velocity field (3.1). Owing to translational invariance along the thread, we assume that the polymeric stress is 
$z$-independent, but allow for an arbitrary dependence in the radial direction. We expect stresses to be of the same magnitude as the capillary pressure, which scales like 
$\unicode[STIX]{x1D6FE}/h_{thr}$. We thus introduce the scaling
$$\begin{eqnarray}\unicode[STIX]{x1D70E}_{zz}=\frac{\unicode[STIX]{x1D6FE}}{\ell }\overline{\unicode[STIX]{x1D70E}}_{0}(\overline{r})+C,\quad \overline{r}=\frac{r}{\ell },\end{eqnarray}$$ where 
$\overline{\unicode[STIX]{x1D70E}}_{0}$ is the dimensionless axial stress inside the thread. Inserting (3.3) into (2.2), the axial component of the terms on the left are,
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{zz}}{\unicode[STIX]{x2202}t}+v_{r}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{zz}}{\unicode[STIX]{x2202}r}=\frac{\dot{\unicode[STIX]{x1D716}}\unicode[STIX]{x1D6FE}}{2\ell }\overline{\unicode[STIX]{x1D70E}}_{0}+\frac{\dot{\unicode[STIX]{x1D716}}\unicode[STIX]{x1D6FE}}{2\ell }\overline{r}\overline{\unicode[STIX]{x1D70E}}_{0}^{\prime }-\frac{\dot{\unicode[STIX]{x1D716}}\unicode[STIX]{x1D6FE}}{2\ell }\overline{r}\overline{\unicode[STIX]{x1D70E}}_{0}^{\prime }=\frac{\dot{\unicode[STIX]{x1D716}}\unicode[STIX]{x1D6FE}}{2\ell }\overline{\unicode[STIX]{x1D70E}}_{0},\end{eqnarray}$$and thus the axial component of (2.2) becomes
$$\begin{eqnarray}\frac{\dot{\unicode[STIX]{x1D716}}\unicode[STIX]{x1D6FE}}{2\ell }\overline{\unicode[STIX]{x1D70E}}_{0}=\left(2\dot{\unicode[STIX]{x1D716}}-\frac{1}{\unicode[STIX]{x1D706}}\right)\left(\frac{\unicode[STIX]{x1D6FE}\overline{\unicode[STIX]{x1D70E}}_{0}}{\ell }+C\right)+\frac{2\dot{\unicode[STIX]{x1D716}}\unicode[STIX]{x1D702}_{p}}{\unicode[STIX]{x1D706}}.\end{eqnarray}$$ For this to be a solution, we have 
$\dot{\unicode[STIX]{x1D716}}=2/(3\unicode[STIX]{x1D706})$ (Bazilevskii et al. Reference Bazilevskii, Entov, Rozhkov and Oliver1990; Entov & Hinch Reference Entov and Hinch1997) and 
$C=-4\unicode[STIX]{x1D702}_{p}/\unicode[STIX]{x1D706}$, so that the solution in the thread becomes
$$\begin{eqnarray}\unicode[STIX]{x1D70E}_{zz}=\frac{\unicode[STIX]{x1D6FE}}{\ell }\overline{\unicode[STIX]{x1D70E}}_{0}(\overline{r})-\frac{4\unicode[STIX]{x1D702}_{p}}{\unicode[STIX]{x1D706}}.\end{eqnarray}$$ Here, 
$\overline{\unicode[STIX]{x1D70E}}_{0}(\overline{r})$ inside the thread can have an arbitrary radial profile, whereas usually a uniform profile has been assumed. The exact form of the profile will be determined by the initial (non-universal) dynamics of the thread. Analysing the other two components in a similar manner, we find to leading order
$$\begin{eqnarray}\unicode[STIX]{x1D70E}_{rz}=\frac{\unicode[STIX]{x1D6FE}}{R_{0}^{3}}\overline{\unicode[STIX]{x1D70E}}_{rz}^{(0)}(\overline{r})\ell ^{2},\quad \overline{\unicode[STIX]{x1D70E}}_{rr}=-\frac{2\unicode[STIX]{x1D702}_{p}}{5\unicode[STIX]{x1D706}}.\end{eqnarray}$$ Thus, inside the thread, the axial stress 
$\unicode[STIX]{x1D70E}_{zz}$ dominates over the remaining components.
3.2 Similarity solution for the corner region
We now describe the transition region at the end of the thread, assuming that the thread thickness 
$h_{thr}\propto \ell$ is the only length scale,
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}h(z,t)=\ell \overline{h}(\overline{z}),\quad p(z,r,t)={\displaystyle \frac{\unicode[STIX]{x1D6FE}}{\ell }}\overline{p}(\overline{z},\overline{r}),\quad \\ \boldsymbol{v}(z,r,t)=v_{0}\overline{\boldsymbol{v}}(\overline{z},\overline{r}),\quad \unicode[STIX]{x1D748}(z,r,t)={\displaystyle \frac{\unicode[STIX]{x1D6FE}}{\ell }}\overline{\unicode[STIX]{x1D748}}(\overline{z},\overline{r}),\end{array}\right\}\end{eqnarray}$$ where 
$\overline{z}=z/\ell$ and 
$\overline{r}=r/\ell$. Then, to leading order as 
$\ell \rightarrow 0$, the similarity equations become
$$\begin{eqnarray}\displaystyle & \displaystyle \overline{\unicode[STIX]{x1D735}}\overline{p}=\overline{\unicode[STIX]{x1D735}}\boldsymbol{\cdot }\overline{\unicode[STIX]{x1D748}}_{p}+\overline{v}_{0}\overline{\unicode[STIX]{x0394}}\overline{\boldsymbol{v}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle (\overline{\boldsymbol{v}}\boldsymbol{\cdot }\overline{\unicode[STIX]{x1D735}})\overline{\unicode[STIX]{x1D748}}_{p}=(\overline{\unicode[STIX]{x1D735}}\overline{\boldsymbol{v}})^{T}\boldsymbol{\cdot }\overline{\unicode[STIX]{x1D748}}_{p}+\overline{\unicode[STIX]{x1D748}}_{p}\boldsymbol{\cdot }(\overline{\unicode[STIX]{x1D735}}\overline{\boldsymbol{v}}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \overline{\unicode[STIX]{x1D735}}\boldsymbol{\cdot }\overline{\boldsymbol{v}}=0, & \displaystyle\end{eqnarray}$$ where 
$\overline{v}_{0}=v_{0}\unicode[STIX]{x1D702}_{s}/\unicode[STIX]{x1D6FE}=\unicode[STIX]{x1D702}_{s}L/(3\unicode[STIX]{x1D706}\unicode[STIX]{x1D6FE})$, and the stress boundary condition is,
$$\begin{eqnarray}\boldsymbol{n}\boldsymbol{\cdot }(\overline{\unicode[STIX]{x1D748}}_{p}+\overline{v}_{0}\overline{\dot{\unicode[STIX]{x1D738}}})=(\overline{p}-\overline{\unicode[STIX]{x1D705}})\boldsymbol{n}.\end{eqnarray}$$ The dimensionless quantity 
$\overline{v}_{0}$ is a capillary number based on the solvent viscosity. Remarkably, the solvent still appears in the leading-order balance even as polymeric stresses go to infinity, as already observed in Clasen et al. (Reference Clasen, Eggers, Fontelos, Li and McKinley2006).
The system (3.7)–(3.10) does not contain time, but there are matching conditions to be satisfied for 
$\overline{z}\rightarrow \pm \infty$. For 
$\overline{z}\rightarrow -\infty$, the solution has to tend toward the thread solution
$$\begin{eqnarray}\overline{\boldsymbol{v}}\rightarrow \boldsymbol{e}_{z},\quad \overline{h}\rightarrow \overline{h}_{0},\quad \overline{p}\rightarrow 1/\overline{h}_{0},\quad \overline{\unicode[STIX]{x1D748}}_{p}\rightarrow \overline{\unicode[STIX]{x1D70E}}_{0}(\overline{r})\boldsymbol{e}_{z}\otimes \boldsymbol{e}_{z}.\end{eqnarray}$$ The first condition says that, since the scale of the similarity solution 
$h_{thr}$ is much smaller than the length 
$L$ of the thread, the velocity of the similarity solution has to match the velocity at the exit of the thread. The other conditions correspond to thickness and the stress distribution inside the thread. To verify the similarity form (3.6) we present our numerical results in scaled form in figures 2 and 3. We indeed observe a collapse of the data and moreover the numerics confirm the matching conditions (3.11). To obtain another measure of the quality of convergence toward the similarity solution, we show two more critical quantities in figure 4. On the left, we demonstrate that the polymeric stress decays very rapidly as one moves from the thread into the drop. Over a distance of 10 in units of the thread thickness, the stress has decayed by five orders of magnitude. Owing to the dominance of elastic stresses, the thread radius has converged to a virtually constant value over its entire length.
Results of a numerical simulation with parameters as in figure 1, rescaled to similarity variables according to (3.6). The profiles for three different times collapse onto a time-independent similarity solution. Panels (a–d) are the components of the stress tensor, evaluated at the surface, while (e–h) show the components of the velocity, the pressure (again at the surface), as well as the surface profile itself. All profiles collapse onto the time-independent similarity solution defined by (3.6). The profiles (as well as all corresponding results below) have been interpolated with splines as appropriate for their representation as fourth-order polynomials inside our numerical code.
The axial stress 
$\overline{\unicode[STIX]{x1D70E}}_{zz}$ on the surface (a), and a magnified view of the thread radius 
$\overline{h}$ (b), as a function of 
$\overline{z}$, for the last dimensionless time 
$t/\unicode[STIX]{x1D70F}=550$. At 
$\overline{z}=5$, inside the drop, the stress has decayed by five orders of magnitude relative to its value inside the thread. The thread radius, shown on the right, only varies by a relative amount of 10-5 over its entire length.
A convergence study of the axial stress distribution inside the uniform thread, taken at the middle of the thread 
$z=0$. Panels (a,c,e) correspond to 
$S=\unicode[STIX]{x1D702}_{s}/\unicode[STIX]{x1D702}_{0}=0.25$ as in figure 1, and (b,d,f) to 
$S=0.7$; all other parameters are as in figure 1. Panels (a,b) demonstrate convergence toward a self-similar stress distribution, which remains non-constant over the radius, as the dimensionless exit velocity converges toward its asymptotic value (c,d). Panels (e,f) show the self-similar stress profiles 
$\overline{\unicode[STIX]{x1D70E}}_{0}(\overline{r})$. The shape of the profiles, scaled according to (3.6), depends on the solvent ratio 
$S$.
Note that the similarity equations (3.7)–(3.8) are invariant under the transformation 
$\overline{h}\rightarrow c\overline{h}$, 
$\overline{p}\rightarrow \overline{p}/c$, and 
$\overline{\unicode[STIX]{x1D748}}_{p}\rightarrow \overline{\unicode[STIX]{x1D748}}_{p}/c$, and hence the dimensionless thread thickness 
$\overline{h}_{0}$ cannot be calculated as part of the solution. However, the combination 
$\overline{h}\overline{\unicode[STIX]{x1D748}}_{p}$ is invariant, which means the axial stress 
$\unicode[STIX]{x1D70E}_{zz}$ can be calculated in terms of 
$\overline{h}_{0}$, as we will do below.
As implied by (3.4) and (3.11), the stress will in general not be constant across the fibre cross-section, corresponding to a stress distribution 
$\overline{\unicode[STIX]{x1D70E}}_{0}(\overline{r})$, which depends on the radius explicitly. A similar dependence of the radial stress profile on the initial evolution of the filament has been described by Yao & McKinley (Reference Yao and McKinley1998). The presence of a non-trivial radial stress distribution is confirmed in figure 5, although deviations from a constant stress are small, since the dimensionless exit velocity 
$\overline{v}_{0}=0.035$ is small. Indeed, we will see below that, for 
$\overline{v}_{0}=0$, the radial stress distribution is uniform. On the top we show the axial stress on the axis and on the surface, which converge to different values in the long-time limit. While for 
$S=0.25$ the stress is greater on the axis, for 
$S=0.7$ it is the other way around.
On the bottom we show the corresponding profiles across the thread for the two different values of the solvent viscosity. The two shapes are entirely different, and indicate a non-universal dependence of the stress distribution on the initial dynamics of the collapsing bridge.
3.3 Force balance
The similarity equation (3.7) contains no inertia, and hence by Newton’s third law there must be a constant tension force in the liquid bridge, independent of the axial position. As we will see below, this tension has a contribution from both the fluid stress, integrated over the cross-section, and from surface tension Clasen et al. (Reference Clasen, Eggers, Fontelos, Li and McKinley2006). We follow Eggers & Fontelos (Reference Eggers and Fontelos2005) and integrate (3.7), written in the form 
$\overline{\unicode[STIX]{x1D735}}\boldsymbol{\cdot }(\overline{\unicode[STIX]{x1D748}}_{p}+\overline{v}_{0}\overline{\dot{\unicode[STIX]{x1D738}}}-p\unicode[STIX]{x1D644})\equiv \overline{\unicode[STIX]{x1D735}}\boldsymbol{\cdot }\overline{\unicode[STIX]{x1D748}}=0$, over the fluid volume bounded by the planes 
$\overline{z}=\overline{z}_{\pm }$. Using the dynamic boundary conditions in the form 
$\overline{\boldsymbol{n}}\boldsymbol{\cdot }\overline{\unicode[STIX]{x1D748}}=-\overline{\unicode[STIX]{x1D705}}\overline{\boldsymbol{n}}$, we have
$$\begin{eqnarray}0=\oint _{S}\overline{\boldsymbol{n}}\boldsymbol{\cdot }\overline{\unicode[STIX]{x1D748}}\,\text{d}s=-\int _{O}\overline{\unicode[STIX]{x1D705}}\overline{\boldsymbol{n}}\,\text{d}s+\int _{C_{+}}\boldsymbol{e}_{z}\boldsymbol{\cdot }\overline{\unicode[STIX]{x1D748}}\,\text{d}s-\int _{C_{-}}\boldsymbol{e}_{z}\boldsymbol{\cdot }\overline{\unicode[STIX]{x1D748}}\,\text{d}s,\end{eqnarray}$$ where 
$S$ is the closed surface between the planes 
$\overline{z}=\overline{z}_{+}$ and 
$\overline{z}=\overline{z}_{-}$, denoted 
$C\pm$, as well as the jet surface 
$O$. According to Eggers & Fontelos (Reference Eggers and Fontelos2005), the integral over 
$O$ is
$$\begin{eqnarray}\int _{O}\overline{\unicode[STIX]{x1D705}}\overline{\boldsymbol{n}}\,\text{d}s=-\left.\frac{2\unicode[STIX]{x03C0}\overline{h}}{\sqrt{1+\overline{h}_{\overline{z}}^{2}}}\right|_{\overline{z}_{-}}^{\overline{z}_{+}}\boldsymbol{e}_{z},\end{eqnarray}$$ and the integral over the 
$z$-component of 
$C_{+}$ and 
$C_{-}$ can be converted to (using incompressibility to transform the viscous term)
$$\begin{eqnarray}\int _{C_{\pm }}\boldsymbol{e}_{z}\boldsymbol{\cdot }\overline{\unicode[STIX]{x1D748}}\boldsymbol{\cdot }\boldsymbol{e}_{z}\,\text{d}s=2\unicode[STIX]{x03C0}\left[\int _{0}^{\overline{h}}\overline{r}\overline{\unicode[STIX]{x1D70E}}_{zz}^{(p)}\,\text{d}\overline{r}-\int _{0}^{\overline{h}}\overline{r}\overline{p}\,\text{d}\overline{r}-2\overline{v}_{0}\overline{h}\overline{v}_{r}(\overline{h},\overline{z})\right]_{\overline{z}=\overline{z}_{\pm }}.\end{eqnarray}$$ Thus, taking the 
$z$-component of (3.12), the force balance in the 
$z$-direction becomes
$$\begin{eqnarray}0=2\unicode[STIX]{x03C0}\left[\frac{\overline{h}}{\sqrt{1+\overline{h}_{\overline{z}}^{2}}}+\int _{0}^{\overline{h}}\overline{r}\left(\overline{\unicode[STIX]{x1D70E}}_{zz}^{(p)}-\overline{p}\right)\,\text{d}\overline{r}-2\overline{v}_{0}\overline{h}\overline{v}_{r}(\overline{h},\overline{z})\right]_{\overline{z}_{-}}^{\overline{z}_{+}},\end{eqnarray}$$ which expresses the fact that the tension 
$\unicode[STIX]{x03C0}\overline{T}$ in the liquid thread is constant,
$$\begin{eqnarray}\overline{T}=\frac{2\overline{h}}{\sqrt{1+\overline{h}_{\overline{z}}^{2}}}+2\int _{0}^{\overline{h}}\overline{r}(\overline{\unicode[STIX]{x1D70E}}_{zz}^{(p)}-\overline{p})\,\text{d}\overline{r}-4\overline{v}_{0}\overline{h}\overline{v}_{r}(\overline{h},\overline{z}).\end{eqnarray}$$ First, we evaluate the tension inside the thread, where we can use (3.11) (and thus 
$\overline{v}_{r}=0$, 
$\overline{p}=1/\overline{h}_{0}$ and 
$\overline{\unicode[STIX]{x1D70E}}_{zz}^{(p)}=\overline{\unicode[STIX]{x1D70E}}_{0}(\overline{r})$) to find that
$$\begin{eqnarray}\overline{T}=2\int _{0}^{\overline{h}_{0}}\overline{r}\overline{\unicode[STIX]{x1D70E}}_{0}(\overline{r})\,\text{d}\overline{r}+\overline{h}_{0}.\end{eqnarray}$$Note that the solvent term is subdominant inside the thread, and the remaining contributions are from the axial polymeric stress and from a capillary term, both of which are positive quantities, as we will confirm below.
Entov & Hinch (Reference Entov and Hinch1997) attempted to develop a theory of the polymeric pinching by considering a force balance over the uniform, cylindrical part alone, a strategy pursued further in a number of subsequent publications (McKinley & Tripathi Reference McKinley and Tripathi2000; Anna & McKinley Reference Anna and McKinley2001; Torres et al. Reference Torres, Hallmark, Wilson and Hilliou2014; Wagner, Bourouiba & McKinley Reference Wagner, Bourouiba and McKinley2015). Entov & Hinch (Reference Entov and Hinch1997) argued that the tension can be taken to vanish, because ‘the filament is attached to large stagnant drops on stationary end plates’. Clearly, this cannot be a correct approximation, since the right-hand side of (3.14) is strictly positive, and thus the tension must be strictly positive. This means that, for an isolated, inertialess drop (for which the tension would be zero, since no external forces are acting on the drop (Eggers & Fontelos Reference Eggers and Fontelos2005)), (3.14) could never be satisfied. Hence, such an isolated polymeric drop cannot form a pinching thread, just like an isolated drop of Newtonian fluid cannot break up (Eggers & Fontelos Reference Eggers and Fontelos2005).
The problem with Entov & Hinch’s (Reference Entov and Hinch1997) argument is that the physical thread tension 
$T=\unicode[STIX]{x1D6FE}\unicode[STIX]{x03C0}\overline{T}\ell$ goes to zero as the thread thins, so that the tension that needs to be supported by the end drops also vanishes as they relax toward a steady state. A small and vanishing deviation from a purely spherical shape is enough to provide a consistent balance. In fact, the assumption actually adopted by Entov & Hinch (Reference Entov and Hinch1997) is to take 
$\overline{T}=2\overline{h}_{0}$, so that 
$\overline{T}$ cancels with the first term on the right of (3.13) within the thread; as shown in (4.30), the correct expression is 
$\overline{T}=3\overline{h}_{0}$. Assuming a constant polymeric stress 
$\overline{\unicode[STIX]{x1D70E}}_{0}$ over the cross-section of the thread, (3.14) becomes 
$\overline{T}=\overline{h}_{0}^{2}\overline{\unicode[STIX]{x1D70E}}_{0}+\overline{h}_{0}$; combining with 
$\overline{T}=2\overline{h}_{0}$, one finds 
$\overline{\unicode[STIX]{x1D70E}}_{0}=1/\overline{h}_{0}$, which differs by a factor of 2 from the correct expression (1.2), as stated in the Introduction.
To bring (3.13) into a more convenient form for the task of calculating the tension inside the drop we note that
$$\begin{eqnarray}\frac{2\overline{h}}{\sqrt{1+\overline{h}_{\overline{z}}^{2}}}-2\int _{0}^{\overline{h}}\overline{r}\overline{\unicode[STIX]{x1D705}}\,\text{d}\overline{r}=\frac{2\overline{h}}{\sqrt{1+\overline{h}_{\overline{z}}^{2}}}-\overline{h}^{2}\overline{\unicode[STIX]{x1D705}}=\overline{h}^{2}\overline{K},\end{eqnarray}$$where
$$\begin{eqnarray}K=\frac{1}{h(1+h_{z}^{2})^{1/2}}+\frac{h_{zz}}{(1+h_{z}^{2})^{3/2}},\end{eqnarray}$$ and correspondingly for the rescaled quantity 
$\overline{K}$. As a result, (3.13) can be written as
$$\begin{eqnarray}\overline{T}\equiv \left.2\int _{0}^{\overline{h}_{0}}\overline{r}\overline{\unicode[STIX]{x1D70E}}_{0}(\overline{r})\,\text{d}\overline{r}+\overline{h}_{0}=2\int _{0}^{\overline{h}}\overline{r}\overline{\unicode[STIX]{x1D70E}}_{zz}^{(p)}\,\text{d}\overline{r}+2\int _{0}^{\overline{h}}\overline{r}(\overline{\unicode[STIX]{x1D705}}-\overline{p})\,\text{d}\overline{r}-4\overline{v}_{0}\overline{h}\overline{v}_{r}(\overline{h},\overline{z})+\overline{h}^{2}\overline{K}\right|_{\overline{z}=\overline{z}_{+}}.\end{eqnarray}$$ As seen in figures 2–4 and discussed in more detail below, all viscous and elastic stresses decay rapidly inside the drop, and the balance is maintained by surface tension forces alone. As a result, the polymeric stress component 
$\overline{\unicode[STIX]{x1D70E}}_{zz}^{(p)}$ as well as the velocity 
$\overline{v}_{r}$ decay rapidly as 
$\overline{z}\rightarrow \infty$, and the pressure almost equals the capillary pressure: 
$\overline{p}\approx \overline{\unicode[STIX]{x1D705}}$. Only the fourth term on the right of (3.16), coming from surface tension, survives, and equation describing the shape of the capillary region is (Eggers & Fontelos Reference Eggers and Fontelos2015),
$$\begin{eqnarray}\overline{T}=\overline{h}^{2}\overline{K}.\end{eqnarray}$$ This equation describes the cross-over toward the drop, whose dimensionless radius 
$\overline{R}_{d}=R_{d}/h_{thr}$ diverges in the limit. Indeed, in the intermediate region where 
$\overline{h}\ll \overline{R}_{d}$, the only way the right-hand side of (3.17) can be constant is to require that 
$\overline{h}_{\overline{z}}/\overline{h}\rightarrow \text{const.}$, which is achieved for
$$\begin{eqnarray}\overline{h}=He^{a\overline{z}},\quad \overline{z}\rightarrow \infty .\end{eqnarray}$$ This means the right-hand side of (3.17) becomes 
$2/a$ in the limit, and we obtain
$$\begin{eqnarray}\overline{T}=\frac{2}{a}.\end{eqnarray}$$4 The elastic correspondence
At least in the limit that the contribution from the solvent in the momentum balance (3.7) is negligible, one can make use of the fact that, in the similarity equations (3.7)–(3.8), terms containing the relaxation time 
$\unicode[STIX]{x1D706}$ have disappeared. Thus, we anticipate that the same set of equations can be obtained in a purely elastic description, where we take the limit 
$\unicode[STIX]{x1D706}\rightarrow \infty$ such that elastic stresses do not relax. In the absence of inertia, the elastic equations describe the equilibrium between surface tension, which tends to deform the elastic medium, and elasticity, which resists deformation. Thus for 
$\unicode[STIX]{x1D706}\rightarrow \infty$, but keeping 
$G=\unicode[STIX]{x1D702}_{p}/\unicode[STIX]{x1D706}$ fixed, this equilibrium state is described by
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}p=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}_{p}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \overset{\triangledown }{\unicode[STIX]{x1D748}}_{p}=G\dot{\unicode[STIX]{x1D738}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}=0, & \displaystyle\end{eqnarray}$$ using that 
$\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D748}_{p}=0$. One can already see that, in the limit of vanishing elasticity, 
$G\rightarrow 0$, (4.1)–(4.3) are identical to the similarity equation (3.7)–(3.8) for the case where 
$v_{0}=0$ (the stress coming from the solvent can be neglected). So indeed, the similarity solution can be computed from the elastic limit. However, it is advantageous to first retain the elastic term 
$G\dot{\unicode[STIX]{x1D738}}$, as it will allow us to properly introduce the history dependence into the description, and to determine the value of 
$\bar{h}_{0}$ we are looking for.
To proceed, we first show (see Larson (Reference Larson1988) and Snoeijer et al. (Reference Snoeijer, Pandey, Herrada and Eggers2019) for more details), that, instead of solving (4.2) with the constraint (4.3), one can write 
$\unicode[STIX]{x1D748}_{p}$ as the stress tensor of a neo-Hookean elastic solid, which depends quadratically on the deformation from an unstressed reference state. Thus one needs to find a transformation of the undeformed state (e.g. a cylinder), parameterised by variables 
$R,Z$, into a deformed state 
$r=r(R,Z),z=z(R,Z)$ such that elastic and surface tension forces are balanced. The equilibrium state is characterised by the deformation tensor (in Cartesian coordinates)
$$\begin{eqnarray}F_{iK}=\frac{\unicode[STIX]{x2202}x_{i}}{\unicode[STIX]{x2202}X_{K}},\end{eqnarray}$$which satisfies incompressibility
$$\begin{eqnarray}\det \boldsymbol{F}=\frac{r}{R}\left(\frac{\unicode[STIX]{x2202}r}{\unicode[STIX]{x2202}R}\frac{\unicode[STIX]{x2202}z}{\unicode[STIX]{x2202}Z}-\frac{\unicode[STIX]{x2202}r}{\unicode[STIX]{x2202}Z}\frac{\unicode[STIX]{x2202}z}{\unicode[STIX]{x2202}R}\right)=1.\end{eqnarray}$$Here, lower case variables and lower case indices refer to the deformed state, upper case variables and upper case indices refer to the undeformed state.
If we write the stress tensor as 
$\unicode[STIX]{x1D748}_{p}=G(\unicode[STIX]{x1D63C}-\unicode[STIX]{x1D644})$, where 
$\unicode[STIX]{x1D63C}$ is the conformation tensor as in (2.8), then (4.2) simplifies to 
$\overset{\triangledown }{\unicode[STIX]{x1D63C}}=0$. Using the fact that 
$\boldsymbol{v}=\text{d}\boldsymbol{x}/\text{d}t$, where the time derivative is taken at a constant material point 
$\boldsymbol{X}$, one shows that (Morozov & Spagnolie Reference Morozov, Spagnolie and Spagnolie2015)
$$\begin{eqnarray}\overset{\triangledown }{\unicode[STIX]{x1D63C}}=\boldsymbol{F}\left[\frac{\text{d}}{\text{d}t}(\boldsymbol{F}^{-1}\unicode[STIX]{x1D63C}(\boldsymbol{F}^{T})^{-1})\right]\boldsymbol{F}^{T}.\end{eqnarray}$$ From this relation, it is clear that 
$\unicode[STIX]{x1D63C}=\boldsymbol{F}\boldsymbol{F}^{T}$ indeed provides an integral to the equation 
$\overset{\triangledown }{\unicode[STIX]{x1D63C}}=0$. In addition, this solution satisfies the condition that 
$\unicode[STIX]{x1D748}_{p}=0$ in the initial configuration, where 
$\boldsymbol{F}=\unicode[STIX]{x1D644}$.
With the above steps we have managed to integrate (4.2), and have identified 
$\unicode[STIX]{x1D748}_{p}$, satisfying a stress-free initial condition. Since 
$-G\unicode[STIX]{x1D644}$ can be written as part of the pressure 
$p$, the elastic stress tensor becomes
$$\begin{eqnarray}\unicode[STIX]{x1D748}=-p\unicode[STIX]{x1D644}+\unicode[STIX]{x1D748}_{p}\equiv -p\unicode[STIX]{x1D644}+G\boldsymbol{F}\boldsymbol{F}^{T},\end{eqnarray}$$which is the classical result for the ‘true’ stress tensor of a neo-Hookean solid (Suo Reference Suo2013a,Reference Suob). Then, in terms of the deformation, the stress tensor becomes (Negahban Reference Negahban2012)
$$\begin{eqnarray}\unicode[STIX]{x1D748}_{p}=G\left(\begin{array}{@{}ccc@{}}\left({\displaystyle \frac{\unicode[STIX]{x2202}r}{\unicode[STIX]{x2202}R}}\right)^{2}+\left({\displaystyle \frac{\unicode[STIX]{x2202}r}{\unicode[STIX]{x2202}Z}}\right)^{2} & 0 & {\displaystyle \frac{\unicode[STIX]{x2202}r}{\unicode[STIX]{x2202}R}}{\displaystyle \frac{\unicode[STIX]{x2202}z}{\unicode[STIX]{x2202}R}}+{\displaystyle \frac{\unicode[STIX]{x2202}r}{\unicode[STIX]{x2202}Z}}{\displaystyle \frac{\unicode[STIX]{x2202}z}{\unicode[STIX]{x2202}Z}}\\ 0 & \left({\displaystyle \frac{r}{R}}\right)^{2} & 0\\ {\displaystyle \frac{\unicode[STIX]{x2202}r}{\unicode[STIX]{x2202}R}}{\displaystyle \frac{\unicode[STIX]{x2202}z}{\unicode[STIX]{x2202}R}}+{\displaystyle \frac{\unicode[STIX]{x2202}r}{\unicode[STIX]{x2202}Z}}{\displaystyle \frac{\unicode[STIX]{x2202}z}{\unicode[STIX]{x2202}Z}} & 0 & \left({\displaystyle \frac{\unicode[STIX]{x2202}z}{\unicode[STIX]{x2202}Z}}\right)^{2}+\left({\displaystyle \frac{\unicode[STIX]{x2202}z}{\unicode[STIX]{x2202}R}}\right)^{2}\end{array}\right).\end{eqnarray}$$ The true elastic stress satisfies the same equilibrium conditions as those of a fluid (Suo Reference Suo2013b) (but with 
$v_{s}=0$),
$$\begin{eqnarray}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}=0,\quad \left.\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D748}_{p}\right|_{S}=\unicode[STIX]{x1D6FE}(p-\unicode[STIX]{x1D705})\boldsymbol{n},\end{eqnarray}$$ where derivatives are to be taken with respect to current coordinates. If 
$R_{0}$ is the radius of the unperturbed interface as usual, the free surface 
$h(z)$ has the parametric representation
$$\begin{eqnarray}h(z(R_{0},Z))=r(R_{0},Z).\end{eqnarray}$$Transforming the equilibrium condition (4.9) to Lagrangian coordinates, we obtain in cylindrical coordinates,
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{GR}{r}\left[\frac{\unicode[STIX]{x2202}^{2}r}{\unicode[STIX]{x2202}R^{2}}+\frac{1}{R}\frac{\unicode[STIX]{x2202}r}{\unicode[STIX]{x2202}R}-\frac{r}{R^{2}}+\frac{\unicode[STIX]{x2202}^{2}r}{\unicode[STIX]{x2202}Z^{2}}\right]=\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}R}\frac{\unicode[STIX]{x2202}z}{\unicode[STIX]{x2202}Z}-\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}Z}\frac{\unicode[STIX]{x2202}z}{\unicode[STIX]{x2202}R}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{GR}{r}\left[\frac{\unicode[STIX]{x2202}^{2}z}{\unicode[STIX]{x2202}R^{2}}+\frac{1}{R}\frac{\unicode[STIX]{x2202}z}{\unicode[STIX]{x2202}R}+\frac{\unicode[STIX]{x2202}^{2}z}{\unicode[STIX]{x2202}Z^{2}}\right]=\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}Z}\frac{\unicode[STIX]{x2202}r}{\unicode[STIX]{x2202}R}-\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}R}\frac{\unicode[STIX]{x2202}r}{\unicode[STIX]{x2202}Z}, & \displaystyle\end{eqnarray}$$for the radial and axial force balances, respectively. The stress boundary conditions on the surface are
$$\begin{eqnarray}\displaystyle & \displaystyle -\left.\frac{h_{z}^{2}\unicode[STIX]{x1D70E}_{zz}-2h_{z}\unicode[STIX]{x1D70E}_{rz}+\unicode[STIX]{x1D70E}_{rr}}{1+h_{z}^{2}}=\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D705}-p)\right|_{r=h}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \left.h_{z}(\unicode[STIX]{x1D70E}_{rr}-\unicode[STIX]{x1D70E}_{zz})+(1-h_{z}^{2})\unicode[STIX]{x1D70E}_{rz}=0\right|_{r=h}. & \displaystyle\end{eqnarray}$$4.1 Elastic similarity equations
We now analyse the elastic problem in the limit 
$E_{c}=GR_{0}/\unicode[STIX]{x1D6FE}\rightarrow 0$, which corresponds to the similarity solution for the pinching of a polymeric thread. We expect a static balance between elastic stresses and surface tension to be established at a thread thickness which scales like an elasto-capillary length scale 
$\ell _{e}$: 
$z\propto r\propto \ell _{e}$. Here, we have anticipated that, as in the similarity solution (3.6), both Eulerian coordinates scale in the same way. By contrast, the Lagrangian coordinates do not have the same scales for 
$R$ and 
$Z$. Namely, the radial coordinate 
$R$ must be rescaled by the original radius 
$R_{0}$; on account of incompressibility (4.5) we have that 
$Z\propto \ell _{e}^{3}/R_{0}^{2}$. This implies that 
$\text{d}/\text{d}Z\gg \text{d}/\text{d}R$, so that the dominant stress contributions in (4.8) is of the order of 
$\unicode[STIX]{x1D748}_{p}\propto GR_{0}^{4}/\ell _{e}^{4}$. The curvature 
$\unicode[STIX]{x1D705}$ scales as the inverse radius 
$h^{-1}\propto \ell _{e}^{-1}$ of the thread (cf. (2.6)), and thus on account of the stress balance (4.9) 
$GR_{0}^{4}/\ell _{e}^{4}\propto \unicode[STIX]{x1D6FE}/\ell _{e}$. As a result, the elasto-capillary length scale becomes
$$\begin{eqnarray}\ell _{e}=\left(\frac{GR_{0}^{4}}{\unicode[STIX]{x1D6FE}}\right)^{1/3}=R_{0}E_{c}^{1/3},\end{eqnarray}$$which sets the thickness of the thread (Eggers & Fontelos Reference Eggers and Fontelos2015), at which elastic and capillary forces are balanced.
According to the scaling analysis above, we introduce the similarity solution
$$\begin{eqnarray}z=\ell _{e}\unicode[STIX]{x1D719}(\overline{R},\overline{Z}),\quad r=\ell _{e}\unicode[STIX]{x1D713}(\overline{R},\overline{Z}),\quad p=(\unicode[STIX]{x1D6FE}/\ell _{e})\unicode[STIX]{x1D6F1}(\overline{R},\overline{Z}),\end{eqnarray}$$ valid in the limit 
$E_{c}\rightarrow 0$, where 
$\overline{R}=R/R_{0}$ and 
$\overline{Z}=ZR_{0}^{2}/\ell _{e}^{3}$. The rescaled thread thickness 
$\overline{h}$ is defined by 
$h=\ell _{e}\overline{h}$, and 
$z=\ell _{e}\overline{z}$, so the slope 
$h_{z}=\overline{h}_{\overline{z}}$ remains invariant under the scalings, and 
$h_{zz}=\ell _{e}^{-1}\overline{h}_{\overline{z}\overline{z}}$. As a result, the curvature scales as 
$\unicode[STIX]{x1D705}=\ell _{e}^{-1}\overline{\unicode[STIX]{x1D705}}$, where 
$\overline{\unicode[STIX]{x1D705}}$ is defined as in (2.6), but on the basis of 
$\overline{h}(\overline{z})$. Note that, with this scaling, 
$\unicode[STIX]{x1D748}_{p}\propto GE_{c}^{-4/3}$, so that the right-hand side of (4.2) becomes subdominant in the limit. As a result, the elastic similarity equations defined by the above rescaling become identical to (3.7)–(3.8) in the limit that solvent effects are negligible (
$v_{0}=0$), and the similarity profiles (overlined variables defined above) are the similarity profiles defined by (3.6) in that limit. The only way 
$G$ still enters into the problem is through the initial condition, for which we require that the conformation tensor 
$\unicode[STIX]{x1D63C}=\unicode[STIX]{x1D644}$. The rescaling (3.6), on the other hand, still contains a free length scale 
$R_{0}$, which we choose such that the similarity profiles of both the time-dependent problem and of the elastic problem are the same.
The key observation is that the ratio 
$E_{c}=(\ell _{e}/R_{0})^{3}$ between a characteristic axial scale and a radial scale in Lagrangian coordinates becomes very small in the limit. As a result, radial derivatives can be neglected relative to axial ones. It also means that, in the limit 
$E_{c}\rightarrow 0$, fluid elements become stretched out in the axial direction, so that the self-similar region comes from just a single point along the 
$Z$-axis in the reference configuration. The radius 
$R_{0}$ has to be taken as the radius at that point in the reference configuration.
With the scalings (4.16), the similarity equations become, beginning with incompressibility,
$$\begin{eqnarray}\frac{\unicode[STIX]{x1D713}}{\overline{R}}(\unicode[STIX]{x1D713}_{\overline{R}}\unicode[STIX]{x1D719}_{\overline{Z}}-\unicode[STIX]{x1D713}_{\overline{Z}}\unicode[STIX]{x1D719}_{\overline{R}})=1.\end{eqnarray}$$ The equilibrium conditions (4.11) and (4.12) are, to leading order as 
$E_{c}\rightarrow 0$,
$$\begin{eqnarray}\frac{\overline{R}}{\unicode[STIX]{x1D713}}\unicode[STIX]{x1D713}_{\overline{Z}\overline{Z}}=\unicode[STIX]{x1D6F1}_{\overline{R}}\unicode[STIX]{x1D719}_{\overline{Z}}-\unicode[STIX]{x1D6F1}_{\overline{Z}}\unicode[STIX]{x1D719}_{\overline{R}},\quad \frac{\overline{R}}{\unicode[STIX]{x1D713}}\unicode[STIX]{x1D719}_{\overline{Z}\overline{Z}}=\unicode[STIX]{x1D6F1}_{\overline{Z}}\unicode[STIX]{x1D713}_{\overline{R}}-\unicode[STIX]{x1D6F1}_{\overline{R}}\unicode[STIX]{x1D713}_{\overline{Z}}.\end{eqnarray}$$The key point is that the second radial derivatives on the right-hand side have disappeared, and as result the equations are of first order only in the radial coordinates.
The extra stresses are at leading order
$$\begin{eqnarray}\unicode[STIX]{x1D70E}_{zz}=\frac{\unicode[STIX]{x1D6FE}}{\ell _{e}}\unicode[STIX]{x1D719}_{\overline{Z}}^{2},\quad \unicode[STIX]{x1D70E}_{rz}=\frac{\unicode[STIX]{x1D6FE}}{\ell _{e}}\unicode[STIX]{x1D719}_{\overline{Z}}\unicode[STIX]{x1D713}_{\overline{Z}},\quad \unicode[STIX]{x1D70E}_{rr}=\frac{\unicode[STIX]{x1D6FE}}{\ell _{e}}\unicode[STIX]{x1D713}_{\overline{Z}}^{2},\end{eqnarray}$$ and so the boundary conditions (4.13), (4.14), to be evaluated at the surface 
$\overline{R}=1$, become to leading order
$$\begin{eqnarray}\displaystyle & \displaystyle -\frac{h_{z}^{2}\unicode[STIX]{x1D719}_{\overline{Z}}^{2}-2h_{z}\unicode[STIX]{x1D719}_{\overline{Z}}\unicode[STIX]{x1D713}_{\overline{Z}}+\unicode[STIX]{x1D713}_{\overline{Z}}^{2}}{1+h_{z}^{2}}=\overline{\unicode[STIX]{x1D705}}-\unicode[STIX]{x1D6F1}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle h_{z}\left(\unicode[STIX]{x1D713}_{\overline{Z}}^{2}-\unicode[STIX]{x1D719}_{\overline{Z}}^{2}\right)+(1-h_{z}^{2})\unicode[STIX]{x1D719}_{\overline{Z}}\unicode[STIX]{x1D713}_{\overline{Z}}=0. & \displaystyle\end{eqnarray}$$But along the surface we have
$$\begin{eqnarray}\left.\frac{\unicode[STIX]{x1D713}_{\overline{Z}}}{\unicode[STIX]{x1D719}_{\overline{Z}}}\right|_{R=1}=\left.\frac{\unicode[STIX]{x2202}r}{\unicode[STIX]{x2202}z}\right|_{R=1}=h_{z},\end{eqnarray}$$ and so dividing through (4.20) and (4.21) by 
$\unicode[STIX]{x1D719}_{\overline{Z}}^{2}$, the left-hand side of both equations is seen to vanish. As a result, the tangential stress balance (4.21) is satisfied identically in the limit 
$E_{c}\rightarrow 0$, while the kinematic boundary condition as well as the normal stress balance become
$$\begin{eqnarray}\unicode[STIX]{x1D713}(1,\overline{Z})=\overline{h}(\unicode[STIX]{x1D719}(1,\overline{Z})),\quad \unicode[STIX]{x1D6F1}(1,\overline{Z})=\overline{\unicode[STIX]{x1D705}}.\end{eqnarray}$$ Thus, as (4.18) is only of first order in 
$\overline{R}$, there is only one stress boundary condition to be satisfied instead of two.
We now show that the pressure may be eliminated from the leading-order equations. Multiplying the first equation of (4.18) by 
$\unicode[STIX]{x1D719}_{\overline{Z}}$ and using (4.17), we find
$$\begin{eqnarray}\unicode[STIX]{x1D719}_{\overline{Z}}\unicode[STIX]{x1D719}_{\overline{Z}\overline{Z}}=\unicode[STIX]{x1D6F1}_{\overline{Z}}+\frac{\unicode[STIX]{x1D713}}{R}\unicode[STIX]{x1D713}_{\overline{Z}}\left(\unicode[STIX]{x1D6F1}_{\overline{Z}}\unicode[STIX]{x1D719}_{R}-\unicode[STIX]{x1D6F1}_{R}\unicode[STIX]{x1D719}_{\overline{Z}}\right)=\unicode[STIX]{x1D6F1}_{\overline{Z}}-\unicode[STIX]{x1D713}_{\overline{Z}}\unicode[STIX]{x1D713}_{\overline{Z}\overline{Z}},\end{eqnarray}$$ using the second (4.18) in the second step. This can be integrated over 
$\overline{Z}$ to
$$\begin{eqnarray}{\textstyle \frac{1}{2}}\left(\unicode[STIX]{x1D719}_{\overline{Z}}^{2}+\unicode[STIX]{x1D713}_{\overline{Z}}^{2}\right)=\unicode[STIX]{x1D6F1}+C,\end{eqnarray}$$ where the constant of integration 
$C$ is in general still a function of 
$\overline{R}$. In our particular case, we will see that in fact 
$C=0$. The conservation law expressed by (4.23) can be shown to be a special case of conservation laws discovered by Eshelby (Reference Eshelby1975).
4.2 Numerical simulations of the elastic equations
Before we continue with solving the similarity equations, we test the similarity ansatz (4.16) by rescaling full numerical simulations of the elastic equations (4.7)–(4.10). For the numerical treatment, we define the reference state by
$$\begin{eqnarray}R=R_{0}(Z)\unicode[STIX]{x1D702},\quad Z=Z_{0}(\unicode[STIX]{x1D709}),\end{eqnarray}$$ with the elastic domain defined by 
$\unicode[STIX]{x1D702}\in [0,1]$ and 
$\unicode[STIX]{x1D709}\in [0,1]$. The function 
$Z_{0}(\unicode[STIX]{x1D709})$ is computed as part of the numerical solution by imposing the constraint of constant grid spacing in the deformed (current) 
$z$-coordinate along the axis. As a result (cf. figure 7), the function depends on the value of 
$E_{c}$. This accounts for the extreme stretching in the 
$Z$-direction within the thread for small values of 
$E_{c}$, and ensures that enough points remain in the corner region. As can be seen in figure 7, points are concentrated into an increasingly smaller 
$Z$-region, which in turn is stretched out over the entire thread at large deformations.
As an initial condition, we take the free-surface shape
$$\begin{eqnarray}h_{in}(z)=1-\unicode[STIX]{x1D716}\cos (z/2),\end{eqnarray}$$ where 
$\unicode[STIX]{x1D716}=0.005$. We are looking for two unknown functions 
$f$ and 
$g$, where 
$r=r(R,Z)=f(\unicode[STIX]{x1D702},\unicode[STIX]{x1D709})$ and 
$z=z(R,Z)=g(\unicode[STIX]{x1D702},\unicode[STIX]{x1D709})$, as well as the pressure 
$p(\unicode[STIX]{x1D702},\unicode[STIX]{x1D709})$. These three unknowns are found from solving the three equations (4.5) and (4.9). The free surface 
$h(z)$ then is given by the parametric representation 
$h(g(1,\unicode[STIX]{x1D709}))=f_{1}(1,\unicode[STIX]{x1D709})$, from which the curvature 
$\unicode[STIX]{x1D705}$ can be evaluated. Periodic boundary conditions are applied at the ends, as in § 2.3.
The domain is discretised using fourth-order finite differences with 11 equally spaced points in the 
$\unicode[STIX]{x1D702}$-direction and 4001 points in the 
$\unicode[STIX]{x1D709}$-direction. The resulting system of nonlinear equations is solved using a Newton–Raphson technique (Herrada & Montanero Reference Herrada and Montanero2016). We start from the reference state as an initial guess and 
$E_{c}$ sufficiently large (
$E_{c}=100$) to ensure the convergence of the Newton–Raphson iterations. Once we get a solution, we use this solution in a new run with a smaller value of 
$E_{c}$. As seen in figure 6, rescaled stresses as defined by
$$\begin{eqnarray}\overline{\unicode[STIX]{x1D748}}=\frac{\unicode[STIX]{x1D6FE}}{\ell _{e}}\unicode[STIX]{x1D748}\end{eqnarray}$$converge nicely onto a universal similarity solution.
4.3 The thread thickness
Convergence of the elastic solution toward the similarity solution (4.16) for smaller and smaller values of the dimensionless elasto-capillary number 
$E_{c}=GR_{0}/\unicode[STIX]{x1D6FE}$. Stresses have been rescaled according to (4.25). Profiles from the similarity solution as described in § 4.4 are shown as the dashed line for comparison.
The transformation 
$Z_{0}(\unicode[STIX]{x1D709})$ designed to account for stretching in the 
$Z$-direction for three values of 
$E_{c}$.
The interface 
$\overline{h}(\overline{z})$ in similarity variables. The solid line is the rescaled profile obtained from the simulation for 
$E_{c}=10^{-6}$, the dashed line is the solution to the similarity equations. The thread thickness is 
$\overline{h}_{0}$; for large arguments the profile grows exponentially.
In solving the similarity equations, we first concentrate on the free surface profile, shown in figure 8 as a real-space profile in similarity variables. We now show that the dimensionless thread thickness 
$\overline{h}_{0}$ can in fact be calculated without first solving the full similarity equations. Inside the thread, a cylinder of constant dimensionless thickness 
$\overline{R}=1$ (the reference state) is transformed into a thread whose thickness 
$\overline{h}_{0}$ is constant. This elastic problem has as its exact solution the transformation 
$\unicode[STIX]{x1D713}=\overline{h}_{0}\overline{R}$, and thus from incompressibility (cf. (4.17)) 
$\unicode[STIX]{x1D719}_{\overline{Z}}=1/\overline{h}_{0}^{2}$. Then using (4.8), the extensional part of the axial stress becomes 
$\overline{\unicode[STIX]{x1D70E}}_{0}=\unicode[STIX]{x1D719}_{\overline{Z}}^{2}=1/\overline{h}_{0}^{4}$. In particular, the radial stress profile 
$\overline{\unicode[STIX]{x1D70E}}_{0}(\overline{r})$ of § 3 is uniform in the purely elastic case 
$S=\unicode[STIX]{x1D702}_{s}/\unicode[STIX]{x1D702}_{p}=0$. This statement remains true even if the initial condition has a non-trivial shape, as in fact is the case in the elastic simulations to be described below, cf. (4.24). The reason is that, owing to mass conservation, the material inside the thread in the limit of small thread thickness only comes from a tiny region near the point of symmetry. This means that in the normalisation of (4.16) 
$R_{0}$ has to be replaced by 
$1-\unicode[STIX]{x1D716}$, but otherwise the solution stays the same.
Since in the elastic case there is no contribution from the velocity, (3.16) becomes
$$\begin{eqnarray}\overline{T}=2\int _{0}^{\overline{h}}\overline{r}\overline{\unicode[STIX]{x1D70E}}_{zz}\,\text{d}\overline{r}+2\int _{0}^{\overline{h}}\overline{r}(\overline{\unicode[STIX]{x1D705}}-\unicode[STIX]{x1D6F1})\,\text{d}\overline{r}+\overline{h}^{2}\overline{K}.\end{eqnarray}$$Evaluating (4.26) inside the thread of constant thickness, one obtains a finite tension,
$$\begin{eqnarray}\overline{T}=\overline{h}_{0}^{2}/\overline{h}_{0}^{4}+\overline{h}_{0}=\overline{h}_{0}+\frac{1}{\overline{h}_{0}^{2}}>0.\end{eqnarray}$$This corresponds to the earlier conclusion that there must be a positive tension in a liquid bridge (polymeric or Newtonian) in order for pinching to occur (Eggers & Fontelos Reference Eggers and Fontelos2005; Clasen et al. Reference Clasen, Eggers, Fontelos, Li and McKinley2006).
On the other hand, toward the ‘drop’ side of the elastic bridge, where 
$\overline{h}$ becomes large, elastic stresses decay rapidly, as is confirmed by our numerics. As a result, the balance (4.26) is between 
$\overline{T}$ and the second term on the left, which must approach a constant. As a result, it follows that, for large 
$\overline{z}$, the similarity profile must be of the form
$$\begin{eqnarray}\overline{h}=He^{a\overline{z}},\quad a=\frac{2}{\overline{T}}=\frac{2\overline{h}_{0}^{2}}{1+\overline{h}_{0}^{3}}.\end{eqnarray}$$ Using (4.28), both radial and axial contributions to the mean curvature (2.6) scale like 
$\text{e}^{-2a\overline{z}}$, and cancel to leading order, resulting in 
$\overline{\unicode[STIX]{x1D705}}\sim \text{e}^{-4a\overline{z}}$. The prefactor 
$H$ can be normalised to unity by shifting the coordinate system. Thus, putting 
$H=1$ is a way of fixing the origin, ensuring a unique solution.
Applying the conservation law (4.23) to the thread surface 
$\overline{R}=1$, we have 
$\unicode[STIX]{x1D6F1}=\overline{\unicode[STIX]{x1D705}}$, and thus
$$\begin{eqnarray}\left.{\textstyle \frac{1}{2}}\left(\unicode[STIX]{x1D719}_{\overline{Z}}^{2}+\unicode[STIX]{x1D713}_{\overline{Z}}^{2}\right)\right|_{\overline{R}=1}=\overline{\unicode[STIX]{x1D705}}+C.\end{eqnarray}$$ The left-hand side is the same as 
$\text{tr}(\overline{\unicode[STIX]{x1D748}})/2$, which vanishes in the limit 
$\overline{z}\rightarrow \infty$. Since the curvature vanishes as well for 
$\overline{z}\rightarrow \infty$, we have 
$C=0$. On the other hand, for 
$\overline{z}\rightarrow -\infty$, 
$\unicode[STIX]{x1D719}_{\overline{Z}}=1/\overline{h}_{0}^{2}$ and 
$\unicode[STIX]{x1D713}_{\overline{Z}}=0$, while 
$\overline{\unicode[STIX]{x1D705}}=1/\overline{h}_{0}$. As a result, (4.29) yields 
$1/(2\overline{h}_{0}^{4})=1/\overline{h}_{0}$, and in summary
$$\begin{eqnarray}\overline{h}_{0}=2^{-1/3},\quad \overline{\unicode[STIX]{x1D70E}}_{0}=2^{4/3},\quad \overline{T}={\textstyle \frac{3}{2^{1/3}}},\quad a={\textstyle \frac{2^{4/3}}{3}}.\end{eqnarray}$$Thus, we have calculated the dimensionless thread thickness without explicitly solving the similarity equations. Remarkably, these results are identical to those found using lubrication theory (Entov & Yarin Reference Entov and Yarin1984; Clasen et al. Reference Clasen, Eggers, Fontelos, Li and McKinley2006; Eggers & Fontelos Reference Eggers and Fontelos2015), although the similarity solution does not satisfy the slenderness (small slopes) requirement for lubrication to be valid. The reason is that (4.30) only relies on the conservation laws (4.29) and (4.26), which are preserved in the lubrication limit.
4.4 The similarity solution
It remains to calculate the actual form of the profile and of the deformation field inside. We can eliminate 
$\unicode[STIX]{x1D6F1}$ from (4.18), which yields after simplification using (4.17)
$$\begin{eqnarray}\unicode[STIX]{x1D719}_{\overline{R}}\unicode[STIX]{x1D719}_{\overline{Z}\overline{Z}}+\unicode[STIX]{x1D713}_{\overline{R}}\unicode[STIX]{x1D713}_{\overline{Z}\overline{Z}}=\unicode[STIX]{x1D719}_{\overline{Z}}\unicode[STIX]{x1D719}_{\overline{Z}\overline{R}}+\unicode[STIX]{x1D713}_{\overline{Z}}\unicode[STIX]{x1D713}_{\overline{Z}\overline{R}}.\end{eqnarray}$$This can be solved together with the incompressibility constraint (4.17),
$$\begin{eqnarray}\unicode[STIX]{x1D713}\left(\unicode[STIX]{x1D713}_{\overline{R}}\unicode[STIX]{x1D719}_{\overline{Z}}-\unicode[STIX]{x1D713}_{\overline{Z}}\unicode[STIX]{x1D719}_{\overline{R}}\right)=\overline{R}.\end{eqnarray}$$ On the surface 
$\overline{R}=1$, we have (4.22), which in view of (4.29) takes the form
$$\begin{eqnarray}\unicode[STIX]{x1D713}(1,\overline{Z})=\overline{h}(\unicode[STIX]{x1D719}(1,\overline{Z})),\quad \left.{\textstyle \frac{1}{2}}\left(\unicode[STIX]{x1D719}_{\overline{Z}}^{2}+\unicode[STIX]{x1D713}_{\overline{Z}}^{2}\right)\right|_{R=1}=\overline{\unicode[STIX]{x1D705}}.\end{eqnarray}$$ Differentiating the first equation of (4.32), it follows at constant 
$\overline{R}=1$ that 
$\unicode[STIX]{x1D713}_{\overline{Z}}=\unicode[STIX]{x1D719}_{\overline{Z}}\overline{h}_{\overline{z}}$, and hence the second equation can also be rewritten
$$\begin{eqnarray}\frac{1}{2}\unicode[STIX]{x1D719}_{\overline{Z}}^{2}(1,\overline{Z})=\frac{\overline{\unicode[STIX]{x1D705}}}{1+\overline{h}_{\overline{z}}^{2}}.\end{eqnarray}$$ For 
$\overline{Z}\rightarrow -\infty$ the boundary condition is
$$\begin{eqnarray}\unicode[STIX]{x1D713}=2^{-1/3}\overline{R},\quad \unicode[STIX]{x1D719}_{\overline{Z}}=2^{2/3},\end{eqnarray}$$ clearly there can be an arbitrary shift in 
$\overline{Z}$. For 
$\overline{Z}\rightarrow \infty$ we must have
$$\begin{eqnarray}\overline{h}=\unicode[STIX]{x1D713}(\overline{R}=1,\overline{Z})=\text{e}^{a\unicode[STIX]{x1D719}(\overline{R}=1,\overline{Z})},\end{eqnarray}$$ where the prefactor has been chosen 
$H=1$. For the other elastic fields the boundary condition is that the pressure and all the elastic stresses decay for 
$\overline{Z}\rightarrow \infty$.
To obtain the similarity solution numerically, the domain 
$[-L\leqslant \bar{Z}\leqslant L]\times [0\leqslant \bar{R}\leqslant 1]$ was discretised using 
$n_{\bar{R}}$ and 
$n_{\bar{Z}}$ Chebyshev collocation points in the 
$\bar{R}$ and 
$\bar{Z}$ directions, respectively. Here 
$L$ is the value at which the 
$\bar{Z}$ coordinate is truncated. Since strong axial gradients are expected near the origin 
$\bar{Z}=0$ the grid is stretched around that location using the stretching function
$$\begin{eqnarray}\bar{Z}=L\frac{\text{arctanh}(\unicode[STIX]{x1D6FD}s)}{\text{arctanh}(\unicode[STIX]{x1D6FD})}.\end{eqnarray}$$ Here, 
$s$ is in the original Chebyshev interval 
$(-1\leqslant s\leqslant 1)$, and 
$\unicode[STIX]{x1D6FD}$ a stretching parameter.
Equations (4.17) and (4.18) with boundary conditions given by (4.32) and (4.34) are discretised in the numerical domain. At 
$\bar{Z}=L$, the soft boundary condition 
$\unicode[STIX]{x1D719}_{\bar{Z}\bar{Z}}=\unicode[STIX]{x1D713}_{\bar{Z}\bar{Z}}=0$ is used. The resulting discrete system of nonlinear equations have been solved using the MATLAB function fsolve. As an initial guess for the nonlinear solver, the numerical solution for the elastic simulation with 
$E_{c}=1\times 10^{-6}$, rescaled and interpolated into the similarity mesh, was used. Results presented have been obtained using 
$n_{\bar{R}}=5$, 
$n_{\bar{Z}}=450$, 
$\unicode[STIX]{x1D6FD}=0.995$ and 
$L=150$. A larger computational domain by imposing a 
$L>150$ does not alter the results.
The resulting similarity profiles are shown in figures 6–9. In figure 6 we compare representative stress profiles obtained from solutions of the elastic problem with those calculated from the similarity solution; in figure 8 the profile 
$\overline{h}(\overline{z})$ is compared in greater detail. In figure 9(a) we show that the solution to the similarity equations has converged in a large domain, using two different methods. In particular, according to (3.18) and (4.30), 
$\overline{h}_{\overline{z}}/\overline{h}$ converges to 
$a=2^{4/3}/3$, as confirmed in figure 9.
Apart from a full numerical solution of the two-dimensional similarity equations, we also show the result of a different method of solution, explained in more detail in appendix B. In view of the weak radial dependence, we expand the solution into a power series (B 1) in the radial direction. The results of both methods of solution agree very well. In figure 9(b), we show contours of the self-similar pressure 
$\unicode[STIX]{x1D6F1}$ in the 
$[\unicode[STIX]{x1D719},\unicode[STIX]{x1D713}]$ domain. Once more, agreement with contours obtained by rescaling the simulation for 
$E_{c}=10^{-6}$ agree very well with the full similarity solution.
Convergence onto the similarity solution. In (a), the ratio 
$\overline{h}_{\overline{z}}/\overline{h}$ approaches the limit 
$a=2^{4/3}/3$, over a wide range of 
$\overline{Z}$. In (b), we show contours of the function 
$\unicode[STIX]{x1D6F1}$ in the 
$[\unicode[STIX]{x1D719},\unicode[STIX]{x1D713}]$ domain, showing the transition region between the thread and the drop region. Red dashed lines correspond to contours obtained by rescaling the pressure in the elastic simulation with 
$E_{c}=10^{-6}$. Although the radial dependence is weak, the profile is not one-dimensional.
Comparison of the similarity solution as described in § 4.4 with the rescaled profiles (3.6) obtained from the time-dependent simulation shown in figure 1 for 
$De=60$. The time-dependent simulation has been rescaled according to (3.6), with an effective 
$R_{0}$ chosen such that 
$\overline{\unicode[STIX]{x1D70E}}_{zz}=2^{4/3}$ inside the thread.
As a final check, we return to the original problem of the exponential thinning of a viscoelastic liquid bridge, as shown in figure 10. Since the capillary number 
$\overline{v}_{0}=0.035$ is small for our simulations, solvent corrections are expected to be small. Indeed, there is good agreement between our similarity theory based on elastic effects alone, and the full viscoelastic simulations. As we reiterate in the discussion below, there is fading memory in the viscoelastic simulation, and hence there is one adjustable parameter in the comparison, which we choose so that 
$\overline{\unicode[STIX]{x1D70E}}_{zz}=2^{4/3}$ inside the thread; 
$\overline{h}=2/\overline{\unicode[STIX]{x1D70E}}_{zz}$ is then predicted without adjustable parameters. This is an effect of 
$De$ being finite; for the value of 
$De=60$ being used in the present simulation, relaxation effects are almost negligible, as discussed below.
5 Discussion
The aim of this paper has been to find similarity solutions describing the time-dependent pinch-off of a polymeric thread. In the process, we found that this similarity solution also describes another problem, the collapse of a soft elastic bridge under surface tension. This fact had been noticed before (Clasen et al. Reference Clasen, Eggers, Fontelos, Li and McKinley2006), but is generalised here to the full axisymmetric equations. In doing so, we make use of the very general correspondence between fluid flow and nonlinear elasticity (Snoeijer et al. Reference Snoeijer, Pandey, Herrada and Eggers2019).
A particular use of our theory is that we can now calculate the radius of the thread that has formed in the middle of the collapsed bridge. In the purely elastic case, returning to dimensional variables, (4.30) implies that the thread thickness is
$$\begin{eqnarray}h_{thr}=R_{0}\left(\frac{E_{c}}{2}\right)^{1/3},\end{eqnarray}$$ where 
$R_{0}$ is the initial radius of the elastic bridge. In the case of very long relaxation times 
$\unicode[STIX]{x1D706}$, such that an initial thread is formed on short times without the polymer having a chance to relax, followed by exponential relaxation in the limit of long times, (5.1) can be combined with (1.1) to give (Clasen et al. Reference Clasen, Eggers, Fontelos, Li and McKinley2006) 
$h_{0}=R_{0}(E_{c}/2)^{1/3}$, and so
$$\begin{eqnarray}h_{thr}=R_{0}\left(\frac{E_{c}}{2}\right)^{1/3}\exp \left(-\frac{t}{3\unicode[STIX]{x1D706}}\right),\quad De\gg 1.\end{eqnarray}$$ This is illustrated in figure 1(a), where the thread radius is shown for 
$De=60$ and for 
$De=\infty$, such that no relaxation takes place. The initial evolution is almost identical for the two cases, and the thread radius at the beginning of the exponential relaxation (blue line) is just slightly below the asymptotic value without any relaxation at all (red line).
The present work provides an example of a similarity solution in three dimensions, which does not reduce to an effective one-dimensional theory (Chen & Steen Reference Chen and Steen1997; Day, Hinch & Lister Reference Day, Hinch and Lister1998; Cohen et al. Reference Cohen, Brenner, Eggers and Nagel1999; Zhang & Lister Reference Zhang and Lister1999; Eggers & Courrech du Pont Reference Eggers and Courrech du Pont2009; Eggers Reference Eggers2018). It is to be expected that non-trivial similarity solutions in higher dimensions abound, but not many have been found, owing to technical difficulties in obtaining them. Another, and perhaps related feature concerns the universality of similarity solutions. In a one-dimensional approximation, and taking into account effects of the solvent viscosity, one obtains similarity solutions which depend on the value of the parameter 
$\overline{v}_{0}$ alone (Clasen et al. Reference Clasen, Eggers, Fontelos, Li and McKinley2006).
Here, we provide evidence (cf. figure 5) that the three-dimensional, axisymmetric solution to the Oldroyd-B model is far less universal, and that there exists a similarity solution for a given stress distribution 
$\overline{\unicode[STIX]{x1D70E}}_{0}(\overline{r})$ inside the thread. To answer this question conclusively, one would have to solve the full similarity equations (3.7)–(3.8) for a prescribed profile 
$\overline{\unicode[STIX]{x1D70E}}_{0}(\overline{r})$, which we have not yet attempted. In the one-dimensional description of Clasen et al. (Reference Clasen, Eggers, Fontelos, Li and McKinley2006), the equivalent of the dimensionless coefficient 
$\overline{h}_{0}$ (cf. (4.30)) was seen already to depend on 
$\overline{v}_{0}$. In the present three-dimensional, axisymmetric case and for finite solvent viscosity 
$S>0$, we expect it to depend on 
$\overline{\unicode[STIX]{x1D70E}}_{0}(\overline{r})$ as well. In the absence of a solvent (
$S=0$) one recovers a uniform stress profile over the cross-section of the thread, and, strictly speaking, (5.2) only applies to that case. It is however worth while keeping in mind that these corrections coming from the solvent viscosity are small in practice. For example, for the experiments by Clasen et al. (Reference Clasen, Eggers, Fontelos, Li and McKinley2006), using a quite elevated solvent viscosity, 
$\overline{v}_{0}=0.031$; for the experiments by Sattler et al. (Reference Sattler, Wagner and Eggers2008), using water as a solvent, 
$\overline{v}_{0}=1.1\times 10^{-3}$.
Finally, it is worth commenting on the experimental situation, since the Oldroyd-B model considered here may be popular for its conceptual simplicity, but cannot be considered a first-principles description of polymeric flow. While an exponential thinning regime of the polymer thread is observed almost universally, indicating exponential stretching of polymer strands, there are experimental indications of a more abrupt change in polymer conformation (Ingremeaux & Kellay Reference Ingremeaux and Kellay2013).
As for the shape of the self-similar profile, Turkoz et al. (Reference Turkoz, Lopez-Herrera, Eggers, Arnold and Deike2018) find that the experimental free-surface profile grows significantly faster toward the drop when compared to full numerical simulations of the three-dimensional Oldroyd-B equations. This means that the theoretical result 
$a=2^{4/3}/3$ for the growth exponent in the free-surface profile (4.28) underpredicts the experimentally observed value by approximately a factor of two. A possible explanation could be that one has to take the finite extensibility of polymers into account, as it is done in the FENE-P model (Bird et al. Reference Bird, Armstrong and Hassager1987; Morozov & Spagnolie Reference Morozov, Spagnolie and Spagnolie2015).
Now that a fully quantitative prediction of the Oldroyd-B fluid is finally in place, and given the time that has passed since the original experiments (Clasen et al. Reference Clasen, Eggers, Fontelos, Li and McKinley2006), it would be well worth repeating the measurements for a different system. Another avenue to pursue is to allow for polymer models with more free parameters, such as a spectrum of relaxation times. It might well be that different regions between the drop and the thread are governed by different time scales.
Acknowledgements
We gratefully acknowledge illuminating discussions with M. Fontelos during the early stages of this research. J.E. acknowledges the support of Leverhulme Trust International Academic Fellowship IAF-2017-010, and is grateful to H. Stone and his group for welcoming him during the academic year 2017–2018. M.A.H. thanks the Ministerio de Economía y competitividad for partial support under the project no. DPI2016-78887-C3-1-R. J.H.S. acknowledges support from NWO through VICI grant no. 680-47-632.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Force balance in Lagrangian coordinates
In order to solve the similarity equations by an expansion in the radial variable, described in appendix B below, it is useful to implement the force balance (4.26) in Lagrangian coordinates. To this end, within the leading-order balance (4.19), we use the 
$z$-component of (3.12)
$$\begin{eqnarray}\int _{O}\overline{\unicode[STIX]{x1D705}}\boldsymbol{n}\,\text{d}s=-\left.\frac{2\unicode[STIX]{x03C0}\overline{h}}{\sqrt{1+\overline{h}_{\overline{z}}^{2}}}\right|_{\overline{Z}_{-}}^{\overline{Z}_{+}}=\int _{C_{\pm }}\boldsymbol{n}_{\pm }\boldsymbol{\cdot }\overline{\unicode[STIX]{x1D748}}\boldsymbol{\cdot }\boldsymbol{e}_{z}\,\text{d}s.\end{eqnarray}$$ However, now 
$O$ is the surface of the thread between two fixed Lagrangian coordinates 
$\overline{Z}_{\pm }$, and 
$C_{\pm }$ are the surfaces
$$\begin{eqnarray}\overline{z}=\unicode[STIX]{x1D719}(\overline{R},\overline{Z}_{\pm }),\quad \overline{r}=\unicode[STIX]{x1D713}(\overline{R},\overline{Z}_{\pm }),\end{eqnarray}$$ parameterised by 
$\overline{R}$ between 0 and 1; here, 
$\boldsymbol{n}_{\pm }$ are the outward normals to these surfaces. Inside the thread, 
$\overline{z}$ is constant, and hence 
$\boldsymbol{n}_{-}=-\boldsymbol{e}_{z}$, while on the right
$$\begin{eqnarray}\boldsymbol{n}_{+}=(\unicode[STIX]{x1D713}_{\overline{R}},-\unicode[STIX]{x1D719}_{\overline{R}})/\sqrt{\unicode[STIX]{x1D719}_{\overline{R}}^{2}+\unicode[STIX]{x1D713}_{\overline{R}}^{2}},\end{eqnarray}$$which follows from the parameterisation (A 2).
The integral over the extra part of the stress is
$$\begin{eqnarray}\displaystyle \int _{C_{+}}\boldsymbol{n}_{+}\boldsymbol{\cdot }\overline{\unicode[STIX]{x1D748}}_{p}\boldsymbol{\cdot }\boldsymbol{e}_{z}\,\text{d}s & = & \displaystyle 2\unicode[STIX]{x03C0}\int _{0}^{1}\boldsymbol{n}_{+}\boldsymbol{\cdot }\overline{\unicode[STIX]{x1D748}}_{p}\boldsymbol{\cdot }\boldsymbol{e}_{z}\unicode[STIX]{x1D713}\sqrt{\unicode[STIX]{x1D719}_{\overline{R}}^{2}+\unicode[STIX]{x1D713}_{\overline{R}}^{2}}\,\text{d}\overline{R}\nonumber\\ \displaystyle & = & \displaystyle 2\unicode[STIX]{x03C0}\int _{0}^{1}\left(\unicode[STIX]{x1D713}_{\overline{R}}\unicode[STIX]{x1D719}_{\overline{Z}}^{2}-\unicode[STIX]{x1D719}_{\overline{R}}\unicode[STIX]{x1D713}_{\overline{Z}}\unicode[STIX]{x1D719}_{\overline{Z}}\right)\unicode[STIX]{x1D713}\,\text{d}\overline{R}=2\unicode[STIX]{x03C0}\int _{0}^{1}\unicode[STIX]{x1D719}_{\overline{Z}}\overline{R}\,\text{d}\overline{R},\nonumber\end{eqnarray}$$having used incompressibility (4.17) in the last step. The pressure contribution is
$$\begin{eqnarray}-\int _{C_{+}}\unicode[STIX]{x1D6F1}\boldsymbol{n}_{+}\boldsymbol{\cdot }\boldsymbol{e}_{z}\,\text{d}s=-2\unicode[STIX]{x03C0}\int _{0}^{1}\unicode[STIX]{x1D6F1}\unicode[STIX]{x1D713}_{\overline{R}}\unicode[STIX]{x1D713}\,\text{d}\overline{R}.\end{eqnarray}$$Inside the thread,
$$\begin{eqnarray}\boldsymbol{n}_{-}\boldsymbol{\cdot }\overline{\unicode[STIX]{x1D748}}\boldsymbol{\cdot }\boldsymbol{e}_{z}=-\unicode[STIX]{x1D748}_{zz}={\textstyle \frac{1}{2}}\left(\unicode[STIX]{x1D713}_{\overline{Z}}^{2}+\unicode[STIX]{x1D719}_{\overline{Z}}^{2}\right)=-2^{1/3},\end{eqnarray}$$ having used (4.34) and (4.30). Thus the integral over the cross-section 
$\unicode[STIX]{x03C0}\overline{h}_{0}^{2}$ yields
$$\begin{eqnarray}\int _{C_{-}}\boldsymbol{n}_{\pm }\boldsymbol{\cdot }\overline{\unicode[STIX]{x1D748}}\boldsymbol{\cdot }\boldsymbol{e}_{z}\,\text{d}s=-\unicode[STIX]{x03C0}2^{-1/3}.\end{eqnarray}$$ Finally, since inside the thread 
$\overline{h}=\overline{h}_{0}$, (A 1) yields
$$\begin{eqnarray}\frac{2\overline{h}}{\sqrt{1+\overline{h}_{\overline{z}}^{2}}}+2\int _{0}^{1}\unicode[STIX]{x1D719}_{\overline{Z}}\overline{R}\,\text{d}\overline{R}-\int _{0}^{1}(\unicode[STIX]{x1D719}_{\overline{Z}}^{2}+\unicode[STIX]{x1D713}_{\overline{Z}}^{2})\unicode[STIX]{x1D713}_{\overline{R}}\unicode[STIX]{x1D713}\,\text{d}\overline{R}=\overline{T}=\frac{3}{2^{1/3}},\end{eqnarray}$$ where we inserted the expression (4.23) for the pressure, with 
$C=0$.
Appendix B. Solution by radial expansion
The similarity profiles only have a weak dependence in the 
$\overline{R}$-direction, and are thus well represented by a polynomial in 
$\overline{R}$. Thus we seek a solution by expanding in the 
$\overline{R}$-coordinate according to
$$\begin{eqnarray}\unicode[STIX]{x1D719}=\unicode[STIX]{x1D719}_{0}+\unicode[STIX]{x1D719}_{2}\overline{R}^{2}+\cdots \,,\quad \unicode[STIX]{x1D713}=\unicode[STIX]{x1D713}_{1}\overline{R}+\unicode[STIX]{x1D713}_{3}\overline{R}^{3}+\cdots \,.\end{eqnarray}$$First, from incompressibility we have
$$\begin{eqnarray}\unicode[STIX]{x1D713}_{1}=\frac{1}{\unicode[STIX]{x1D719}_{0}^{\prime 1/2}},\quad \unicode[STIX]{x1D713}_{3}=-\frac{\unicode[STIX]{x1D719}_{0}^{\prime }\unicode[STIX]{x1D719}_{2}^{\prime }+\unicode[STIX]{x1D719}_{0}^{\prime \prime }\unicode[STIX]{x1D719}_{2}}{4\unicode[STIX]{x1D719}_{0}^{\prime 5/2}}.\end{eqnarray}$$ Expanding (4.31) in a power series in 
$\overline{R}$ we find at first order that
$$\begin{eqnarray}-4\unicode[STIX]{x1D719}_{0}^{\prime 4}\unicode[STIX]{x1D719}_{2}^{\prime }+4\unicode[STIX]{x1D719}_{0}^{\prime 3}\unicode[STIX]{x1D719}_{0}^{\prime \prime }\unicode[STIX]{x1D719}_{2}-\unicode[STIX]{x1D719}_{0}^{\prime }\unicode[STIX]{x1D719}_{0}^{\prime \prime \prime }+\unicode[STIX]{x1D719}_{0}^{^{\prime \prime }2}=0,\end{eqnarray}$$which is valid exactly, as an equation between Taylor coefficients.
The second equation is (4.33), which can be satisfied approximately by truncating the expansion (B 1) after the first two terms. Hence, we obtain
$$\begin{eqnarray}\frac{\left(\unicode[STIX]{x1D719}_{0}^{\prime }+\unicode[STIX]{x1D719}_{2}^{\prime }\right)^{2}}{2}=\frac{\overline{\unicode[STIX]{x1D705}}}{1+\overline{h}_{\overline{z}}^{2}},\end{eqnarray}$$where
$$\begin{eqnarray}\overline{h}=\unicode[STIX]{x1D713}_{1}+\unicode[STIX]{x1D713}_{3},\quad \overline{h}_{\overline{z}}=\frac{\unicode[STIX]{x1D713}_{1}^{\prime }+\unicode[STIX]{x1D713}_{3}^{\prime }}{\unicode[STIX]{x1D719}_{0}^{\prime }+\unicode[STIX]{x1D719}_{2}^{\prime }}\quad \overline{h}_{\overline{z}\overline{z}}=\frac{\unicode[STIX]{x1D713}_{1}^{\prime \prime }+\unicode[STIX]{x1D713}_{3}^{\prime \prime }}{\left(\unicode[STIX]{x1D719}_{0}^{\prime }+\unicode[STIX]{x1D719}_{2}^{\prime }\right)^{2}}-\frac{\left(\unicode[STIX]{x1D713}_{1}^{\prime }+\unicode[STIX]{x1D713}_{3}^{\prime }\right)\left(\unicode[STIX]{x1D719}_{0}^{\prime \prime }+\unicode[STIX]{x1D719}_{2}^{\prime \prime }\right)}{\left(\unicode[STIX]{x1D719}_{0}^{\prime }+\unicode[STIX]{x1D719}_{2}^{\prime }\right)^{3}},\end{eqnarray}$$ and 
$\unicode[STIX]{x1D713}_{1},\unicode[STIX]{x1D713}_{3}$ is expressible through 
$\unicode[STIX]{x1D719}_{0}$ and 
$\unicode[STIX]{x1D719}_{2}$ using (B 2). We have to solve the system of (B 3), (B 4) for the two variables 
$\unicode[STIX]{x1D719}_{0}$ and 
$\unicode[STIX]{x1D719}_{2}$ with conditions 
$\unicode[STIX]{x1D719}_{0}^{\prime }=2^{2/3}$ and 
$\unicode[STIX]{x1D719}_{2}=0$ for 
$\overline{Z}\rightarrow -\infty$.
A problem with (B 4) is that it is of quite high order. A much better alternative is to use (A 3), because it only requires first derivatives 
$\overline{h}_{\overline{z}}$, and because it implements the force balance exactly; inserting (B 1), the integrals can be done easily. The approximated solution was obtained with the same numerical procedure used to obtain the solution to the full similarity equation.
