Wetting and dewetting processes in the axial retraction of liquid filaments

We study the hydrodynamic mechanisms involved in the motion of the contact line formed at the end region of a liquid filament laying on a planar and horizontal substrate. Since the flow develops under partially wetting conditions, the tip of the filament recedes and forms a bulged region (head), that subsequently develops a neck region behind it. Later on, the neck breaks up leading to a separated drop, while the rest of the filament restarts the sequence. One main feature of this flow is that the whole dynamics and final drop shapes are strongly influenced by the hysteresis of the contact angle typical in most of the liquid/substrate systems. The time evolution till breakup is studied experimentally and pictured in terms of a hybrid wettability theory which involves the Cox-Voinov hydrodynamic approach combined with the molecular kinetic theory developed by Blake. The parameters of this theory are determined for our liquid/substrate system (silicon oil / coated glass). The experimental results of the retracting filament are described in terms of a simple heuristic model, and also compared with numerical simulations of the full Navier-Stokes equations. This study is of special interest in the context of pulsed laser induce dewetting (PLiD).


I. INTRODUCTION
Controlling the size, shape and placement of nanoparticles is a main concern in many processes of nanostructures fabrication. For instance, the incorporation of plasmonic nanoparticles into photovoltaic devices can strongly increase their efficiency [1,2], since the surface plasmon resonance between metallic nanostructures depends on the size and spacing [3,4] of the particles. Other applications include biodiagnostics and sensing, where Au nanoparticles bind to specific DNA markers, permitting their detection [5]. The potential applications for organized metallic nanostructures are wide ranging and include Raman spectroscopy [6,7], catalysis [8], photonics [9], and spintronics [10]. One strategy to create and organize structures at the nanoscale is to make use of the inherent self-assembly mechanisms of a material. The combination of natural instabilities related to the physical properties of liquid metals, such as low viscosity and high surface energy, with the technological ability to lithographically pattern nanoscale features, has created an open field of research. These techniques require the study of the liquid-state dewetting dynamics [11] leading to liquid instabilities [12], with the goal of directing the assembly to precise, coordinated nanostructures in one [13,14] and two [15,16] dimensions. Within this context, the phenomena of wetting and dewetting of a shaped liquid film (metallic or not) covering a solid substrate are at the basis of all these processes and applications. However, even if these phenomena have been studied for the last two decades [17,18], several of the main issues related to them are still under discussion. Therefore, one of the purposes of this paper is to provide a clear description of the wetting and dewetting mechanisms, including effects of contact angle hysteresis, and their implications in the formation of liquid structures.
We consider in detail the dynamics of a long liquid filament sitting on a horizontal substrate under partial wetting conditions. In particular, we focus our attention on the bulged region at its ends that is formed due to the filament retraction along its longitudinal direction. This process not only involves a dewetting motion, but also a wetting one in the transverse direction. Due to these advancing and receding motions at the ends, the contact line becomes curved in the region connecting the bulge (head) with the rest of the filament. This perturbation triggers an instability that leads to the appearance of a neck with rapidly decreasing width [19,20]. The breakup of this neck leads to two fluid portions which evolve separately. On one side, the head separated from the filament recedes until a final drop shape is achieved. The main features of this drop, such as its non-circular footprint, was thoroughly studied in previous works [19,21]. On the other side, the remaining filament has a new shaped end as a result of the breakup process. The following receding longitudinal motion of this new end of the filament leads to another bulged region and a neck formation in a successive fashion. A group of pictures of this process, as obtained in this work for a silicon oil on a coated glass, is shown in Fig. 1(a). For comparison, we also include in Fig. 1(b) the analogous process observed at nanoscale [13,16] for a filament of width 293 nm and thickness 141 nm, which was generated in the context of pulsed laser induce dewetting (PLiD) (see e.g. [16] and references cited there).
The main goal of this work is to describe and understand the physical mechanisms involved in the evolution of the head and its retraction. The flow before and after the breakup is studied experimentally by means of optical techniques, which allow to measure contact line positions, contact angles, thickness profiles and contact line contours. These data are then compared with a simple heuristic model as well as numerical simulations, which solve the Navier-Stokes equations with appropriate boundary conditions at the contact line. Note that the modeling of this flow implies Similar process observed at nanoscale for a melted Ni filament [16].
the description of the contact line movement, and the subsequent pressure and velocity fields in a fluid domain confined by the substrate and the fluid free surface. Thus, the main issue in this modeling is how to determine the contact line motion, where both the stress and viscous dissipation rate diverge in the continuum description if slip is not allowed.
Here, we consider a small amount of slip within the Navier boundary condition, and consistently add a relationship between the velocity and the contact angle. Since these relations must account for wetting and dewetting, we choose the combined approach proposed in [22], which accounts for both the hydrodynamic of the macroscopic region [23] and the molecular kinetic [24,25] for the contact line motion. The relevant physical parameters in this approach, which leads to a constitutive relationship for the contact line dynamics, are determined here by an additional study of the hysteretic behavior of a sessile drop.
In Section II, we obtain the main wettability properties of the substrate and liquid (silicon oil) system to be used in the filament experiments described in Section III. The experimental data are compared both with a simple heuristic model (Section IV) and numerical simulations of the Navier-Stokes equations (Section V). Section VI is devoted to a summary and main conclusions of the work.

