Breakup Dynamics of Capillary Bridges on Hydrophobic Stripes

The breakup dynamics of a capillary bridge on a hydrophobic stripe between two hydrophilic stripes is studied experimentally and numerically. The capillary bridge is formed from an evaporating water droplet wetting three neighboring stripes of a chemically patterned surface. The simulations are based on the Volume-of-Fluid (VOF) method implemented in Free Surface 3D (FS3D). In order to construct physically realistic initial data for the VOF simulation, Surface Evolver is employed to calculate an initial configuration consistent with experiments. Numerical instabilities at the contact line are reduced by a novel adaptation of the Navier-slip boundary condition. By considering the breakup process in phase space, the breakup dynamics can be evaluated without the uncertainty in determining the precise breakup time. It is found that within an intermediate inviscid regime, the breakup dynamics follows a $t^{2/3}$-scaling, indicating that the breakup process is dominated by the balance of inertial and capillary forces. For smaller bridge widths, the breakup velocity reaches a plateau, which is due to viscous forces becoming more important. In the final stage of breakup, the capillary bridge forms a liquid thread that breaks up consistent with the Rayleigh-Plateau instability. The critical wavelength is identical to the distance between the tips of two liquid cones between which the thread is arranged. The existence of satellite droplets in a regular pattern indicates that the primary breakup process is followed by self-similar secondary breakups.


