# Thermally Driven Imbibition and Drainage Induced by Terraced Nanostructures

###### Abstract

Theoretical analysis and fully atomistic molecular dynamics simulations reveal a Brownian ratchet mechanism by which thermal fluctuations drive the net displacement of immiscible liquids confined in channels or pores with micro- or nanoscale dimensions. The thermally-driven displacement is induced by surface nanostructures with directional asymmetry and can occur against the direction of action of wetting or capillary forces. Mean displacement rates in molecular dynamics simulations are predicted via analytical solution of a Smoluchowski diffusion equation for the position probability density. The proposed physical mechanisms and derived analytical expressions can be applied to engineer surface nanostructures for controlling the dynamics of diverse wetting processes such as capillary filling, wicking, and imbibition in micro- or nanoscale systems.

## I Introduction: Nanoscale Wetting and Brownian Ratchets

Advances in nanofabrication and characterization techniques have enabled the engineering of nanostructured surfaces with geometric features as small as a few nanometers Rosei (2004); Tawfick et al. (2012); Checco et al. (2014). At nanoscales, the interplay between intermolecular forces, Brownian motion, and surface structure can give rise to complex interfacial phenomena that are challenging for the application of conventional, continuum-based and deterministic, models Schoch et al. (2008); Rauscher and Dietrich (2008); Bocquet and Charlaix (2010); Snoeijer and Andreotti (2013); Colosqui et al. (2013a). For example, nanoscale surface structures can induce energy barriers that lead to wetting processes governed by thermally-activated transitions between metastable states Davidovitch et al. (2005); Kaz et al. (2012); Kavousanakis et al. (2013); Colosqui et al. (2013b); Razavi et al. (2014); Duvivier et al. (2013). These thermally-activated transitions can result in directed transport of fluids and solutes when there is directional asymmetry of the energy barriers induced by the physicochemical structure of the confining surfaces Chinappi et al. (2006); Fu et al. (2007); Zuo et al. (2009); Sparreboom et al. (2010). Analogous mechanisms for rectification of thermal motion into directed transport underlie fundamental biological processes such as selective charge transport in ion channels or translocation of proteins across cellular membranes. Physical systems where thermal fluctuations are able to drive net directional motion, while performing work against “load” or resistance forces, are known as thermal ratchets or Brownian motors and have been extensively studied in the framework of statistical physics Astumian and Bier (1994); Reimann (2002); Hänggi and Marchesoni (2009).

Thermal ratchets can operate without thermal or chemical gradients provided that the system has not reached all necessary conditions for thermodynamic equilibrium Reimann (2002); Hänggi and Marchesoni (2009). A variety of novel nano/microfluidic devices perform as thermal ratchets to accomplish the handling, separation, and detection of diverse solutes (e.g., DNA, macromolecules, ionic species) and/or colloidal particles with an unprecedented precision Eijkel and Van Den Berg (2005); Han et al. (2008); Napoli et al. (2010); Bernate and Drazer (2012). These devices usually work with single-phase fluid solvents and must combine external electromagnetic fields, electrolyte solutes in proper concentration, and formation of electric double layers in order to induce energy landscapes with directional asymmetry (i.e., ratchet potentials). A different class of ratchet systems involving multiphase fluids has been demonstrated to produce “self-propulsion” of micro- or millimeter-sized droplets by combining micro/nanostructured surfaces, thermal/chemical gradients, and/or mechanical vibration Mahadevan et al. (2004); Prakash et al. (2008); Noblin et al. (2009); Mettu and Chaudhury (2010); Linke et al. (2006); Quéré (2013). Self-propulsion mechanisms in these multiphase systems are attributed to diverse dynamic phenomena, such as capillarity and contact angle hysteresis Noblin et al. (2009); Mettu and Chaudhury (2010), or evaporation flows and the Leidenfrost effect Linke et al. (2006); Quéré (2013), where thermal fluctuations play a secondary role.

There is a class of multiphase (two fluid) system that can perform as a thermal ratchet under isothermal and incompressible conditions, with or without the presence of electrolyte solutes and net surface charge. In this class of system the thermal ratchet mechanism is enabled by surface nanostructures that induce surface energy barriers with directional asymmetry. The particular configuration considered in this work, illustrated in Fig. 1a, consists of two macroscopically immiscible liquids (fluid-1 and fluid-2) confined in a slit-shaped channel or pore of height , length , and width . The surfaces confining the fluids are chemically homogeneous and neutrally charged. One of the surfaces has a terraced structure with regular tread length and riser height [cf. Fig. 1a] of nanometric dimensions. Similar terraced structures have been synthesized on crystalline substrates via diverse nanofabrication techniques such as wet etching, high-temperature annealing, and deposition of epitaxial films Rosei (2004); Kawasaki et al. (1994); Luttrell et al. (2009); Guisinger et al. (2009). The studied terraced structure with steps reduces the local height of the channel according to for (here, is the floor function and is the coordinate in the longitudinal direction). In the presence of an interface between two immiscible fluids, the interplay between thermal motion and surface energy barriers induced by the nanoscale structure can drive imbibition and filling/drainage processes in micro/nanoscale channels or pores for a range of wettability conditions unanticipated by conventional wetting models.