II. WETTABILITY OF THE SYSTEM SUBSTRATE/LIQUID
In order to perform the experiments, we employ a similar setup to that reported in [21], i.e. we use a substrate that is partially wetted by our working fluid, namely a silicon oil (polydimethylsiloxane, PDMS). The substrate is a microscope slide (glass) which is coated with a fluorinated solution (EGC-1700 of 3M) using a Chemat Dip Coater under controlled speed. The PDMS partially wets the substrate, since the coating lowers the surface energy of the glass. In order to have reproducible wetting properties, we found that it is necessary to follow a protocol to get rid of the remaining solvent in the coating. Thus, the coated substrates are heated in an oven at 30 • C for 30 minutes, and left on a close and dry vessel for about 3 days.
In order to characterize the wettability of the PDMS, we employ a technique which advances or recedes the contact line by changing the volume of a sessile drop initially placed on the substrate [26]. This is accomplished by injecting or withdrawing liquid through of a needle (connected to a syringe) in contact with the apex of the drop (see inset in Fig. 2(a)). The syringe volume is controlled by an Automated Dispensing System of a Rame-Hart Model 250 goniometer, and the measurements are based on the analysis of the axisymmetric drop shape profile (ADSA-P technique).
We distinguish here two scenarios: one in which the injection/withdrawal process is produced by pulses (between which the drop is allowed to reach a static shape), and another one where this process is practically continuous. Thus, in the first case we focus on the measurement of the static contact angle, θ e , while in the second case we are concerned with the dynamic angle, θ, and its relationship with the drop contact line velocity.

A. Static wettability
We start with a drop of volume V 0 = 25µl and use injection/withdrawal pulses of 1 s duration at a flow rate Q = ±1µl/s. This is done by means of an automated control of the piston motion in the syringe. The pulses are separated by time intervals of about 15 s, so that the drop has time enough to relax to its static shape after every volume variation, ∆V = 1µl.
Initially at t = 0, we have a sessile drop of volume V 0 with θ e = 56 • and null contact line displacement (δx = 0). As the volume is increased, it spreads and reaches intermediate rest states on a dry surface (not previously covered by the PDMS) till the maximum volume, V max = 40µl, is achieved. This is shown by the blue symbols in Fig. 2(a) for θ e (full circles) and δx (hollow circles), respectively. Note that θ e remains practically the same for all these volumes, while δx increases. This stage is not part of the hysteresis cycle, since only contact line motions on a previously wetted surface are considered in its determination.
After this preliminary stage, the needle withdraws liquid until V diminishes to V min = 10 µl (see the portion of the cycle from A to B in Fig. 2(a): clockwise for θ, and counterclockwise for δx, as indicated by the arrows). Finally, the increasing volume stage (from B to A) starts at V min till V max is reached again. This process can be repeated indefinitely, but only the first hysteresis cycle is shown in Fig. 2(a) for brevity (the full black circles stand for θ e , while the hollow red ones correspond to δx). In order to fully understand the meaning of these curves in terms of the local behavior at the contact line, we plot θ e versus δx in Fig. 2(b), and ignore V since it is not a relevant parameter for this analysis. Here, we see that the front starts to have a forward displacement (advance) at θ a = 52 • to achieve a new static position, so that this angle is called the (static) advancing contact angle. Similarly, we define the (static) receding contact angle, θ r = 46 • , as the value of θ e at which the front must have a backwards displacement (recede) to achieve equilibrium. As it will be shown in the following sections, these two limiting angles are related to the dynamic contact angle behavior. Another relevant couple of angles is given by the maximum and minimum values of θ e , between which static drop states are possible. This range is (θ min , θ max ) = (40 • , 55 • ).