Introduction
Wetting of patterned surfaces is omni-present in nature. The Lotus effect [1] or the fog harvesting of the Stenocara desert beetle in the Namib Desert [2] are only two examples. Also when it comes to engineering applications like ink-jet printing [3] or water harvesting [4], the understanding of wetting behavior on (chemically) patterned surfaces is of crucial importance. In the present study the focus lies on chemically patterned striped surfaces which are wetted by water droplets with radii of the order of magnitude of the stripe width. For such kind of systems, static wetting behavior has been studied by using energy minimization techniques [5][6][7], lattice Boltzmann simulations [8,9] and experiments [10]. Also inertia-driven spreading [9] and the splitting [11][12][13] of impinging droplets on such kind of surfaces has been investigated in some detail. When it comes to evaporation, simulations using the phase field [14] method have been performed more recently. Hartmann and Hardt [15] showed that an evaporating droplet wetting two hydrophilic stripes, with a hydrophobic one in between them, is stable as long as the pressure inside the liquid forming the capillary bridge above the hydrophobic stripe can be balanced in the liquid above the two hydrophilic stripes. This is the case until a certain width of the capillary bridge is reached. When the width of the bridge decreases further due to evaporation, it breaks up. This breakup process is the subject of the present article. Qualitatively, this process shows similarities to the collapse of a soap film between two circular rings, as it was investigated by Chen and Steen [16] using numerical calculations. The authors found two breakup regimes, which are both dominated by the balance between capillary and inertial forces. The change in regime is due to a geometric transition during breakup. The first scaling they found is a classical result for the evolution of the minimum width d of a "free liquid bridge" (i.e. without contact to a substrate) obtained from dimensional analysis [17] d(τ ) ∝ C στ 2 ρ where σ is the surface tension coefficient, ρ is the fluid density and τ = t 0 − t is the time t before the bridge breaks up at t 0 . While the prefactor C has been believed to be a universal constant with a value close to 0.7 for quite some time (see, e.g., [18]), recent studies show that this is not the case [19]. Instead, the experimentally observable value of C may depend on the initial state of the system as well as on viscous effects [19].
When viscous forces come into play, the appropriate dimensionless number characterizing the breakup of a free capillary bridge of a liquid is the Ohnesorge number which is the square root of the ratio of the viscous length scale µ 2 /(ρσ) (with µ being the viscosity) and the system's characteristic length scale which, according to Li and Sprittles [20], is given by the minimum width d. For water in air, the Ohnesorge number is given as Oh(d) ≈ 890 · 10 −6 kg/(ms) 997 kg/(m 3 ) · 72 · 10 −3 N/m ·d · 10 −6 m where d =d · µm. Clearly, the limiting cases Oh → 0 and Oh → ∞ characterize the inviscid and viscous breakup regime, respectively. The inviscid breakup regime is already discussed above for the soap film. Within the viscous regime, it is viscous forces that play the dominant role and d scales linearly with time τ , d ∼ Kτ . The authors of [20] presented a phase diagram, in which, besides the above-mentioned two regimes, also a third regime, the viscous inertial regime for intermediate values of Oh, can be identified in the parameter space that is spanned by d and Oh. Here, as well as in the viscous regime, d ∼ Kτ is valid but with a different constant K. For water as fluid at usual lab conditions and typical bridge widths as they appear in the present study, Oh is always larger than 6 · 10 −3 and, according to [20], all of the above-mentioned regimes are possible.
Note that [16] as well as [20] deal with a free capillary bridge without substrate contact. This is not the case in the present problem, where the capillary bridge is in contact with the substrate. For liquids in contact with a substrate, Bostwick and Steen [21] provided a review, where they focused on the stability of constrained capillary surfaces in general. The same authors theoretically investigated the instability of static rivulets in [22] and considered varicose (symmetric) and sinuous (anti-symmetric) modes for pinned and free contact lines. For symmetric modes, they confirmed the results by Davis [23]: For all contact angles, breakup can happen in the case of a free contact line, while for pinned contact lines, contact angles must be greater than 90 • . Other studies exist for unpinned (e.g. [24,25]) and pinned (e.g. [26] ) contact lines, as well. In case of anti-symmetric modes, a static rivulet is always stable in the case of a pinned contact line and unstable only for contact angles larger than 90 • .
The above-mentioned articles deal with the stability of rivulets. To the best knowledge of the authors, no investigations of the breakup and dynamics of capillary bridges in contact with surfaces have been published so far. The aim of this study is to examine the corresponding breakup dynamics compare it to a "free" liquid bridge.
these wafers are silanized in a low-pressure chemical vapor deposition process, for which 1H,1H,2H,2H-Perfluorodecyltrichlorosilane (PFDTS, CAS: 78560-44-8, abcr GmbH, Germany) is used. In the following step, the photo resist is rinsed off using acetone and isopropanol. Subsequently, the substrate is dried in a nitrogen stream. This leads to a pattern of hydrophilic and hydrophobic stripes. The substrates are then stored until experiments are performed. Different hydrophilic stripe widths w phil and different ratios α of hydrophobic and hydrophilic stripe widths are used, where α = w phob /w phil . Experiments are performed by placing de-ionized water droplets (18.2 MΩ·cm at 25 • C, Milli-Q) with a pipette onto the stripe pattern. The volume is chosen in a way that two hydrophilic stripes with one hydrophobic stripe in between them are wetted. The water then evaporates in the lab environment (temperature ∼25 • C) until a critical width of the bridge on the hydrophobic stripe is reached. The critical width marks the transition to a configuration beyond which no stationary wetting state of the droplet exists (see [15] for a discussion of the stability of the liquid bridge). We did not control the lab humidity since the final breakup process occurs at a small time scale (which we expect to be in the order of the capillary time scale τ , see equation (1)), and therefore the influence of evaporation can be neglected. The final time span (∼ 0.13 s) before breakup, i.e. when the capillary bridge on the hydrophobic stripe decays, is recorded with a high-speed camera (Photron FASTCAM SA-1.1) in top view mode. In order to obtain the values necessary for fixing the initial condition for the numerical simulations, e.g. the wetted length and the contact angle on the hydrophilic stripe in the moment of breakup, a second high-speed camera (Photron FASTCAM SA-X2) is used, which records the droplet from the side. The cameras are synchronized via dedicated software (Photron FASTCAM Viewer) and triggered externally. A frame rate of 75,000 frames per second (fps) is used. Each camera is connected to a macro objective (Navitar-12X). Illumination is performed via backlight in side view and co-axially in top view using two cold-light sources (VOLPI intra LED 5). A silicon wafer that acts as a mirror is placed below the glass substrate in order to achieve good illumination conditions and record the breakup process at a sufficient frame rate and magnification. We will also report results on the very last stages of bridge breakup. In these cases, the capillary bridge is recorded in bottom view mode with the Photron camera attached to a microscopy body. A more detailed description of this experimental setup, used materials and substrate preparation steps can be found in a previously published article [15]. Contact angles on the silanized glass wafers are measured using the evaporation method, in which the receding contact angle can be measured in the constant-contact angle mode during evaporation [27].

