Continuum concepts in nanoscale capillary impregnation

The design of tailor-made nanofluidic devices requires an extension of macroscale hydrodynamic theories for capillary impregnation (CI). Large-scale molecular dynamics (MD) simulations of a simple capillary pump consisting of a nanoscale gold slit attached to a liquid propane reservoir reveal that distinct finite-size-effects impact CI on the nanoscale. A continuum theory is derived that captures these finite-size-effects: properties of a prewetting monolayer, a related Navier-slip, nanoscopic contact angles and wall-induced oscillatory pressure fluctuations enter as non-heuristic atomistic input into the derivation of extended lubrication equations that exactly reproduce the capillary rise dynamics and menisci from our MD simulations. It turns out that impregnation in bare nanochannels can be significantly accelerated by the strong slip induced by a spreading precursor. As expected, this effect is absent in micron sized channels, where our extended continuum theory predicts a capillary dynamics that is already insensitive to all nanoscopic details at the contact line.


Introduction
Capillary impregnation (CI) underlies many traditional microfluidic applications like sap rise in trees [1], wicking assays for virus detection [2], catalyst formation [3] or oil recovery [4]. Modern lab-on-chip applications [5] such as separation of biomolecules [6], drug delivery [7], single molecule detection [8] or fluidic computing [9,10] foster the interest in CI as a robust and efficient driving force pumping autonomous nano-and microfluidic devices [11]- [13]. A detailed theoretical description of the capillary dynamics of liquids in nanosized gaps is essential for tailoring such nano-and microcapillary pumps [14].
Also other fundamental questions concerning capillary bridges in atomic force microscopy [15], lubricant transport in nanotribology [16], liquids inside carbon nanotubes [17] or water transport in fuel cells [18] would benefit from a deeper understanding of nanoscale capillary impregnation (NCI). In this context, molecular dynamics (MD) [19]- [21] provides deep insights into the basic microscopic processes. Unfortunately, these atomistic simulations are still too time-consuming to be routinely used by scientific and technological practitioners. This raises the fundamental question: how far can established continuum concepts for a quantitative description of capillary phenomena be carried over to scales where finite-size effects are expected?
The continuum theory of CI dates back to Laplace [22] who discovered that the curvature κ of a meniscus determines the driving capillary pressure p c = γ κ. Here, γ denotes surface tension and κ depends on the contact angle θ e , defined as the inclination of the equilibrium liquid-vapor interface with the solid boundary. A simple argument using the balance of p c and the pressure p p to drive a Poiseuille channel flow results in the celebrated Washburn equation [23]. For instance, a liquid sheet of length L(t) in a slit of width b needs p p = −12ηLL/b 2 leading to L(t) = γ bt cos θ e /(3η) 3 with η the viscosity of the liquid. During rapid early stage impregnation, θ e in equation (1) must be replaced by an apparent dynamic contact angle θ d [24], since capillary numbers Ca =Lη/γ are large then. Stating phenomenological constitutive laws θ d (Ca) is at the heart of the multiscale description of dynamic wetting phenomena [25], since molecular details at the contact line become important. Even for macroscale CI, θ d (Ca) must be inferred from asymptotic matching of hydrodynamic solutions on larger length scales supplemented by ad hoc nanoscale contact line models (see [26] and references therein). For nanosystems the situation is even more complicated, since finite-size-effects [27] can violate the assumed separability into continuum and atomistic scales. And even if a separation ansatz is viable, it is unlikely that established heuristic contact line models describe NCI quantitatively. Precision on the nanoscale requires more accurate atomistically derived boundary conditions [28,29]. Such questions can only be answered by a direct quantitative comparison of large scale reference MD simulations with corresponding continuum predictions for specific nontrivial nanosystems [30]. Here, we chose the transient impregnation dynamics in a nanocapillary pump consisting of liquid propane in a gold slit pore as a realistic test case. It is demonstrated that a quantitative continuum description of NCI is possible. We identify fluid-wall slip, mass flows induced by nanoscopic precursor films, nanoscopic contact angles and wall-induced oscillatory pressure fluctuations as necessary atomistic information that can be quantified by separate auxiliary MD simulations. Based on these details, we derive lubrication equations (LEs) that suffice to predict the menisci and rise of the capillary column.