B. Dynamic wettability
Now, we consider the relationship between the dynamic contact angle, θ, and the contact line velocity, v cl . Here, we inject/withdraw the liquid at the same flow rate, Q, as before, but with longer and larger pulses. Except for the first pulse (at t = 0), which has a duration of ∆t = 15 s with a volume increase of ∆V = 15 µl, all the other ones last ∆t = 30 s with a volume variation of ∆V = ±30 µl (see dashed lines in Fig. 3). Unlike the static case, we continuously measure the dynamic contact angle, θ(t), and the corresponding contact line displacement,∆x(t), as the drop spreads/contracts and relaxes to equilibrium for about 200 s. The values of θ and ∆x as a function of time are plotted in Fig. 3. Note that during the injection and withdrawal stages θ varies in phase (without delay) with V , while the effect of ∆V on the contact line displacement, ∆x, is delayed. This is due to the fact that during these stages the drop volume variation is quickly absorbed by a modification of the contact angle, while the contact line does not move till θ reaches the critical contact angle (θ a for advancing, and θ r for receding). We derive numerically the data of ∆x versus t to obtain the contact line velocity, v cl (t). The hollow circles in Fig. 4 show θ as a function of v cl , using t as a parameter. Interestingly, the θ range for which v cl = 0 is practically coincident with the hysteresis range of θ e , namely (θ r , θ a ) = (46 • , 52 • ) shown in Fig. 2(b). This agreement confirms that the methodology used to determine both types of angles, namely static and dynamic, is appropriate to describe the wettability of the substrate by the PDMS.
In order to fully describe the dewetting stage, we must enlarge the velocity range for v cl < 0. In fact, it is known that the dewetting velocity has a maximum absolute value at which the dynamic contact angle is zero [22,27], and our results in Fig. 4 for the dynamic cycle of the sessile drop cannot show this phenomenon.
With this goal, we perform a different experiment which allows to obtain the maximum velocity of dewetting, referred here to as v max (> 0). In this new setup, we plunge a coated substrate into a deep pool of PDMS, and withdraw it at a speed, V s , greater than v max . This is achieved when, for a certain V s , the liquid is entrained off the pool, and forms a film on the substrate. Once the substrate has stopped, the contact line of this film starts to slide down the surface due to gravity. This dewetting motion proceeds at the maximum possible velocity, v max , with zero contact angle [28]. We perform experiments with different withdrawal speeds (0.033 < V s < 0.5 cm/s), and confirmed that the value of v max is independent of V s , as expected. The contact lines on both sides of the coated glass are seen on the digital image due to the transparency of the substrate. They are not exactly coincident, thus introducing some error in the measurement, which is at most 7% and increases as the front dewets. Finally, by measuring the contact line position as a function of time, we are able to obtain the linear regression expression v max = (5.73 ± 0.24) × 10 −3 cm/s.
This value of v max is shown in Fig. 4 for θ = 0, and it is of great importance to determine the parameters of the physical mechanism that describes the θ-v cl relationship, as will be shown in the following section.

C. Wettability model
In order to describe the motion of the contact line over the substrate (xy-plane), we employ two constitutive relations to account for the liquid-solid interaction. Firstly, we overcome the stress singularity by relaxing the no slip boundary condition at the substrate through the Navier formulation (see e.g. [29]) , where is the slip length. This parameter embraces somehow the main features related with the intermolecular forces in that region. Secondly, we use a relationship between the dynamic contact angle, θ, and the contact line velocity, v cl , which is normal to the contact line itself. Since we need to describe not only advancing motions (wetting, v cl > 0), but also receding (dewetting, v cl < 0) ones, we define the dynamic contact angle from a combination of the hydrodynamic and molecular kinetics models [24,25,30] as first proposed in [22]. For the former, we have the Cox-Voinov relationship [23,31] where θ m is the microscopic contact angle, and L, are macroscopic and microscopic length scales, respectively. Usually, the former is assumed to be a constant given by Young's law, and independent of v cl . Here, we take L as the capillary length, i.e.
where g is the gravity. Thus, we leave as a parameter to be determined by adjusting the experimental data from the final modeling. Note, however, that Eq. (3) with θ m = const. is not appropriate to adjust well the experimental data in Fig. 4, which show a sharp jump at v cl = 0 due to hysteresis effects. Then, we must resort to an alternative approach. In fact, several authors [23,31] admit the possibility of a non-hydrodynamic velocity dependence of θ m . Here, we consider that it is given by the Blake molecular-kinetics model [25] in the form, where θ 0 is the (microscopic) equilibrium contact angle (v cl (θ 0 ) = 0). The other parameters, which are of molecular origin, are given by [25] v 0 = 2κλ, where κ is the frequency of molecular displacement at equilibrium, λ is the average length of each molecular displacement (or distance between adsortion sites), T is the temperature, and k is the Boltzmann's constant. Since we do not focus our attention in the kinetic process itself, we consider v 0 and Γ as fitting parameters of our modeling. On the other hand, since the surface energies are usually associated with the cosine of a contact angle, as in Young's law, we consider the cosine of the equilibrium angle, θ 0 , as an average of the cosines of two characteristic angles of the hysteresis range. If we take the static advancing and static receding contact angles, θ a and θ r , respectively, we have which yields θ 0 = 49.07 • . We could also have taken the pair (θ min , θ max ) instead of (θ r , θ a ), in which case the angle would be θ * 0 = 47.95 • . However, as we will discuss later, this choice does not lead to any improvement in the comparison with the experimental results.
Note that both Eqs. (3) and (5) account for the hysteresis of the contact line, since they provide the dynamic contact angle, θ, as a function of v cl . In order to determine the parameters v 0 , Γ, and , we fit the data from the dynamic cycling contact angle measurements (circles in Fig. 4), and the value (v max , 0) (rhombic symbol in Fig. 4, and Eq. (1)) with the approximating function This combined Cox-Voinov-Blake model (for brevity CV+B, or hybrid model) is solved with the following iterative procedure. For a given pair ( (1) ) as given by a square minimum method using Eq. (8) with the dynamic cycling measurements. By using this pair, we determine a new value of v The best fitting curve given by Eq. (8) with these parameters is shown in Fig. 4 (see blue line). The red line is the corresponding curve using θ * 0 instead of θ 0 . Even if the difference between them is small and is mainly noticeable in the dewetting region, we will use in what follows the parameters in Eq. (9) for the blue curve in Fig. 4. These parameters allow us to estimate the molecular kinetic parameters for our system coated glass/PDMS/air as given by Eq. (6). Thus, we find Interestingly, λ is of the same order as those reported in [22], but κ is much smaller than for those systems with glycerol, likely due to the higher viscosity of PDMS. At the same time, is at least ten times larger here than in [22].