Stationary States from Energy Minimization
To be able to simulate the process of droplet breakup, a precise geometric description of liquid distribution at the beginning of the simulation is indispensable. Hence, the interface geometry has to be extracted from the experimental observations and transformed into a phase fraction field that can be processed by FS3D. For this purpose, a triangulated surface is produced using Surface Evolver [28], a tool for calculating minimal surfaces by minimizing the free energy functional subject to a prescribed volume. In the present case, the interfacial tensions σ i are constant, while the specific interfacial areas |A i | can change to minimize the free energy E. The subscripts sl, s and l denote solid-liquid, solid and liquid. As input data, the wetted length of the hydrophilic stripe as well as the contact angle on it, as observed in the experiments, is used. These values are summarized in table 1. Furthermore, the contact angle  Figure 1: Initial geometry for the numerical simulations obtained from Surface Evolver (blue). Due to symmetry, only 1/4th of the droplet is calculated. Red areas denote hydrophobic, green areas hydrophilic stripes. The wetted length l phil , the contact angle on the hydrophilic stripe Θ phil , the hydrophilic stripe width w phil and α, the ratio between the hydrophobic and the hydrophilic stripe width, are held constant during the evolution. Used values can be found in TABLE 1.
on the hydrophobic stripe is assumed to be the receding contact angle on a silanized glass wafer. All of these three parameters are kept constant during the whole minimization procedure. Due to symmetry considerations and to reduce calculation time, only a quarter droplet is calculated, with one mirror plane being parallel to the stripes in the middle of the hydrophobic stripe and the other one being perpendicular to the stripes, cutting the hydrophobic bridge in the middle. At these mirror walls, the contact angle is set to be 90 • . The corresponding surface is indicated with blue color in FIG. 1. To find a geometry with minimum stable bridge width, the algorithm is as follows: 1. Start calculation with a certain volume.
2. Evolve surface and calculate c v,7 , the coefficient of variation of the energy of the last 7 iterations, as well as the bridge width d. During the evolution of the surface, the mesh is successively refined until the cell size is smaller than 1/100 · (w phil + w phob ). Convergence can only be achieved when this is true. In steps 3 and 4, the volume is de-or increased following the principle of nested intervals. The result of the above-described algorithm is the triangulated surface mesh with the smallest stable volume. This geometry is exported as a triangulated surface mesh in STL format and subsequently transformed into a volume-fraction field for FS3D.

Continuum Mechanical Model
The continuum mechanical model is based on the incompressible two-phase Navier Stokes equations. The conservation equations for momentum and mass in the bulk phases read To simplify the model, it is assumed that no mass is transferred across the liquid-vapor interface. This can be expected to be a good approximation on the timescale of the breakup process which is much smaller than the evaporation timescale. Together with the assumption of no slip at the liquid-vapor interface, we obtain continuity of the velocity field, i.e. v = 0 on Σ(t), x − hn Σ ) denotes the jump of a physical quantity across the interface. The interfacial transmission condition for momentum reads where p is the pressure and S = η(∇v + ∇v T ) is the viscous stress tensor. The right-hand side is the surface tension force with σ > 0 the surface tension coefficient and κ the mean curvature. For simplicity, the surface tension is assumed to be constant. The motion of the interface is coupled to the bulk flow through the kinematic boundary condition where V Σ is the interface normal velocity. Formally, the latter condition can be reformulated as the advection equation where is the indicator function for the liquid phase. The transport equation (8) forms the basis for the Volumeof-Fluid interface capturing method (see Section 3.3). In order to regularize the moving contact line singularity [29], the Navier slip boundary condition is applied for the velocity at the solid wall, i.e.
v · n ∂Ω = 0, βv + (Sn ∂Ω ) = 0 on ∂Ω \ Γ(t), where L = η/β > 0 is the so-called slip length. The wettability of the solid is modeled through the contact-angle boundary condition where the contact angle may be a function of position to model regions of different wettability. Here, for simplicity, we apply a fixed contact angle in both regions, i.e.