The atomistic model
The propane molecules in our large scale MD simulations were described within the united atom picture following [31] supplemented by angle-bending potentials according to [32]. The gold walls were modeled by an embedded atom method potential [33], whereas the interaction potential between the gold atoms and propane segments was taken from [34]. Note that these interatomic interactions are essentially the same as in a previous study of another nanocapillary phenomenon, namely the Rayleigh instability of a propane nanojet leaving a 6 nm gold nozzle [30].
Transient impregnation simulations (see figure 1(a) for the MD setup) are performed as follows. An initial propane reservoir containing 152 000 molecules is equilibrated to T p = 150 K for a period of 0.5 ns employing a Langevin thermostat [35]. The solid channel is formed by two gold walls. Each channel wall contains 70 000 gold atoms organized in two fixed layers and two inner layers that are also Langevin-thermostated to T p . The channel axis is oriented in the x-direction, its inner width is chosen as b = 20.4 nm in the z-direction and periodic boundary conditions are applied in the y-direction with 12.24 nm periodicity. At t = 0 ps, thermalization of the reservoir is switched off, the propane is allowed to enter the gold channel and the evolution of the capillary column is monitored for a further 3 ns.

Results of the transient impregnation simulations
In the following, we contrast two different simulations. One whose further impregnation proceeds on a bare gold surface (indexed as 'Au' and visualized with red color in the following) and another one with a complete propane monolayer precoverage ('Au/ML', blue). The strong propane-gold interaction results in θ e = 0 • for both simulations. For Au, the additional free energy release upon formation of new solid-liquid interface area [25] drives a thin precursor film spreading on the bare gold ( figure 1(a)), in contrast to Au/ML with an semi-infinite monolayer that permanently adheres to the walls. Indeed, if and how these different nanoscopic details influence the continuum description of NCI will be of central importance in this paper.
The upper graph of figure 1(b) shows the evolution of L(t), defined as the x-distance between pore entrance and center of the meniscus. For both prewetting conditions, a rapid pore filling is observed. Interestingly, the spreading monolayer significantly accelerates the impregnation. It is generally assumed that the presence of any prewetting film should favor the 'slip' of a contact line, and thus enhance the dynamics (see for example [36]). Here, the opposite fact is found, which is a very interesting novelty and therefore this paper will also lay a special focus on the explanation of this effect.
The filling dynamics can be divided into distinct stages. During the first 200 ps, the initially flat meniscus relaxes toward a curved shape leading to negative L(t). This time is in accord with the fundamental period τ 0 = ρ(b/2) 3 /γ of the capillary wave spectrum [37] 5 . The subsequent dynamics is dominated by the balance of capillary and inertial forces; a regime characterized by an almost linear rise of L(t), extending over a period of approximately τ 1 = ρ(b/2) 2 /η ≈ 0.2 ns, the time needed to establish Poiseuille flow [24]. Thus, after τ 0 + τ 1 ≈ 0.4 ns one expects a balance of capillarity and viscosity and therefore the argument underlying equation (1) should hold. Thus, the further discussion is restricted to t 0.5 ns.
The Washburn equation (1) greatly overestimates the simulated dynamics (compare green and blue curves in figure 1(b)). To explain this, we associate the dynamic contact angle θ d directly with the actual driving pressure, i.e. θ d is defined as the inclination angle of the channel walls with a circle segment fitted to the central region of the meniscus 6 . The evolution of θ d in figure 1(b) explains the failure of Washburn's equation by the fact that θ d θ e , since the curvature (i.e. driving pressure) is much smaller. Further, it becomes clear why propane on Au/ML is slower than on Au: θ d at equal times is consistently larger. Although this resembles the established continuum description [24], this is not sufficient to become quantitative for our specific examples.