III. DYNAMICS OF A LIQUID FILAMENT ON THE SUBSTRATE: EXPERIMENTS
The main dynamical process that we analyze here is the axial dewetting motion of a liquid filament placed on a horizontal and previously coated glass substrate [19,21,32]. We generate the filament from a vertical jet of PDMS flowing out from a small nozzle at the bottom of a vessel filled with PDMS. It is captured on the substrate by horizontally flipping the substrate in a rigid frame. Quickly afterwards, the frame is rotated 90 degrees and a second horizontal flip is performed so that a new (auxiliary) filament, perpendicular to the first one, is also captured on the substrate. This other filament crosses the first one near one of its ends, so that the latter adopts a reproducible and characteristic rounded shape, whose axial dewetting will be observed and measured in detail. Finally, the substrate is rapidly placed on a horizontal position. All these movements take about 5 seconds, which is a very short time interval compared to the time scale of the experiment itself due to the high viscosity of the PDMS.
Thus, we obtain a fluid filament of uniform width, w, with parallel and straight contact lines, so that the initial configuration has a constant cross sectional area along its axis. We calibrate the system by relating the fluid height in the vessel and the jet diameter for a given nozzle. Both the jet diameter and the corresponding filament widths (for a contact angle equal to θ a ) could be varied from 0.3 to 1.3 mm and from 0.1 to 1.0 mm, respectively.

A. Axial dewetting between breakups
A typical dewetting process of the filament end is shown in Fig. 1. Initially, it is rounded and has a contact angle of around 25 • so that the tip certainly recedes according to the wettability framework described above (see e.g. Fig. 4). This motion induces the development of a head, whose width and thickness grow with time. The evolution proceeds until the filament tip stops and, after a while, a narrower region (neck) starts forming in the filament somewhere beyond the bulged region. Finally, the neck breaks up, so that the portion of fluid between the tip and the neck gives place to a static drop. The characteristic shape of this sessile drop with non-circular footprint was previously studied in [21], but here we will focus on the description of the dynamical process previous to the rupture. This dewetting and breakup mechanism repeats itself starting from a new end of the filament, and finishes with the formation of a series of similar drops. Here, we show results for the tracking of the first four heads at the left end of this filament.
The experimental setup allows to observe, though not simultaneously, both top and side views of the filament evolution. In Fig. 5, we show the filament profiles for two similar, albeit not strictly equal, widths. From this type of profiles, we are able to extract the position of the tip, the contact angle there, as well as the thickness and width of the head. These parameters are later compared with the results from an heuristic model and numerical simulations of the full hydrodynamic equations.

B. Evolution of fronts, contact angles, and thicknesses
It should be noted that the first drop formed from the end of the filament is a consequence of a different breakup process (numbered here as head 0) than the following ones. This is due to the fact that its preceding breakup occurs at crossings with a transversal auxiliary filament, unlike the rest of the drops. Since the initial conditions of the formation process determines important features of the dynamics, such as the distance travelled by the retracing tip till the following breakup, we shall restrict ourselves to the most repetitive case. Therefore, we exclude it in our analysis.
In Fig. 6 we show with symbols the time evolution of the axial position of the tip (front), x f , the dynamic contact angle there, θ x , and the maximum thickness of the head for the first four heads (1 to 4) formed at the left end of the filament of width w a (see Fig. 5(a)). The origin of x f -coordinates for successive heads is the point where each breakup took place. The almost perfect superposition of the data for all three quantities (x f , θ x , and h 0 ) indicates that the dewetting process is repeatable after each breakup.
From the data presented in Fig. 6, we generate the corresponding θ x -v cl relation, as shown by symbols in Fig. 7. Interestingly, these spontaneous filament dewetting data are in good agreement with the previously obtained relationship θ(v cl ) for the drop under the cyclic procedure of forced wetting and dewetting motions (solid line in Fig. 7). This fact confirms that the previously obtained relationship certainly describes a local phenomenon, and is independent of the geometry of the flow.
The measured evolution of the footprint of a filament can be seen in Fig. 8 for two different filament widths (see also Fig. 5(b)). The results show that all footprints are inside an angular sector whose borders are always tangent at  some point (x 1 , y 1 ) of the contact line. The aperture of this sector is fixed, so that the slopes of the straight lines are constant during the axial retraction process, except when the neck behind the head becomes very narrow. The straight lines satisfy the expression y = ±(0.114 x + 0.031) cm, which means that they intersect at (x * , 0) = (−0.2719, 0) cm. We take this point as the origin of coordinates for the following discussions, so that this envelope can be described as straight lines with slopes α = 6.5 • . The physical meaning of this result will become apparent from the analysis of the full simulations in Section V. In Fig. 9 we show the tip position, x f , the maximum width at the head, w head , and the minimum width at the neck, w neck , as a function of time for the filament width w b = 0.0665 cm. As expected, the curve for x f is very similar to that in Fig. 6(a). In particular, that for early times, say t < 100 s, the data for both widths are practically coincident, the main difference being for later times when x f approaches the corresponding limiting value. On the other hand, Fig. 9(b) shows w head and w neck , which are obtained as the maximum and minimum values of the width along the  filament. Note that the growth rate of the head slightly diminishes when the neck starts developing (t ≈ 40 s).