Volume-of-Fluid Discretization
In the present study, the two-phase flow solver Free Surface 3D (FS3D), originally developed by Rieber [30], is employed to solve the incompressible two-phase Navier Stokes equations. Here we only discuss the most important aspects of the discretization and refer to [30][31][32] for details.
The two-phase Navier Stokes equations in the Continuum Surface Force (CSF) formulation [33] are discretized using the finite-volume approach on a fixed Cartesian grid. Within the CSF formulation, the effect of surface tension is modeled as a singular source term in the Navier Stokes equation and the equations (5) and (7) are replaced by where δ Σ denotes the surface delta distribution on Σ(t).
The phase indicator function χ is replaced by the discrete volume fraction in each computational cell, which is used to track the location of the interface. Integration of (8) over a control volume V 0 yields the transport equation for the volume fraction, i.e.
Within the geometrical Volume-of-Fluid method, the flux on the right-hand side of (12) is approximated with geometrical methods. Based on a reconstructed interface geometry, the numerical fluxes for the volume fraction are computed using an operator-splitting method [34], which decomposes the transport problem in a series of one-dimensional transport problems along the coordinate axis. The interface is locally reconstructed as a plane in each cell (also known as PLIC) [35], where the Youngs method [36] is used to estimate the interface normal vector based on the volume fraction field. At the contact line we employ a three-dimensional variant of the Boundary Youngs reconstruction method developed in [37].
Within the CSF formulation, the two-phase flow is treated as a single fluid where the density ρ is volume averaged according to Here the subscripts l and g refer to the liquid and gas phases respectively. The averaging of the viscosity in interface cells involves both arithmetic and harmonic averaging as discussed in [38] (see also [39]).
The time integration is based on an explicit Euler method, where the pressure-velocity coupling is realized by Chorin's projection method which leads to an elliptic equation for the pressure to perform a projection onto the space of solenoidal velocity fields. Note that the grids for velocity and pressure are staggered to enhance the stability of the method [40].
The surface-tension force is discretized with the balanced CSF method introduced in [41]. A height-function representation of the interface is constructed in order to approximate the mean curvature. It has been demonstrated that this method is able to significantly reduce spurious currents at the interface away from the boundary. Following the approach of [42,43], the height function is also used to indirectly enforce the contact-angle boundary condition. The idea is to linearly extrapolate the height function at the contact line into a ghost cell layer, where the slope of the extrapolated interface is determined by the prescribed contact angle. As a result, the approximated value of the mean curvature is altered leading to a "numerical force" that drives the interface towards the desired contact angle. A drawback of that method is that it may create spurious currents at the contact line.
In order to ensure stability of the numerical method, the time step is chosen according to the stability criterion ∆t = min{(∆t) σ , (∆t) µ , (∆t) u } with the timescales (see [44] for a similar criterion) given by No-slip and Navier Slip boundary conditions: To allow for a motion of the contact line one can either make use of the numerical slip inherent to the method or prescribe the (staggered) Navier slip condition. The numerical slip is a property of the advection algorithm, which uses face-centered values to transport the volume fraction field. In fact, the velocity boundary condition is indirectly enforced within the Finite Volume method using the concept of "ghost cells".
(i) In the classical approach, the no-slip boundary condition is enforced by extrapolating the tangential velocity field into the layer of ghost cells next to the physical boundary, such that the velocity interpolates to zero exactly at the boundary (see FIG. 3(a)). This velocity gradient is subsequently counteracted by the discrete realization of the viscous forces applied next to the boundary. This approach enforces the no-slip condition only in the limit of the mesh size going to zero. Therefore, the face-centered velocity at the boundary cell layer may be non-zero leading to a motion of the contact line. However, the extent of numerical slip decreases with increasing mesh resolution, leading to a significant mesh dependence of the solution. This numerical effect has been first described in the context of VOF methods in [45].
(ii) The Navier slip boundary condition is enforced following the same approach by setting a different velocity in the ghost-cell layer, such that the tangential velocity interpolates to zero at a fixed distance L > 0 from the physical boundary (see FIG. 3(b)). It has been demonstrated in the literature that this may lead to mesh-convergent results if the resolution of the computational mesh is well below the slip length L (see, e.g., [46], [47]). However, the physically expected slip length is on the scale of nanometers [48]. This length scale is not accessible with the present numerical method without massive computational costs. Moreover, note that for L > 0 the "counter velocity" in the ghost cell is always less than or equal to the counter velocity in the case of the numerical realization of no-slip.
"Staggered" slip: In the numerical simulations using the no-slip and Navier slip boundary conditions, we observed an instability of the numerical method at the contact line. As a result, unphysical velocities are created that trigger capillary waves propagating along the droplet. The origin of this instability is not yet fully understood. However, by introducing a modification of the Navier slip boundary condition we are able to dampen these unphysical velocities at the contact line significantly. The idea is to apply the slip length with respect to a virtual boundary located at the position of the face-centered velocities (see FIG. 3(c)). This allows to apply a slip length below the mesh resolution leading to larger "counter velocities" than with the standard Navier slip or no-slip conditions. Through (discrete) viscous forces this leads to a stronger damping or numerical dissipation at the contact line. In fact, for a staggered slip lengthL smaller than half the cell size, the staggered slip can be interpreted as standard Navier slip with a negative slip length −∆x/2 > L > 0, whileL = ∆x/2 corresponds to L = 0 (numerical slip). Note that a small staggered slip length may lead to very small numerical timesteps since the velocity in the ghost cells may become very large. Therefore, one cannot choose arbitrary smallL.   Conversion of a surface mesh from Surface Evolver to VOF data using OpenFOAM: A volume mesh is produced in OpenFOAM that has the same mesh resolution as the mesh in FS3D. The surface mesh obtained from Surface Evolver is then intersected with the OpenFOAM mesh using the Surface-Cell Mesh Intersection (SMCI) algorithm [49]. The SMCI algorithm resolves each vertex of the experimental surface mesh as a node of an octree data structure, which subdivides the bounding-box of the experimental surface mesh into octants. The octree spatial subdivision is then used to find mesh points closest to the triangles of the experimental surface mesh with logarithmic complexity in terms of the number of vertices in the surface mesh. Signed distances between volume mesh points and nearest surface mesh triangles are computed: positive inside the fluid phase, negative outside. This sign information is then propagated throughout the volume mesh resulting in the inside-outside information for each cell of the volume mesh. This defines each cell of the mesh as inside (χ = 1), outside (χ = 0), or intersected 0 < χ < 1. Centers of the intersected cells are then used together with their sizes to form enclosing spheres. The enclosing sphere is then used to query the octree structure of the surface mesh for all triangles of the surface mesh that are intersected with the enclosing sphere. Finally, each intersected cell is triangulated into tetrahedra to handle possible non-convexity of cells, which are then intersected with the triangles in the enclosing sphere. This returns a highly accurate geometrical value of 0 < χ < 1 for the intersected cells. The accuracy of 0 < χ < 1 of course depends on the resolution of the surface mesh delivered from Surface Evolver.
Numerical setup for FS3D: In order to save computational resources, we make use of the symmetry of the problem and simulate only a quarter of the droplet, while applying symmetry boundary conditions at the respective symmetry boundaries. At the outer boundaries of the domain, we apply no-slip for the velocity, which turns out to be irrelevant for the breakup dynamics (compared to, e.g., fixed pressure outflow boundary conditions). The computational domain has a size of 1000 µm ×2000 µm ×500 µm and is subdivided into an equidistant mesh of 2N × 4N × N cells, where N varies between 48 and 128. Therefore, the maximum resolution of the mesh is 500 µm /128 ≈ 3.9 µm in each direction. Unfortunately, adaptive mesh-refinement is currently not available for our solver, which limits the applicability of the numerical simulations to problems with a certain minimum length scale (see Section 4). The physical parameters used in the numerical simulation are listed in Table 2