Derivation of a generalized Washburn equation by a lubrication treatment of the transient impregnation
One would expect that the differences between Au and Au/ML are largely due to the difference in θ d (t). However, as will be seen soon, for Au, a position-dependent slip-length λ(x) impacts the overall flow profile directly. Spatial variations of the Navier slip length are frequently discussed and applied in the context of moving contact lines [38]- [41]. This λ(x) boundary condition at the gold wall can be introduced into a lubrication treatment of the Stokes equation [37]. The resulting parabolic velocity field obeys the Navier-slip law as well as additional zero-stress boundary conditions at the meniscus and at the channel center (exploiting the symmetry; green line labeled h(x) in figure 1(a)). The flow through an area 5 Separate MD simulations are performed to obtain surface tension γ = 0.0150 ± 0.0008 N m −1 and viscosity η = 0.00030 ± 0.00002 N s m −2 for the T p = 150 K, ρ = 0.6 g cm −3 propane model fluid. 6 For the transient impregnation simulations, the meniscus has been determined by dividing the slit into a stack of 50 slabs oriented along the channel axis. The location of the liquid-vacuum interface has been defined by the length of the liquid column in each slab. A circle segment has been fitted to the inner region of the meniscus covering 80% of the width of the slit. Note, that this fit is excellent for t > 0.5 ns, whereas there are minor deviations in the initial phase of impregnation. The dynamic contact angle has been defined as the angle between the circle segment and the mathematical plane spanned by the inner-most layer of gold atoms. This definition establishes a one-to-one correspondence between the curvature of the meniscus (i.e. the driving Laplace pressure) and the dynamic contact angle. In the case of the steady state plug simulations, a different method for the determination of the meniscus has been chosen. Here, extensive temporal averaging has been possible and therefore the location where the steady state density (see, e.g. inset of figure 3(c)) dropped to 50% of its bulk value has been used to define h(z). The resulting higher resolution allows for an inspection of the meniscus oscillations close to the gold walls. between the wall and the curve h(x) can be calculated as Note that for later purposes, an optional wall velocity U was introduced. Evaluating equation (3) for U = 0 and p = p c = −2γ cos θ d /b (realizing that Q(t) =Lb/2) yields a generalized Washburn equation, The notation θ d (Ca(t)) indicates that θ d at time t is related to a Ca set by the current velocitẏ L(t) of the liquid column.

Determination of constitutive laws from atomistic auxiliary simulations
The constitutive laws θ d (Ca) and also λ(x) are crucial for a quantitative prediction of L(t).
Fortunately, it turns out that both laws can be determined from auxiliary steady state simulations of a finite driven propane plug with a stationary meniscus (see inset of figure 2). In these simulations, a b = 20.4 nm gold slit is considered. Each wall of the slit consists of 156 000 atoms and extends 145 nm in the x-direction. The channel is filled with a 40 nm long propane plug preceded by precursor monolayers covering the walls of the empty part of the slit. For the Au case, the monolayers are allowed to spread freely across the surface, whereas for Au/ML propane molecules in the vicinity of the rightmost egdes of the slit channel are frozen in order to emulate an infinite extension of the monolayer to the right. In total, the fluidic part consists of 73 000 and 81 000 molecules in the Au and the Au/ML case, respectively.