IV. DESCRIPTION OF THE AXIAL DEWETTING: HEURISTIC MODEL
In order to develop a simple model that accounts for the main features of the dynamic process of dewetting, we consider that the shape of the head at the end of the filament can be approximated as a part of an ellipsoidal cap (see Fig. 10). This simple geometrical configuration continuously connects to the rest of the filament of circular cross section. The semi-diameters along x, y and z axes are a, b and c, respectively, which we consider as constants during the whole evolution. The distance between the center of the ellipsoid and the substrate in the vertical direction is z 0 , and the cap has length 2x 0 along x-axis and width 2y 0 in the y direction, as shown in Fig. 10(b). The key idea is that the head can be emulated with the part of this ellipsoid that is above the substrate level. This cap volume, V , varies as the dewetting motion of the tip proceeds. The dynamics can be described as the superposition of both a vertical and a horizontal motion of the ellipsoid, with velocities dz 0 /dt and U , respectively. Since the ellipsoid is not only emerging but also moving horizontally, we describe it in terms of the coordinates of its center, (x c , 0, z c ), and thus, it is described by the equation Then, the volume of the ellipsoidal cap is given by Due to the horizontal retraction of the tip, there is a mass accretion in the cap from the filament that is being swept by the receding head. Therefore, the cap volume V (z 0 ) increase can be calculated in first approximation by where A 0 , the area of the cross section of the filament of width w and transverse contact angle θ a , is given by Here, we define x 0 as the distance between the tip and the center position, x c , i.e. x 0 = x f − x c . Note that the vertical motion of the ellipsoid yields a variation of x 0 . Therefore, the contact line velocity along x-axis is the result of the superposition of the bulk horizontal displacement and the variation of where U = dx c /dt and v f = v cl (θ x ) = dx f /dt is the front (tip) velocity as a function of the dynamic contact angle at the tip, θ x , whose functional form is given implicitly by Eq. (8). The vertical motion of the ellipsoid and the horizontal displacement are not independent since so that there is an increase of the volume in the head as well as a variation of the contact angle, θ x . Upon writing dV /dt = V (z 0 )dz 0 /dt, Eqs. (13) and (15) lead to whose integration allows a comparison with the time evolution of most of the parameters measured in the experiments. In order to find the values of the constants a and c of the ellipsoid, we consider Eq. (11) at the tip position (y = 0), and its slope (contact angle) where For given values of x 0 , θ x , and h 0 (= c − z 0 ), Eqs. (18) where ξ = 2h 0 /x 0 must satisfy the condition ξ < tan θ x . Here, we consider x 0 , θ x , and h 0 at t = 0. On the other hand, to calculate b we consider Eq. (11) at the point of the contact line where the width of the head is maximum (x = x c ), and its slope (contact angle) where By using Eqs. (18)- (20) and (22)-(24), we define the following ratios: Therefore, in order to determine b in terms of a, any of these equalities could be used. Since θ x and θ y cannot be measured simultaneously, we discard the ratio of their tangents. Moreover, even if one could attempt to produce an yz-plane view, any measurement of θ y would be blocked by the drops of the previous ruptures. Instead, we consider the ratio y 0 /x 0 , which can be measured from top view observations. From experiments with three different widths (w b , w c and w d = 0.128 cm) and four ruptures for each width, we obtain an average value of the ratio x 0 /y 0 , as φ = 0.599 ± 0.08. Finally, we obtain the function v f (θ x ) inverting θ(v cl ) evaluated at the tip, see Eq. (8), and solve Eq. (17) numerically in time. Thus, we have where θ x (x 0 ) is given by Eq. (19), and x 0,i = x 0 (t = 0). Once x 0 versus t is found, we can calculate the contact angle θ x (x 0 (t)), the maximum thickness h 0 (t), and the tip position where the integral I(x 0 (t)) is We show in Fig. 6 a comparison of the time evolution as predicted by the model with the experimental data. As a first attempt, we determine the constants a, b, and c by taking the values of x 0 , h 0 , and θ x corresponding to the first frame capture after the breakup (t = 0). Thus, we obtain the curves shown as dashed magenta lines in Fig. 6 (t i = t 1 ). In spite of the simplicity of the model, based on a rough approximation of the free surface shape (assumed as an ellipsoidal cap), the agreement with experiments is remarkably good, except for the maximum thickness of the head, h 0 (t). Note, however, that immediately after the breakup the shape of the head is strongly changing from a spike to a rounded contour (see Fig. 8). Thus, during a short time of the order of a few seconds, the head does not achieve an ellipsoidal shape, and then, a too early selection of the initial values of x 0 , h 0 , and θ x can lead to an inaccurate determination of a, b, and c. Moreover, since these initial profiles are quite flat at the top, with h 0 being close to the filament thickness, it is difficult to determine x 0 with sufficient accuracy.
Better agreement with experiments could be expected if the initial values are taken from a later head whose shape has settled closer to that of an ellipsoidal cap. Therefore, we use instead the second frame (≈ 20 s later) as starting time, and obtain new results from the model (see solid magenta lines in Fig. 6, t i = t 2 ). We observe that the agreement with experiments is greatly improved for x f (t) and θ x (t), and even the h 0 (t) curve better approximates the experimental data.
Another interesting feature of the flow that can be extracted from the model is the maximum displacement of the receding tip, x f,max . In fact, according Eq. (8), the tip stops when θ x reaches the value θ 0 . Thus, the maximum value of x 0 is given by Eq. (19) in the form which yields a maximum value for the tip position, x f,max , as given by Eqs. (27) and (28) for x 0 = x 0,max . The values of x f,max for the first four heads in Fig. 6 are shown in Table I. We note that the ratio F = x f,max /w is practically constant in experiments. The corresponding theoretical F as obtained from our simple model is F theo = 3.82, which is in good agreement with the average experimental value. This result is in accordance with one of the features of the model, namely that all the distance parameters in Eq. (27), such as x 0 , h 0 , a, b and c, are proportional to w.  The use of the model for top view experiments is a bit more involved due to the lack of information about the contact angle at the tip, which is essential to the dynamics of the system. However, taking into account the proportionality of the coefficients with w (see Table I), it is possible to infer the values of the semi-axes of the ellipsoid by using the information from a single side view experiment, such as the case above. Thus, by multiplying the previous values of a, b, and c by the ratio w a /w b , we find the model results shown in Fig. 9. Even with this rather indirect method to determine the parameters, the model is able to predict reasonably well the evolution of the tip position, x f (t). Instead, it fails to properly describe w max (t), and yields higher values than expected. This result is consistent with the fact that h 0 (t) is underestimated in the side view case (see Fig 6(b)), since it implies an increase of w max to balance the mass flow swept by the receding tip. Note that the top view experiments do not provide any out of plane information, and therefore the input to the model dynamics is really scarce, so that the fact that it can pretty well describe at least x f (t) is remarkable.
These results suggest that the shape of the head can be roughly considered an ellipsoidal cap during most the receding stage. However, our model does not include any information about the transversal wetting, so it does not have the necessary ingredients to accurately predict the evolution of its width, w max (t), and maximum thickness, h 0 (t). On the other hand, the experimental evidence shows that by considering x 0 , h 0 , and θ x at successive times, the corresponding values of a, b, and c vary with time. In fact, while a slightly oscillates around a typical value (so that one could assume an average constant), b and c steadily increase with time. Thus, the differences in h 0 (t) observed between the model of constant semi-axes and the experiments can be explained in terms of the increasing length of transverse semi-axes of the actual approximating ellipsoid.