Qualitative Comparison between Experiments and Numerics
Since the focus of the present study is on the dynamics of the capillary bridge on the hydrophobic stripe, this part of the droplet is shown in more detail in FIG. 4. It shows the bridge at three different instants in time τ in experiment and simulation for w phil = 500 µm and α = 1. In the top row, the black and dark gray regions represent the liquid (with a reflection of the light in the middle of the capillary bridge), while the liquid is white in the bottom row. The images of both rows also show some liquid that is wetting the hydrophilic stripes in the top and bottom part of each frame. Qualitatively, in both experiment and simulation the capillary bridge develops from a catenoid type to a narrower shape. Then a liquid thread is formed which is getting constricted at two points. The shapes, as well as the minimal widths of the bridge d are in good agreement between experiment and simulation at each time. In the final instants before breakup (τ = 0.02667 ms), the minimal diameter of the bridge, as well as the bulge which forms in the middle, seems to be larger in the simulations. This is probably a numerical issue since d approaches the minimum cell size. As a result, the bridge is no longer properly resolved. Therefore, the formation of the liquid thread in the middle of the hydrophobic stripe, its bulging and the subsequent constriction at the upper and lower end of the thread is not part of this comparison and will be addressed in Section 4.4. During breakup, the parameter R, which is the maximum width of the capillary bridge at the boundary between the hydrophobic and the hydrophilic stripe, stays constant over time. The absolute value of R differs slightly between simulation and experiment. This might be due to different contact angles on the hydrophobic stripe: Owing to unavoidable contaminations, the receding contact angle in the final stages of evaporation can be significantly smaller than the measured receding contact angle (see e.g. [27] and the supplementary of [15]).