## Ii Theoretical Description: Thermally Driven Wetting

Analytical descriptions of thermally-driven wetting processes must consider that atoms or molecules in a liquid-fluid interface undergo thermal motion. We will analyze the case of unidirectional motion described by the average position of all atoms of the first fluid species (fluid-1) that lie at the front liquid-liquid interface [cf. Fig. 1a]. Adopting the average interface position to describe the dynamics of the confined molecular fluids implies projecting the (multidimensional) system energy landscape onto a one-dimensional profile along a “reaction coordinate” . The sequence of random displacements of the front interface position can be statistically described by the conditional probability density ; here, is the average interface position observed at a time . The stationary probability density is prescribed by the free energy profile and the thermal energy ; here, is the corresponding partition function, is the Boltzmann constant and is the system temperature. Assuming overdamped Brownian dynamics, the time evolution of the probability density is governed by the Smoluchowski diffusion equation

(1) |

where is the local friction coefficient or resistivity (i.e., the inverse of the mobility). For the studied conditions we consider a linear friction force that is mainly due to hydrodynamic effects and thus

(2) |

where is a drag coefficient, is the shear viscosity of the confined fluids, and is the total number terraces in the structure.
For analytical simplicity, we consider in Eq. 2 the case that both fluid-1 and fluid-2 have the same viscosity ; expressions for can be readily derived using similar hydrodynamic arguments.
Analytical estimates of the drag coefficient in Eq. 2 can be obtained by making simplifying assumptions about the modes of translation of the interface and the hydrodynamic velocity profiles induced.
A drag coefficient is obtained by naively assuming that a linear flow profile (i.e., Couette flow between two flat surfaces) develops after the sudden displacement of the contact line on one wall, while the contact line on the opposite wall remains stationary.

For isothermal and incompressible conditions and assuming sharp interfaces, the free energy profile is determined by surface energy contributions

(3) |

Here, is the interfacial energy of the liquid-liquid interface, is the Young contact angle measured on the fluid-1 phase, and is the area of the interface between the fluid-1 and solid phases. In the studied configuration the distance between the front and rear liquid-liquid interfaces is larger than the length of the terraced structure, which allows us to simplify the analysis. While the front interface moves inside each terrace there is a linear change in the interfacial area that results from conservation of volume, and thus we have

(4) |

for ().

The energy profile (Eq. 3) exhibits sharp energy increments when moving in the negative -direction across the edge of a terrace at position ().
Hence, periodic energy barriers at regular steps hinder backward random displacements as the liquid-liquid interface undergoes thermal motion along the terraced surface.
The time to cross over an energy barrier induced by nanoscale surface features can be predicted via Kramers theory of thermally-activated transitions, as documented in prior work Colosqui et al. (2013b); Razavi et al. (2014).
For the present analysis, it suffices to recognize that a mean time must elapse before observing a backward displacement of the interface over a terrace edge.
Therefore, over a time interval the presence of the energy barrier can be treated as a reflective boundary condition for Eq. 1 imposed at the edge of each terrace.
Furthermore, the friction coefficient (Eq. 2) and the wetting or capillary force (Eq. 4) remain constant within each of the terraces.
Hence, we have and for , which facilitates the analytical solution of Eq. 1.
In order to solve Eq. 1 within each terrace for which the edge lies at , it is convenient to introduce the dimensionless position
and time
.
In addition, we introduce the dimensionless force parameter given by the ratio of the work performed by the wetting force to the thermal energy within each terrace.
For a reflective boundary at and initial condition the analytical solution of Eq. 1 is given by Smoluchowski (1916); Lamm and Schulten (1983) , where

(5) |

is a dimensionless function determined by the variables and for a given force parameter . The mean displacement of the front liquid-liquid interface according to Eq. 5 is , where

(6) |

is the dimensionless mean displacement within the -th terrace. Taking the upper integration limit to infinity in Eq. 6 is an approximation valid for finite times . From Eq. 6 we can estimate the mean time

(7) |

at which (). Here, is the inverse of the dimensionless mean displacement function in Eq. 6 and is the time at which .