7
In order to establish a situation with a constant capillary number we chose a setup as shown in figure 2. In the rest frame of the wall the liquid is driven by an expanding parabolic numerical membrane in the x-direction. The numerical membrane consists of a 4 nm long plug of initially frozen propane molecules which are driven by a predescribed external velocity depending on their z coordinate. During expansion of the membrane fluid in the center of the slit is displaced, whereas the first layer at the walls is allowed to stick. The expansion of the parabola is chosen such that the same amount of liquid is displaced as with a piston which is moved with a spatially constant velocity in the x-direction. In this way, a steady state situation is established and the liquid is driven by a constant capillary number. This setup is chosen in order to generate the Poiseuille flow characteristic for a semi-infinite plug. Note that first attempts using a piston instead of a parabolic membrane introduced considerable disturbances in the first and second monolayer adhering to the wall.
ForvariouscapillarynumbersCa=Uη/γ,theconstitutivelawθ d (Ca) is extracted from thestationarymenisci(seefootnote6).Infigure2,thedynamiccontactangleisplottedas the third power over the capillary number. Note, that the classical Cox-Voinov theory [42,43] predicts a finite slope for small capillary numbers -in clear contrast to our results. We find that our θ d (Ca) compare approximately to more recent results by Eggers and Stone [26] who derived logarithmic corrections at small Ca. Indeed the graphs can neatly be fitted to their functional form with fit parameters α = 5.97 and β = 0.55 for Au and α = 3.78 and β = 0.46 for Au/ML. Note that the values for α and β deviate somewhat from those derived in [26]. It would be interesting to understand the origin of these discrepancies. However, such an investigation would certainly go beyond the scope of this work that uses equation (5) merely as a convenient interpolation of the atomistically derived θ d (Ca). Let us instead focus on the main finding reported in figure 2. The dynamic contact angle θ d for Au is significantly smaller than for Au/ML. Of course, this in combination with equation (1) provides an immediate qualitative explanation of the faster filling in the bare pore (figure 1). However, before we return to the full quantitative analysis of L(t), we have to determine the other constitutive law λ(x). This will also elucidate the mechanisms underlying the observed relation θ Au d < θ Au/ML d . An inspection of the stationary fluid velocities u(x, z) reveals marked differences between Au/ML ( figure 3(a)) and Au ( figure 3(b)). For Au/ML, the liquid close to the walls is transported in the negative x-direction with almost constant velocity, whereas for bare gold the liquid creeps away from the contact line along the positive x-direction. The appearance of the velocity profiles is similar to continuum predictions [44]: a self-similar caterpillar motion persists up to the mathematical point where the meniscus hits the boundary and the fluid stresses diverge if regular wall-stick boundary conditions are assumed [38]. The density distribution ρ(x, z) depicted as 2D plot in figure 3(c) exhibits a pronounced molecular ordering. Surprisingly, this layering extends right to the contact line forming terraces that resist the high shear stresses at the turning point of the flow field. This unexpected structural robustness facilitates the definition of an x-dependent slip length: figures 3(c) and (d) display the x-component of the velocities u x (x, L1 − L3) in the first three molecular layers (L1, L2 and L3) for Au/ML and bare Au, respectively. The pronounced differences between the curves in both panels are a direct consequence of the different simulated downstream boundary conditions. In the Au/ML case, the L1 molecules are pushed into the meniscus by the sticking monolayer at the right-hand side of the channel. Therefore, the L1 curve in figure 3(c) is close to the U = −10 m s −1 wall velocity. On the other hand, L1 propane on the bare Au spreads in downstream direction leading to a reduced L1 velocity at the triple line ( figure 3(d)). Of course, these differences have a crucial influence on L2 and L3 velocities.
We now employ a finite difference approximation of the Navier slip law equation (2) in order to evaluate the slip length λ(x). The deep valley between L1 and L2 in the inset of figure 3(c) indicates that L1 is separated from the hydrodynamics of the bulk fluid. Therefore, the λ(x) values for Au/ML and Au were computed from the respective L2 and L3 velocities. The x-dependent slip length in figure 4 suggests that stick boundary conditions (i.e. λ = 0) can be assumed for Au/ML, whereas considerable slip with a strong increase of λ close to the contact line occurs for Au. For the following lubrication calculations, the slip law for the Au case is approximated by Here, m λ = 0.44 denotes the slope and λ(0) = 1.64 nm the slip length in the minimum of the capillary meniscus at x = 0.