A. Basic equations and method
We numerically simulate the evolution of one end of the filament, and assume that the other one is so far away that it does not affect the other. Thus, we consider that at t = 0 the filament starts at x = 0 with the shape of a cylindrical cap of length L 0 , width w ( L 0 ), and with transverse equilibrium contact angle, θ a , along both parallel contact lines. The fluid region around the origin (x 0, |y| w/2) remains practically at rest, while the main flow develops far away from there, close to the end region (x L 0 ). Therefore, the boundary condition at x = 0 is that for a no-slip wall.
In order to emulate the actual filament, which ends at a rounded shape due to the breakup process that took place before (t < 0), we here assume that this shaped end can be approximated by an additional cap of length x 0 at the end of the cylinder (x = L). Thus, the fluid domain is composed by a cylindrical cap of length L 0 plus an ellipsoidal cap of length x 0 . The length x 0 is chosen in such a way that the filament tip has a given contact angle θ x,i , resulting from the breakup stage (see Fig. 11(a)).
The time evolution of this liquid filament of total length L = L 0 + x 0 is obtained by numerically solving the dimensionless Navier-Stokes equation where the last term stands for the gravity force. Here, the scales for the position x = (x, y, z), time t, velocity v = (u, v, w), and pressure p are the capillary length a c , t c = µa c /γ,γ/µ, and γ/a c , respectively. Therefore, the Laplace number is La = ργa c /µ 2 . In our experiments we have a c = 1.49 mm and La = 0.006, so that inertial effects are practically irrelevant. The x and y-axes are assigned along and across the original filament, respectively. Besides, the normal stress at the free surface accounts for the Laplace pressure in the form wheren = (n x , n y , n z ) andτ are the versors standing for the normal and tangential directions to the free surface. Since, the surrounding fluid (e.g., air) is passive, we assume that the tangential stress is zero, i.e. Σ τ = 0. As regards to the boundary condition at the contact line, the dynamic contact angle, θ, is given by the dimensionless contact line velocity, Ca = µv cl /γ, according to the CV+B model in the form (see Eq. (3)) where Ca 0 = µv 0 /γ, andˆ = /a c . The contact line velocity is calculated from the velocity field as where (N x , N y ) = (n x , n y )/ n 2 x + n 2 y is the versor normal to the contact line. Note that this condition introduces a high nonlinearity to the problem, since the solution itself, namely the velocity field at z = 0, yields the corresponding contact angle.
We use a Finite Element technique in a domain which deforms with the moving fluid interface by using the Arbitrary Lagrangian-Eulerian (ALE) formulation [33][34][35][36]. The interface displacement is smoothly propagated throughout the domain mesh using the Winslow smoothing algorithm [37,38]. The main advantage of this technique compared to others such as the Level Set or Phase Field techniques is that the fluid interface is and remains sharp [39]. The main drawback, on the other hand, is that the mesh connectivity must remain the same, which precludes the modeling of situations for which the topology might change (rupture of the filament). The default mesh used throughout is unstructured, and typically has 3 × 10 4 triangular elements (linear elements for both velocity and pressure). The mesh nodes are constrained to the plane of the boundary they belong to for all but the free surface.