A few comments aobout the derived expressions are in order. Analytical integration in Eq. 6 is feasible and thus the dimensionless time to traverse the -th terrace () can be obtained by solving the implicit equation . The inverse function in Eq. 7 is uniquely determined by the dimensionless force parameter ; the function diverges at and decays for . An approximate explicit expression can be obtained from a first-order Taylor expansion of Eq. 5 about ; this approximation yields less than 10% error for . For a neutral contact angle the wetting force vanishes and thus for . Moreover, under neutral wetting conditions Eq. 6 predicts a mean displacement characteristic of diffusive processes, with an effective diffusivity . Notably, forward liquid displacements against wetting forces are expected to occur for contact angles provided that .

## Iii Molecular Dynamics Simulations

In order to verify analytical predictions and underlying assumptions we perform fully atomistic molecular dynamics (MD) to simulate the dynamics arising from couplings between Brownian motion, hydrodynamic effects, and wetting forces. The MD techniques employed in this work are extensively described in the literature Rapaport (1995); Frenkel and Smit (2002); Koplik et al. (1989); Koplik and Banavar (1995) and prior work by the authors Colosqui et al. (2013b); Razavi et al. (2014). The simulated system [see Fig. 1a] comprises two monatomic liquids (atomic species ) labeled as fluid-1 and fluid-2, and a crystalline solid (atomic species ). Atomic interactions are modeled by pairwise Lennard-Jones potentials , where is a characteristic interaction energy, is a symmetric attraction coefficient between species (), is the interatomic distance between any two atoms, and is the diameter of the repulsive core, which roughly correspond to the atomic diameter. Following conventional techniques for computational efficiency, atomic interactions are not computed for . The time scale of atomic displacements is and is in the order of picoseconds for simple molecular liquids; the atomic mass of the liquid species is set to be equal for both modeled liquids. The equations of motion are integrated using a fifth-order prediction-correction algorithm with a time step . A Nose-Hoover thermostat Rapaport (1995); Frenkel and Smit (2002); Nosé (1984) models the interaction with a thermal bath, regulating the system temperature to a prescribed value. In all MD simulations in this work the prescribed temperature is , the mean number density is and the shear viscosity is for both liquids. Solid atoms form a face-cubic-centered (fcc) lattice with uniform spacing ; typical values of for crystalline solids range between 0.2 and 0.5 nm. The set of values employed for the attraction coefficients (, , , and ) render two macroscopically immiscible fluids with a liquid-liquid interfacial tension and Young contact angles near neutral wetting conditions .

## Iv Results: Theory and Molecular Dynamics

The studied geometric configuration [cf. Fig. 1a] consists of a slit nanoscale channel of height = 30–60, length = 70, and width = 10; periodic boundary conditions are applied in the - and -direction. The bottom wall has different surface structures with = 5–11 terraces of length = 3–9 and riser height = 1–2. The volume confined above the terraced structure is . The fluid-1 phase fills a fraction of the volume as the front liquid interface moves within [cf. Fig. 1a]. In our MD simulations, the number density is constant and the filling fraction is directly computed from the number of atoms of fluid-1, , and fluid-2, , occupying the volume . At initialization, the volume is fully occupied by the fluid-2 phase and thus .

The filling fraction for six MD realizations
^{1}^{1}1Reported average quantities correspond to ensemble average over 10 to 12 MD simulations for the same initial macroscopic conditions and different initial conditions for the atomic velocities.
and the ensemble average are reported in Fig. 1b for one of the studied structures (, ) and neutral wetting conditions ().
As seen in Figs. 1b, the mean filling fraction uniformly increases with time while individual MD realizations exhibit rapid transitions between “long-lived” metastable states at specific values ().
The imbibition process occurs within a range of Young contact angles [cf. Fig 1c] for which capillary forces can be neutral or negative ().
Displacement beyond the last terrace edge, where , is prevented by a large energy barrier
.

Theoretical predictions from Eq. 7 for mean displacements, , and filling fractions, , are in close agreement with MD simulations reported in Figs. 2–3 when using = 4–5.5 in Eq. 2 for different structures. As seen in Fig. 2a, the mean filling rate is approximately constant under neutral wetting conditions for the different studied structures. Notably, mean displacements collapse to a unique curve when scaling by the terrace length and diffusion time as showed in Fig. 2b; here, is the resistivity in Eq. 2 for . As shown in Fig. 2c, increasing the channel height for a given length enhances the interface displacement rate by decreasing the local resistivity . The effect of varying the contact angle and thus the wetting force (Eq. 7) can be seen in Fig. 2d. For the mean displacement is enhanced by positive wetting forces . As predicted, positive displacements still can be observed for when ; the most adverse wetting force corresponds to (, ) in Fig. 2d.