A lubrication treatment of the steady state impregnation
In order to find the mechanism underlying the faster wetting of the bare gold, let us sharpen our discussion. Could we just start with the microscopic information contained in figures 3(c) and (d) and build continuum models for both prewetting conditions that explain the stationary menisci? We display their shapes h(x) for various capillary numbers in figure 5(a) and (b) for Au and Au/ML, respectively. Here, h(x) is defined as the 50% bulk density iso-contour in the stationary density profiles ρ(x, z) and therefore the molecular layering of the liquid is reflected in h(x). For small Ca, the lubrication treatment yields a relatively simple equation for h(x), as the driving capillary pressure p = p c = γ κ is now directly related to the local curvature κ = −h /(1 + h 2 ) 3/2 . For the stationary situation depicted in figure 2, the flow in equation (3) can be evaluated (in analogy to [45]) far away from the meniscus as Q = (U wrf 1L + U )h 1L , implying Here, h 1L = 0.43 nm denotes the height and Ca wrf 1L = U wrf 1L η γ the capillary number of the first monolayer, related with its speed U wrf 1L in the wall rest frame. Note that prewetting conditions are reflected in different λ(x) and Ca wrf 1L : for Au/ML Ca wrf 1L = 0 and for bare Au we found Ca wrf 1L ≈ 0.227. Strikingly, equation (6) provides a virtually quantitative description of the atomistic meniscus shapes for all capillary numbers Ca < 0.1 (see black dashed lines in figures 4(a) and (b) and note that a breakdown of the lubrication approximation is expected for Ca 0.1 [46]). The solution of equation (6) depends crucially on the starting height h 0 and slope h 0 (the differential equation has been solved by shooting). In figures 5(a) and (b), the dashed solution starts from the second monolayer (i.e. h 0 = h 2L = 0.86 nm) with h 0 = 0 corresponding to θ e of propane on gold. The choice h 0 was justified by an analysis of the MD menisci exhibiting almost vanishing slope of the second monolayer terrace (see green curve in figure 5(c)). Of course, another possible boundary condition to the LE could have been h 0 = h 1L and h 0 = 0. However, in this case the LE solution clearly failed to predict the MD menisci (see black dotted line in figure 5(c)). As mentioned before, the first monolayer is too strongly adsorbed to the wall. At the triple line it can be viewed as a separate effectively solid system.

Influence of the structural disjoining pressure
One should also mention that equation (6) fails to reproduce the oscillations of h(x) for x < 2 nm. Of course, the present LE carries no information about the underlying molecular layering mechanisms. However, simply augmenting the pressure term in equation (6) with an additional structural disjoining pressure [47] s (h) (i.e. p = p c + s ) removes these discrepancies. We calculated s (h) for all Ca < 0.1 in the Au/ML case from the integrated equation (6) s Note that the structural disjoining pressure is closely related to the solvation forces [48] acting in a liquid between two walls separated by a distance h + h 1L . The solvation pressure in liquid propane between two gold walls was obtained from MD simulations analogous to [48]. In order to be comparable with the structural disjoining pressure (the propane adsorbed to one gold wall only), the interaction between one wall and the fluid was reduced tenfold to mimic a free surface. These solvation forces (black dots in the the inset of figure 5(c)) have the same order of magnitude and the same oscillatory behavior as the structural disjoining pressure.
Inserting a s (h) that has been obtained from the inversion of the Au/ML Ca = 0.049 meniscus in the extended pressure for equation (6), we were able to reproduce atomistic menisci for other capillary numbers (compare for instance the black solid line with the green curve in figure 5(c) for Ca = 0.073).