Breakup Dynamics in Phase Space
Following the literature, the breakup dynamics is usually described via the minimum width d as a function of the time τ before the breakup event, i.e. τ = t 0 − t, where t 0 is the breakup time. For the inviscid breakup of a free capillary bridge, it can be shown by means of an asymptotic analysis that this function follows a power law (1). However, the precise time of the breakup event is hard to determine both in experiments (due to finite spatial and temporal resolution) and in the simulation. Note that the choice of t 0 can have a large effect on the effective exponent that is extracted from the data. In [50] it has been reported that the same set of data appears to approach power laws with an exponent ranging from 0.6 to 0.8 depending on the choice of t 0 . In the Volume-of-Fluid simulation, the actual breakup is usually mesh-dependent since it is ultimately performed by the interface reconstruction algorithm. Moreover, since the breakup process involves very small length scales, it cannot be fully resolved by the numerics. Therefore, the numerical results can only be considered meaningful down to a certain length scale determined by the computational mesh. In the present study, this length scale is approximately 40 µm. Consequently, the breakup time cannot be extracted from the numerics in a meaningful way without extrapolating the data.
We therefore propose a different approach to describe the breakup dynamics which is independent of the choice of the breakup time (as for the method introduced in [19]). Given the minimum width of the capillary bridge as a function of physical time t, we consider the breakup speed, i.e. the time derivative V = −ḋ as a function of the minimum width itself, i.e. V (d) = −ḋ(d). The latter quantity is well-defined since the minimum width d is a monotonically decreasing function (over the time scale of the breakup process). Obviously, the function V (d) is invariant with respect to shifts in the time coordinate.
For the power law (1) describing the inviscid regime we have Hence the power law (1) translates to More generally, one can easily show the relation for an arbitrary power law (α, β > 0). For the material parameters of water in air (T = 298 K, p = 100 kPa, see Table 2), the relations (1) and (13) take the form where τ =τ ms and d =d µm. In the present paper, these relations will be referred to as "inviscid theory".
Note that the above method requires to differentiate potentially noisy data with respect to time. This issue has been addressed by filtering out high-frequency oscillations in the experimental values of V (d) by locally fitting a straight line to the data (using six neighboring points). Despite this difficulty, the method allows studying the breakup dynamics in detail without the uncertainty in choosing t 0 . In the following, we will report both V (d) and d(t 0 − t) for completeness.