Results reported in Fig. 2 confirm that a thermal ratchet mechanism mediated by surface energy barriers can lead to thermally-driven transport of immiscible liquids in slit channels or pores. This phenomenon is expected to occur provided that the time to cross over each energy barrier is larger than the time (Eq. 7) to traverse the -th terrace. Indeed, MD simulations for large terrace length and small riser height [see Fig 3a] report a decay in the mean displacement rate after a time . Increasing the energy barrier by increasing the riser height to prevents the observed decay in the displacement rate [see Fig 3a]. As illustrated in Fig 3b, the mean displacement rate is expected to vanish for , after which unbiased Brownian motion persists over a time . For Brownian motion with a mean square displacement causes the liquid-liquid interface to diffuse to the next terrace edge, after which thermal motion becomes biased, , and the cycle repeats.

## V Conclusions

Theoretical description in the framework of statistical physics predicts that surface nanostructures with directional asymmetry can induce nontrivial wetting processes that are beyond the reach of conventional continuum-based models (e.g., Lucas-Washburn equation). Fully-atomistic simulations showing close agreement with theoretical predictions document the thermally-induced displacement of immiscible fluids confined in a slit channel or pore with a nanoscale terraced structure. For the studied nanoscale systems under neutral wetting conditions, a water-air interface would exhibit mean displacement rates between 0.1 and 1 m/s. Moreover, the observed fluid displacement can oppose the action of wetting or capillary forces. Hence, the studied mechanism can induce the wetting and dewetting of nanostructured pores or capillaries under unexpected wettability conditions. The maximum contact angle for which thermally driven imbibition or drainage occurs against the action of capillary forces is largely determined by the geomtry of the channel and terraced structure.

The analytical approach presented in this work indicates the necessity to consider the thermal motion of liquid-fluid interfaces in order to predict novel phenomena and associated useful effects. The proposed thermal ratchet mechanism enabled by engineered surface nanostructures could enable novel nanofluidic devices for passive handling and separation. Micro/nanostructured surfaces for superhydrophobic or superoleophobic behavior could exploit the thermally driven drainage of micropores and cavities to enhance their performance. In the presence of electrolyte solutes, the studied thermal ratchet mechanism could be employed to direct the transport of charges, which can potentially enable microfluidic applications for energy harvesting.

###### Acknowledgements.

We thank A. Checco and M. Sbragaglia for useful discussions. We acknowledge support from the SEED Grant Program by Brookhaven National Laboratory (BNL) and Stony Brook University. This work employed computational resources at the BNL Center for Functional Nanomaterials supported by The U.S. DOE under Contract No. DE-SC0012704.## References