B. Comparison with experiments
In order to compare the numerical results with the experiments, we perform simulations with the wettability parameters as given by Eq. (9). Although the calculation cannot be continued beyond the breakup moment, the simulation is useful enough to account for the dewetting and wetting processes between breakups. Thus, the initial condition used here tries to emulate the scenario just after an actual breakup, and the resulting fluid motion is simulated till just before the next breakup. Fig. 11 shows the time evolution of the filament of width w a . Here, we use L 0 = 8 to define the cylinder length (we have checked that for L 0 > 6, the results do not depend on L at all), and θ x,i = 25 • at the tip, which determines the ellipsoidal shape that emulates the filament end (or nose) just after the rupture.
The solid lines in Fig. 12 show the time dependence of the tip position, x f , tip contact angle, θ x , and maximum head thickness, h 0 , as given by the numerical simulations using the values of the constants θ 0 and v max from Eqs. (7) and (1), respectively. Regarding Fig. 12(a), the simulations yield smaller values of both the slope of x f (t) for early times and the maximum tip position for late times respect to the experimental data. Note also, that the numerical values of θ x and h 0 lie below the experimental points (see Figs. 12(b) and (c)).
Clearly, the degree of agreement between simulations and experiments is highly dependent on the parameters in Eq. (32). The weakest part of this fitting process is the inclusion of a single point to account for an important part of the dewetting region of this constitutive relation, namely v max , which is obtained by performing a quite different experiment than that used to generate the rest of the data (see Fig. 4). Since there is an uncertainty in the value of v max due to the experimental error (of the order of 7%), we explore the possibility that its inaccuracy can be a source of these differences. Moreover, although there is consensus on the hypothesis that the contact angle is zero at v max [28], one cannot guarantee that this is actually the case. If a different contact angle corresponds to this velocity, then the resulting θ(v cl ) curve allows both higher angles and dewetting velocities. The simulations show that considering a 7% increase in v max strongly improves the comparison of the slope of x f (t) for early times, while the saturation value is not practically not affected. On the other hand, there is an additional uncertainty on how the value of θ 0 is determined, since the only requirement is that it falls inside the hysteresis range, namely (θ a , θ d ). In this case, the numerical results show that the saturation value of x f (t) favorably compares with the experiments if θ 0 is increased a few degrees with respect to the value given in Eq. (7).
Thus, in order to better approximate the experimental data, we perform a simulation by slightly modifying both parameters to the values v max = 1.07v max and θ = 50.57 • , and using consistent modified values of the coefficients in Eq. (9) to fit the experimental data in Fig. 4: Γ = 95.4553, = 0.0008302a c = 1.24 × 10 −4 cm, v 0 = 6.2121 × 10 −7 cm/s. The results of the simulation are shown by the dashed lines in Fig. 12. We observe a clear improvement of the comparison for all three quantities shown in that figure, in particular for the evolution of the contact angle at the tip, θ x .
A similar tendency to approach the experimental data is evident in the thickness profiles when using the alternative values v max and θ 0 (see Fig. 13). In order to focus our attention on the shape of the profiles, both the experimental and numerical ones have been shifted so that the position of the tip is always at x = 0. Although, there are some departures both near the maximum and the neck, the latter are the most significant ones.
Since the differences between experiments and numerics are mainly quantitative and reasonable small, the latter are still very useful to learn about the hydrodynamics of the dewetting process at the filament end. In fact, the numerical results in Fig. 14(a) show the behaviour of the velocity and pressure fields in the head and neck regions for times close to breakup at the intersection of the free surface with the plane y = 0. At the tip, the pressure is relatively high and induces the dewetting motion (u < 0). Inside the head, the pressure is practically uniform and increases towards the neck, where it reaches a maximum value. This peak drives a flow out of the neck region towards the head (u < 0) and the filament (u > 0). Fig. 14(b) shows the flow field in a horizontal plane close to the substrate. This vectorial field velocity complements the picture already described and indicates that the change of behaviour from receding to advancing occurs at different values of x whether one observes the problem at the top of the free surface or near the contact line.
In order to describe the evolution of the shape of the head, we resort now to an experiment observed from the top, see Fig. 5(b), which allows us to register the evolution of the whole contact line. In Fig. 15, we compare the tip position, x f (t), the head width, w head , and the neck width, w neck , with the experimental data. The comparison of x f (t) yield similar conclusions with respect to the effects of using the modified values (v max , θ 0 ) to those of the side view case in Fig. 12(a). However, unlike the previous case we observe here that the slope of x f (t) for early times is a bit smaller than the experimental results. Also, Fig. 15(b) shows a greater departure of the widths of both head and neck from the experiments, as well as very slight improvements when using the primed parameters. We believe that all these differences can be attributed to the fact that the ellipsoidal shape assumed at t = 0 for the filament end after breakup might not be the best choice to describe the longitudinal thickness profile. Unfortunately we are not able to observe simultaneously both the thickness profile and the shape of the footprint to completely elucidate this issue. In Fig. 16, we compare the contours of the footprints at different times by using (v max , θ 0 ) (solid lines) and (v max , θ 0 ) (dashed lines), as we have done for the thickness profiles in Fig. 13. Here, we also focus our attention on the shape of these footprints, so that both the experimental and numerical contours have been shifted so that the position of the tip is always at x = 0. Unlike the comparison made for the thickness profiles, the modified parameters yield here a slight improvement only in the prediction of the neck shape, while the shape of the rest of the contact line is barely affected.
Interestingly, the simulations also show the existence of an angular sector that embraces the evolution of the head contact line, similarly to what is observed in the experiments (see Fig. 7(b)). This feature is shown in Fig. 17(a). In fact, the red thick lines in that figure correspond to the envelope of the parametric family of curves (x(s, t), y(s, t)) representing the contact line in the simulations. Here, the envelope is calculated by numerically solving the condition: Note this curve can be approximated by a straight line, as the one observed in the experiments and whose slope turned out to be independent of w. Additional simulations show that this is actually the case, though the constant α is 7.69 • here. In accordance with Fig. 16, this angle is not modified by using the set of parameters (v max , θ 0 ). The difference with the experimental value of α is consistent with the departures of x f and w head mentioned above.
On the other hand, we calculate the curve joining the contact line points where the local contact angle is equal to θ 0 , i.e. where the contact line velocity is zero. The interesting fact is that this other curve turns out to be identical to the above envelope. Thus, we conclude that the existence of the envelope is related to the points of the contact line which do not have normal velocity, but only tangential one. This effect is illustrated by the black points in Fig. 17(b), where the contact line velocity field is shown by vectors. Therefore, these points divide the portion of the contact line that is receding (dewetting) from that portion which is advancing (wetting) in the normal direction to the contact line. This result allows to understand the nature of the envelope that embraces the evolving footprint of the head and provides the physical underlying reason of its existence.