Quantitative Comparison between Experiments and Simulations
We first consider the case α = 1. FIG. 5 shows the experimental data for the breakup speed as a function of the minimum width together with the numerical data for different choices of the hydrophobic contact angle θ phob ∈ {80 • , 90 • , 102 • , 110 • , 120 • } and a fixed slip length ofL = 500 nm. Note that the measured hydrophobic contact angle from the experiment is 102 • . After a short initial phase, the speed from the numerical simulation appears to be consistent with the inviscid breakup theory (V ∝ d −1/2 ), where the constant C may depend on θ phob . In particular, the cases θ phob = 110 • and θ phob = 120 • show a significantly smaller value of the constant C. The breakup speed is monotonically decreasing with θ phob . The best agreement with the experimental data is found for θ phob = 102 • (as expected from experiment) and with the inviscid theory for C ≈ 0.82 (in the region 80 µm d 120 µm). At approximately 70 µm, the speed in the experiment reaches an almost constant value of approximately 570 µm /ms before it starts to decrease further towards the breakup. This trend, which can also be observed in the numerical data, can be understood by a transition to a regime where viscous effects become important and (1) is no longer valid (see slip-length study below). The breakup speed for small minimum widths d 60 µm depends only weakly on θ phob .   6 shows a comparison for the minimum bridge width as a function of time between the simulation (θ phob = 102 • ,L = 500 nm) and the experiment. Here the breakup time t 0 for the simulation has been chosen to best fit the inviscid theory while ignoring numerical data below 40 µm. The experimental value of the breakup time is determined from the pictures taken by the high speed camera. The first image where the capillary bridge is pinched off defines the breakup time t 0 . The horizontal error bar is defined by the used recording frequency (75,000 Hz) which corresponds to ∆t = 1.33 · 10 −5 s. In the vertical direction, the error bar represents the standard deviation obtained from 5 experiments. The numerical results are in good quantitative agreement to the experimental data and the inviscid theory for minimum widths in the region between about 80 µm and 140 µm. Below this length scale, viscous effects become important, and both the experiment and the simulation start to deviate from the inviscid theory (see below). Moreover, the dynamics at the onset of breakup (d 140 µm) in both the experiment and the simulation is significantly slowed down compared to the inviscid theory. In the experiments, this is because the bridge width is decreasing due to evaporation. This continues until the critical width is reached, indicated by the maximum of the simulation curve in FIG. 6 corresponding to the configuration computed by Surface Evolver. To quantify the influence of the numerical discretization, FIG. 7(a) shows numerical data for the breakup speed for different mesh resolutions (8N 3 computational cells where N ∈ {48, 64, 96, 128}). The results appear to be reasonably mesh-independent for d larger than approximately 40 µm, where the simulation on the coarsest mesh (N = 48, ∆x ≈ 10.4 µm) breaks down. The finest mesh (N = 128, ∆x ≈ 3.9 µm) delivers reasonable results down to approximately 10 µm. In each case this corresponds to only 4 computational cells within the bridge width. Since mesh convergence cannot be assured on this scale, the numerical data below approximately 40 µm should not be used to draw physical conclusions in a quantitative way.
As shown in FIG. 7(b), whereL is varied from 250 nm to 2 µm, the choice of the (staggered) slip length shows only a minor influence on the breakup dynamics in the early stage of breakup. This behavior is to be expected in the inviscid regime (low Oh), where viscous dissipation due to (partial) slip is irrelevant for the breakup process. A monotonic increase in the breakup speed withL is visible for a minimum bridge width below approximately 80 µm (Oh larger than 0.017). The experimental breakup speed for small minimum widths is slightly larger than the numerical value even for a very large slip length of 2 µm.
The breakup dynamics for α = 0.5 is reported in FIG. 8. After a short initial phase, a good quantitative agreement between the experimental data and the numerical data is found for θ phob = 102 • . In the region 60 µm d 100 µm, the dynamics is close to the prediction of the inviscid theory with C = 0.8. At a bridge width of approximately 60 µm, the experimental breakup speed reaches a plateau at approximately 630 µm /ms, which is captured accurately by the numerics. At approximately 50 µm, the speed in the experiment decreases further. Qualitatively, this effect is also present in the numerical simulation.
The results for α = 1.5 are displayed in FIG. 9. The experimental data in the inviscid regime (60 µm d 120 µm) are in good agreement with the inviscid theory with a slightly larger constant of C ≈ 0.87.    Hence a non-universality of the constant C, as discussed in [19], is indicated by the experimental data. The experimental value lies within the range that has been reported previously in the literature (0.45 C 0.97 according to [19]). The qualitative character of the experimental curve is similar to the previous cases. A plateau is reached at about 80 µm before the breakup speed decreases further, approaching breakup. The breakup speed data from the numerical simulations are in best agreement with the inviscid theory for a slightly smaller constant C ≈ 0.79 (compared to the cases α = 0.5 and α = 1.0). Moreover, the plateau in the numerical simulation is significantly lower than in the experiment. Therefore, the quantitative agreement between experimental and numerical data is weak in this case, while the qualitative behavior is similar.
The reasons for the differences may lie in the experiments as well as in the numerics. The numerical model assumes a fixed contact angle and neglects both contact angle hysteresis as well as dynamic effects. Following FIG. 9, the hydrophobic contact angle can influence the value of the constant C in equation (1) significantly. In particular, the breakup speed increases with decreasing contact angle. Moreover, the hydrophobic contact angle θ phob corresponds to the measured receding contact angle. It is well-known that for small values of the base diameter of an evaporating droplet, the measured contact angle is smaller than the receding angle in the constant contact angle mode (see, e.g., [27]). The critical bridge width lies within a range in which this effect might be relevant. Therefore, there remains an uncertainty in the contact angle that is implemented in the numerical model.