- Rosei (2004) F. Rosei, J. Phys. 16, S1373 (2004).
- Tawfick et al. (2012) S. Tawfick, M. De Volder, D. Copic, S. J. Park, C. R. Oliver, E. S. Polsen, M. J. Roberts, and A. J. Hart, Adv. Mater. 24, 1628 (2012).
- Checco et al. (2014) A. Checco, A. Rahman, and C. T. Black, Adv. Mater. 26, 886 (2014).
- Schoch et al. (2008) R. B. Schoch, J. Han, and P. Renaud, Rev. Mod. Phys. 80, 839 (2008).
- Rauscher and Dietrich (2008) M. Rauscher and S. Dietrich, Annu. Rev. Mater. Res. 38, 143 (2008).
- Bocquet and Charlaix (2010) L. Bocquet and E. Charlaix, Chem. Soc. Rev. 39, 1073 (2010).
- Snoeijer and Andreotti (2013) J. H. Snoeijer and B. Andreotti, Annu. Rev. Fluid Mech. 45, 269 (2013).
- Colosqui et al. (2013a) C. E. Colosqui, M. E. Kavousanakis, A. G. Papathanasiou, and I. G. Kevrekidis, Phys. Rev. E 87, 013302 (2013a).
- Davidovitch et al. (2005) B. Davidovitch, E. Moro, and H. A. Stone, Phys. Rev. Lett. 95, 244505 (2005).
- Kaz et al. (2012) D. M. Kaz, R. McGorty, M. Mani, M. P. Brenner, and V. N. Manoharan, Nat. Mater. 11, 138 (2012).
- Kavousanakis et al. (2013) M. E. Kavousanakis, C. E. Colosqui, and A. G. Papathanasiou, Colloids Surf. A 436, 309 (2013).
- Colosqui et al. (2013b) C. E. Colosqui, J. F. Morris, and J. Koplik, Phys. Rev. Lett. 111, 028302 (2013b).
- Razavi et al. (2014) S. Razavi, I. Kretzschmar, J. Koplik, and C. E. Colosqui, J. Chem. Phys. 140, 014904 (2014).
- Duvivier et al. (2013) D. Duvivier, T. D. Blake, and J. De Coninck, Langmuir 29, 10132 (2013).
- Chinappi et al. (2006) M. Chinappi, E. De Angelis, S. Melchionna, C. Casciola, S. Succi, and R. Piva, Phys. Rev. Lett. 97, 144509 (2006).
- Fu et al. (2007) J. Fu, R. B. Schoch, A. L. Stevens, S. R. Tannenbaum, and J. Han, Nat. Nanotechnol. 2, 121 (2007).
- Zuo et al. (2009) G. Zuo, R. Shen, S. Ma, and W. Guo, ACS Nano 4, 205 (2009).
- Sparreboom et al. (2010) W. Sparreboom, A. Van Den Berg, and J. Eijkel, New J. Phys. 12, 015004 (2010).
- Astumian and Bier (1994) R. D. Astumian and M. Bier, Phys. Rev. Lett. 72, 1766 (1994).
- Reimann (2002) P. Reimann, Phys. Rep. 361, 57 (2002).
- Hänggi and Marchesoni (2009) P. Hänggi and F. Marchesoni, Rev. Mod. Phys. 81, 387 (2009).
- Eijkel and Van Den Berg (2005) J. C. Eijkel and A. Van Den Berg, Microfluid. Nanofluidics 1, 249 (2005).
- Han et al. (2008) J. Han, J. Fu, and R. B. Schoch, Lab Chip 8, 23 (2008).
- Napoli et al. (2010) M. Napoli, J. Eijkel, and S. Pennathur, Lab Chip 10, 957 (2010).
- Bernate and Drazer (2012) J. A. Bernate and G. Drazer, Phys. Rev. Lett. 108, 214501 (2012).
- Mahadevan et al. (2004) L. Mahadevan, S. Daniel, and M. Chaudhury, Proc. Natl. Acad. Sci. 101, 23 (2004).
- Prakash et al. (2008) M. Prakash, D. Quéré, and J. W. Bush, Science 320, 931 (2008).
- Noblin et al. (2009) X. Noblin, R. Kofman, and F. Celestini, Phys. Rev. Lett. 102, 194504 (2009).
- Mettu and Chaudhury (2010) S. Mettu and M. K. Chaudhury, Langmuir 26, 8131 (2010).
- Linke et al. (2006) H. Linke, B. Alemán, L. Melling, M. Taormina, M. Francis, C. Dow-Hygelund, V. Narayanan, R. Taylor, and A. Stout, Phys. Rev. Lett. 96, 154502 (2006).
- Quéré (2013) D. Quéré, Annu. Rev. Fluid Mech. 45, 197 (2013).
- Kawasaki et al. (1994) M. Kawasaki, K. Takahashi, T. Maeda, R. Tsuchiya, M. Shinohara, O. Ishiyama, T. Yonezawa, M. Yoshimoto, and H. Koinuma, Science 266, 1540 (1994).
- Luttrell et al. (2009) T. Luttrell, W.-K. Li, X.-Q. Gong, and M. Batzill, Phys. Rev. Lett. 102, 166103 (2009).
- Guisinger et al. (2009) N. P. Guisinger, T. S. Santos, J. R. Guest, T.-Y. Chien, A. Bhattacharya, J. W. Freeland, and M. Bode, ACS Nano 3, 4132 (2009).
- Smoluchowski (1916) M. V. Smoluchowski, Z. Phys. 17, 557 (1916).
- Lamm and Schulten (1983) G. Lamm and K. Schulten, J. Chem. Phys. 78, 2713 (1983).
- Rapaport (1995) D. Rapaport, The art of molecular dynamics simulation, 2nd ed. (Cambridge University, New York, 1995).
- Frenkel and Smit (2002) D. Frenkel and B. Smit, Understanding Molecular Simulations (Academic Press, 2002).
- Koplik et al. (1989) J. Koplik, J. R. Banavar, and J. F. Willemsen, Phys. Fluids 1, 781 (1989).
- Koplik and Banavar (1995) J. Koplik and J. R. Banavar, Annu. Rev. Fluid Mech. 27, 257 (1995).
- Nosé (1984) S. Nosé, J. Chem. Phys. 81, 511 (1984).
- (42) Reported average quantities correspond to ensemble average over 10 to 12 MD simulations for the same initial macroscopic conditions and different initial conditions for the atomic velocities.