Summary of the nanoscale effects and discussion of the observed speedup mechanism for the impregnation in the bare gold channel
In summary, our nanoscale capillary LE needs only four pieces of atomistic information for a quantitative description of the MD menisci: Including the structural disjoining pressure s improves the results but is not essential for the curvature in the middle of the channel. With (a)-(c) alone, we can infer the correct constitutive law θ d (Ca) and finally understand why θ Au d < θ Au/ML d and consequently why the wetting in the bare channel is faster. Starting from the Au/ML where U wrf 1L = 0 and λ(x) = 0, the influence of (a) and (b) can easily be studied by increasing first the monolayer velocity U wrf 1L to its maximum value of 11.4 m s −1 . This is then followed by an increase of λ(x) to its full strength (see figure 5(d)). Interestingly, the effect of U wrf 1L points in the 'wrong' direction, i.e. θ d increases. Only after the introduction of the full slip is this effect overcompensated and the smaller θ d of the bare Au case is recovered. Thus eventually, the faster impregnation in the bare Au pore is traced back to the pronounced slip at the contact line.
Let us cast this mechanism in a simple physical picture. The wall stick boundary conditions for the Au/ML case result in a strong drag acting on the wings of the meniscus leading to a decreased curvature (compared to the equilibrium meniscus) and therefore to a reduced impregnation speed. This has to be contrasted to the Au case where spreading monolayers leak out of the wings of the meniscus and decouples to some extent the bulk of the meniscus from the motion of the gold walls. This leads to significant wall slip that reduces the drag on the wings and brings the meniscus closer to its equilibrium shape eventually accelerating the impregnation.

Performance of the generalized Washburn equation
Now, we are able to address the final question: can our generalized Washburn equation (4) predict the atomistic L(t) curves from figure 1? If we insert the θ d (Ca) and λ(x) that emerged from our stationary simulations into equation (4), the transient L(t) for both the Au and the Au/ML case is reproduced within error bars (black solid lines in figure 1, integration runs for t 0.5 ns). Therefore, we can conclude that a quantitative continuum description of nanoscale CI based on atomistically derived boundary conditions was here indeed achieved for liquid propane. Let us stress once more the chain of arguments that lead to this conclusion. The comprehension runs both ways, down and up the scales: 1. the dynamically measured contact angle (in figure 1) describes impregnation correctly, within the generalized Washburn framework; 2. it can already be inferred from an auxiliary, quasi-stationary simulation (figure 2); 3. a close-up view of the MD directly at the gold wall (figure 3) reveals the active microscopic boundary conditions (slip-laws in figure 4); 4. if these are then combined with essentially a mass flow condition imposed by the precursor, one can reproduce the macroscopic dynamic contact angle fully within a lubrication framework (figure 5).

Final conclusions
Further studies for other solid/liquid combinations are now in order. Our work will be certainly a useful guideline for these future efforts being transferable to other low-viscous fluids, like water. For instance, traditional Washburn calculations significantly overestimated the experimental pump power of the smallest microcapillary pumps considered in [46]. This is in agreement with our transient simulations in figure 1. Therefore, we suggest that instead of the simple equation (1) an extended Washburn treatment (in the spirit of equation (4)) should be used for an accurate design of nanocapillary systems. In addition, our treatment gives a hint of when genuine nanoscale behavior sets in. Figure 5(e) displays how θ d (Ca = 0.07) varies with increasing slit width b, for Au (red dots) and Au/ML (blue squares). Starting at slit width as thick as 0.2 µm the variation in θ d becomes notable. This is an important point considering that channels as thin as 10-20 nm are nowadays routinely fabricated [6,12,49]. Only for larger slit widths the influence of the nanoscopic precursor (U wrf 1L and λ(x)) is considerably diminished and even a solution starting with θ e from the first monolayer (instead from the second) yields the correct macroscopic dynamic contact angle (see blue diamonds). All these nanoscale details which are important for finite-size-effects in nanochannels are already immaterial on the micron-scale, indicating for instance that slip is unlikely to have an influence on the performance of microcapillary systems.