Rayleigh-Plateau Instability
In FIG. 10a   with diameter D is formed between the tips of two cones with opening angle β (τ = 0.04 ms). The thread then gets pinched at the tips of the cones (τ = 0.013 ms) and finally breaks up, leaving a primary droplet in the middle of the hydrophobic stripe (τ = −0.013 ms). Besides the primary droplet, some smaller secondary droplets can be seen. These droplets have their origin in a self-similar breakup process. This behavior is also found in the breakup of a soap film which is spanned between two circular rings [16]. The described breakup process exhibits two correlated features: First of all, cones, which in the twodimensional projection appear as triangles, are formed in each experiment (at all observed values of α) at a certain point in time. This is close to the point where a liquid thread is formed. As can be seen in FIG. 10b, these cones have unique opening angles β for each α. Their structure is independent of the hydrophilic stripe width. The value of β is highly reproducible between different experiments. FIG. 10b also contains results from numerical simulations for w phil = 500 µm. Since the minimal width of the capillary bridge is in the order of the mesh resolution in the final stages of the breakup process, the angles that can be inferred from the simulations are prone to error. Nevertheless, the trend is confirmed even though the values are about 10 • higher.
The liquid thread breaks up in a self-similar process. FIG. 10c shows all data of experiments performed with a hydrophilic stripe width of 500 µm for α = 0.5 − 1.5. The diameter D of the cylindrical thread is plotted versus the parameter λ, being defined as the distance between the tips of the cones in FIG. 10a. Different colors are used for different α's, and the mean values and standard deviations for each α are represented by ellipses around black crosses. Unfortunately, the numerical simulations cannot be used to support the experimental data because the length scale is below the smallest length scale that can be resolved numerically (in our case 40 µm). The values are compared with the classical Rayleigh-Plateau instability, which predicts that a cylindrical liquid jet of diameter D becomes unstable, since it is energetically more favorable to decompose into droplets. D is connected to the critical wavelength λ that belongs to the most unstable mode according to the formula From FIG. 10c it can be seen that, within the experimental error, the final breakup process on the substrate follows the classical Rayleigh-Plateau instability. The measured receding contact angle is 102 • , which is slightly higher than 90 • . For the case of 90 • a liquid jet with unpinned contact lines becomes unstable following the classical Rayleigh-Plateau result (see, e.g., Bostwick and Steen [21]). This also follows from symmetry considerations. Remarkably, the critical wavelength corresponds to the distance between the tips of the cones and the distance between the secondary satellite droplets. The reasonable agreement with the classical Rayleigh-Plateau theory, which is based on an inviscid fluid, also indicates that the influence of viscous stresses that shows up in the final stages of breakup, is not significant for the decay of the thread. The scattering of the data points of each experiment might be due to the uncertainty in the receding contact angle. This is already discussed at the end of section 4.3. In (a), the capillary bridge has evolved to a liquid cylinder between two cones. The final breakup happens at two pinch-off points leaving behind a primary droplet with smaller secondary droplets in the middle of the hydrophobic stripe. In (b), the cones opening angles β are evaluated for different α's and hydrophilic stripe widths. In (c), the diameter D of the cylindrical segment is plotted versus the measured wavelength λ for w phil = 500 µm and different α. Each point represents one experiment. These are condensed in the ellipses, representing the standard deviation around the mean value for each α. The dashed line corresponds to (14).
The mean critical diameter D of the final liquid thread and the critical wavelength increase with increasing α. A part of the data scatter visible in FIG. 10c is due to the spatial resolution of the camera. Since the final breakup process occurs at a relatively small time scale (≈0.04 ms in case of w phil =500 µm and α = 1, see FIG. 10a), a high frame rate is needed for image acquisition. Additionally, the region of interest is smaller than 1 mm. In order to achieve a good image quality with the given light source, image acquisition is limited to the given pixel resolution.

Conclusion
The breakup dynamics of a capillary bridge on a hydrophobic stripe, which forms during evaporation of a droplet, has been studied experimentally and numerically. The droplet wets two hydrophilic stripes, separated by one hydrophobic stripe. Different ratios of the hydrophobic and hydrophilic stripe width α were considered. The geometric Volume-of-Fluid method was employed to solve the three-dimensional two-phase Navier Stokes equations. In order to dampen spurious currents at the contact line, the Boundary Youngs interface reconstruction algorithm [37] was adapted to three dimensions, and a modification of the Navier Slip boundary condition was introduced ("staggered slip"). Physically realistic initial data for the VOF simulations were generated by directly importing geometrical data from Surface Evolver. The data import was achieved using an algorithm within OpenFOAM [49] that converts the surface mesh from Surface Evolver into a volume fraction field.
Different dynamic regimes were observed both in the experiments and the numerical simulations. After an initial stage, the breakup dynamics follows the well-known inviscid relation for a free capillary bridge, i.e.
where the constant C is, however, not universal. In the experiments both C ≈ 0.80 and C ≈ 0.87 were observed. The numerical simulations indicate that C may depend on the hydrophobic contact angle. Moreover, it is found that the dynamics in this regime is insensitive to the slip length, as can be expected from the inviscid limit. For smaller values of d (corresponding to Oh 0.02) viscous effects become relevant as indicated by a constant velocity regime (d ∝ τ ). This is to be expected in the viscous and viscous-inertial regime. For even smaller minimum widths of the capillary bridge, the breakup speed decreases further.
In the final observable stage of the breakup, a liquid thread (a cylinder cut by a plane) is formed between two cones. These have a well-defined opening angle that decreases with increasing α. The liquid thread breaks up in a Rayleigh-Plateau type instability. Remarkably, the critical wavelength of the breakup corresponds to the distance between the tips of the cones.
The introduced phase space picture of the breakup allows for a detailed study of the physics of the process. In particular, the ambiguity in the choice of the breakup time t 0 is removed.