VI. SUMMARY AND CONCLUSIONS
In this paper we are concerned with the dynamics of the axial dewetting of a liquid filament laying on a horizontal substrate under partial wetting conditions. We observe that the retraction motion as well as the consequent transverse spreading strongly depend on the features of the contact angle hysteresis. For this reason, we study in detail the hysteresis cycles both under static and dynamic conditions. In the first case, we are able to determine the hysteresis range as well as the critical angles at which the drop contact line position must change in order to obtain a new static shape as its volume is varied (see Section II A). For the dynamic case, we obtain the constitutive relationship between the contact line velocity and the dynamic contact angle, and include this relation within a hybrid model that takes into account both a hydrodynamic approach and a molecular kinetic description (see Section II B). These results play a fundamental role in understanding the dynamics of the retracting filament, which is the main goal of this work.
The measurements of different characteristics of the evolving filament, such as tip position, contact angle at the tip, head thickness, head width, and neck width as well as contact line shape and vertical thickness profile, allow to perform a complete description of the whole evolution, so that the proposed wettability model can be thoroughly tested. In order to understand the main features of this flow, we develop a heuristic model to account for the time evolution of these parameters. The main outcome of this model is a remarkably good prediction of the retraction motion of the tip and the contact angle there. However, it is unable to properly describe the time evolution of the width of the head, since it does not consider the properties of the contact angle hysteresis in the transverse direction.  Therefore, we also perform a comparison of the experimental data with the results of the numerical solution of the full Navier-Stokes equations, which include the contact line velocity versus the contact angle as given by the hybrid constitutive relationship resulting from Eqs. (3) and (5). The simulations are able to properly describe the flow when using this hybrid law with the coefficients as measured in the characterization of the hysteresis cycle. It should be pointed out that the use of this function θ(v cl ) is an essential ingredient to adequately account for the experimental data. Naturally, one could be tempted to use a more usual and simpler relation as given by Eq. (3) with θ m = θ 0 , known as Cox-Voinov model. However, the simulations show that using this law instead of the hybrid one leads to results that do not compare well with the experiments (see gray lines in Figs. 12 and 15). This fact strongly emphasizes the need to properly model the hysteretic effects when describing the dynamics of a receding filament, which simultaneously involves wetting and dewetting motions of the contact line. As discussed in Section V, the main differences with the simulations observed in some experimentally measured quantities for the retracting filament end can be attributed to an incomplete knowledge of its shape just after the breakup. One of the most interesting facts revealed by the simulations is the physical property of the contact line points which belong to the envelope that encloses the footprints at different times. It is shown that these points have no normal velocity of the contact line, and correspond to the characteristic contact angle of the hysteresis curve, θ 0 . This envelope is also observed in the experiments as two straight lines, in agreement with the simulations.
Our analysis has focused on the retraction mechanism of the filament prior to its breakup. However, the final stages of this process can be affected by the rupture dynamics, whose understanding certainly merits future work. in the numerical simulations.