## Abstract

Hydrophobic coatings, such as octadecyltrichlorosilane or n-alkyl monolayers, enhance the slippage of liquids on solid walls. For a given alkyl chain length, the main structural parameter for homogeneous coatings is the tilt angle between coating molecules and the surface normal. In this paper, *ab initio* calculations are used to calculate the equilibrium configuration of coating molecules, showing that the tilt angle easily changes from 0^{°} to 30^{°} depending on the specific head group binding the solid substrate. These values are used to set up classical molecular dynamics of water slippage over the coatings using different water models (Transferable Intermolecular Potential 3 point (TIP3P), Transferable Intermolecular Potential 4 point (TIP4P) and TIP4P/2005). The slippage is found to be robust with respect to the coating tilting, while a slight dependence on the water model is observed.

## 1. Introduction

The recent progress in micro and nano fluidic device fabrication demands a better comprehension of fluid–solid interactions at the nanoscale where the usual description in terms of Navier–Stokes equations with no-slip boundary conditions at the solid wall is unsuitable. There are two possible reasons for this failure, namely the breakdown of the continuum assumption and fluid slippage at the wall. For simple liquids, the continuum assumption is rather robust [1,2], and departures from continuum model predictions are observed only under strong confinement, when the device dimension approaches the molecular size. For instance, the so-called single file motion regime sets in for flows through nanopores only when the pore becomes of the order of three to four molecule diameters (≃1 nm) [1]. Hence, at least for simple liquids and with the exception of sub-nanometre size, failure of the classical description should be ascribed to the violation of the no-slip boundary condition at the solid wall.

Wall slippage is observed for both smooth and patterned surfaces, see [3] for a review, and it is usually described in terms of the Navier boundary condition *v*_{w}=*L*_{s} d*v*/d*z*, where *v*(*z*) is the velocity profile with *z* being the wall-normal coordinate and *v*_{w} the slip velocity, i.e. the relative fluid–solid velocity at the interface. The parameter *L*_{s} is indicated as the slip length. The slip length encompasses the relevant rheological properties of the interface needed to model the fluid as a continuum, e.g. [4] where liquid polymer melts flowing over polymer brushes are addressed by molecular dynamics (MD).

In this paper, we focus our attention on smooth surfaces coated with self-assembled monolayers (SAMs) of either octadecyltrichlorosilane (OTS) or 18-alkyl chains, as typical examples of promising technologies for wall functionalization. In principle, structural properties of the SAM, such as the tilt angle, could affect the slippage. The study and characterization of SAMs obtained by deposition of alkyl-silanes (such as OTS) or n-alkyl chains on Si or SiO_{2} surfaces have been undertaken mainly from the experimental point of view by using several techniques to evaluate both the film quality (coverage, surface roughness and temperature stability) and its structure (self-assembling properties, thickness, chain orientation and hydrophobic termination groups). The literature provides experimental data concerning structural properties of the deposited films obtained by, e.g. ellipsometry, infrared spectroscopy, X-ray photoelectron spectroscopy (XPS), reflectivity and absorption fine structure (near-edge X-ray absorption fine structure; NEXAFS) [5–9]. It is found that the structural properties of SAM films depend dramatically on the deposition method and chemical environment used to deposit it.

Much less is known about the actual effects that the SAM structure may have on the slippage properties of the functionalized wall. The purpose of the present paper is to describe an investigation recently undertaken in this direction. We first describe *ab initio* calculation simulations aimed at extracting the correct structural parameters. These parameters are then used in a classical MD simulation of a water flow over the SAM-coated surfaces to address the slippage properties of the functionalized wall. A careful examination of the experimental literature on SAM-functionalized Si provides indications concerning the tilt angle *θ*_{tilt} formed by the carbon chain axis and the normal to the Si surface. More specifically, from infrared spectroscopy, NEXAFS or equivalent techniques, the tilt angle for the OTS SAM on (111) Si substrate has been measured in the range [0^{°}≤*θ*_{tilt}≤10^{°}], while the carbon chain axis of n-alkyl SAM on (111) Si surface forms a *θ*_{tilt}≃25^{°} tilt angle [9]. A clear picture is, however, rather difficult to extract from the variety of experimental conditions and measurement techniques. Alternatively, we decided to rely on *ab initio* simulations. So far, to the author’s knowledge, no reliable atomistic model of OTS or n-alkyl SAM deposited on Si or SiO_{2} substrates has appeared in the literature, which mostly describes only in qualitative terms the surface chemical bonding and the OTS/n-alkyl chain structure. Here, preliminary results on reliable atomistic modelling of OTS and the n-alkyl-isolated molecule deposited on the hydrogenated (111) Si (H : Si) substrate by means of total energy *ab initio* calculations are discussed. These data are then used to set up classical MD simulations in order to estimate the slip length *L*_{s}. The measurement of *L*_{s} is performed following the procedure described in previous work [10]. Both cases of OTS and 18-alkyl are considered in order to assess the effect of the SAM orientation on slippage. To evaluate the robustness of the results to the model details, different classical water models are considered.

The paper is organized as follows. In §2, the *ab initio* calculations are presented. Section 3 reports a brief description of the classical MD set up for slip measurement, and the results for the different cases are considered. A comparison with experimental and numerical studies in the literature is discussed in §4.

## 2. Tilting angle: *ab initio* simulations

In this section, we report some preliminary results on reliable atomistic modelling concerning the configuration of two typical coating molecules. We deal with a single molecule of the SAM, either OTS or 18-alkyl, deposited on a hydrogenated (111) Si substrate by means of total energy *ab initio* calculations.

The study of the structural stability of SAMs of CH_{3}(CH_{2})_{17}SiCl_{3} (OTS) on Si (111) surface is performed by first-principles calculations based on density functional theory (DFT; [11]). The silicon surface has been modelled as a three-layer slab of 96 atoms with 32 Si atoms per layer. Periodic boundary conditions (PBCs) have been employed, and the computational cell is completed with a vacuum region of approximately 43 Å, corresponding to 18 Si atomic layers. All the dangling bonds of the Si slab have been passivated with H atoms. The total zero-temperature energy of the system is *E*[*n*(** r**),

*R*_{i}], and the ground state to be determined minimizes

*E*with respect to both the electron density

*n*(

**) and the ion positions**

*r*

*R*_{i}. The equilibrium configuration is found by a two-step iterative procedure. The first step consists of minimizing the energy with respect to the electron density

*n*(

**) for given positions of the ions. In the second step, the ion positions are upgraded, accounting for the local energy gradients ∂**

*r**E*[

*n*(

**),**

*r*

*R*_{i}]/∂

*R*_{i}, the so-called Hellmann–Feynman forces. In more detail, concerning the first step, the total energy of the system is treated, in the frame of the adiabatic approximation, as a functional of

*n*(

**) and is variational with respect the electron density itself at the ground state 2.1where**

*r**v*

_{ext}(

**,**

*r*

*R*_{i}),

*T*[

*n*(

**)] and**

*r**V*

_{ee}[

*n*(

**)] are the external potential owing to the ions, the electron kinetic energy and the electron–electron interaction energy containing the Coulomb repulsive energy and the exchange-correlation energy. The above variational problem is usually recast into an independent electron eigenvalue problem [11], and the corresponding so-called**

*r**Kohn–Sham*equations can be solved self-consistently to find the ground-state electron density and total energy.

The key factor affecting the accuracy of *ab initio* calculations is the exchange-correlation functional that is unknown *a priori*. Here, we have employed the generalized gradient approximation (GGA), using the Perdew–Burke–Ernzerhof (PBE) formula for the exchange-correlation potential and energy functionals [12]. The usual choice for the exchange-correlation potential for semiconductor bulk systems or surfaces (such as the Si surface employed here) is the local density approximation (LDA) that has been derived from a uniform electron gas adapted locally. This approximation, despite being quite rough, is satisfactory for semiconductors, but is not satisfactory for molecules and for systems where long-range interactions are involved. Moreover, it is a general opinion that the reaction energy calculated using LDA is quite distant from the experimental values, especially for molecules. Thus, several GGA schemes have been developed and checked, even though they have not always improved the LDA results. Among all these schemes, the PBE density functional is a good choice to obtain fair accuracy for systems ranging from molecules to solids; hence, it is fit for the present case study. Concerning the external potential, norm-conserving pseudo-potentials have been constructed following the Troullier–Martins scheme [13] that allows the inclusion of oxygen 3-*d* electrons in a soft way, thus not affecting the computational cost too much. The Kohn–Sham equations are solved using a plane-wave basis set with a cut-off of 30 Ry for the wave functions and a 120 Ry cut-off for the electronic charge density, while a (3×3×1) Mokhorst–Pack (MP) *k*-point grid [14] has been used for Brillouin zone sampling and integration. The above computational parameters have been chosen after several accurate convergence tests on the total energy of the Si surface for different *k*-point MP grids as reported in figure 1; subsequently, the energy cut-off has been resized to include the presence of an alkyl-silane molecule. Finally, the Broyden–Fletcher–Goldfarb–Shanno (BFGS; [15]) algorithm has been used to optimize the atomostic configurations. All the calculations have been performed by using the QUANTUM ESPRESSO package [16]. The convergence threshold for the total force has been set to 0.001 (atomic units) and a PBC dipole correction is used to handle the dipole electric field [17].

The above specified *ab initio* theoretical model shows that the OTS stable configuration is nearly normal to the Si surface (figure 1*b*), whereas the carbon chain of an 18-alkyl molecule adhering on a hydrogenated Si (111) surface is about 30^{°} off the surface normal (figure 1*c*).

## 3. Water slippage: classical molecular dynamics

Concerning slippage, the main relevant structural property of the coating is the chemical nature of the terminal group exposed to the liquid. However, for a given terminal group, certain geometrical features may, in principle, affect the liquid–solid interface. Beyond chain length, for uniform coatings, the tilt angle *θ* is expected to be the most important. The results of §2 suggest that *θ* strongly depends on the specific system, and it is extremely difficult to determine from experiments. The zero-temperature *ab initio* simulations discussed so far may not correspond to the conditions encountered at finite temperature. Nevertheless, their indications can be used as starting point to investigate the effect of *θ* on slippage.

Values of 0^{°} for OTS chains and 30^{°} for 18-alkyl chains are used to set up classical MD simulations of a Couette flow in the nano-channel between two planar walls (figure 2*a*). The lower wall is coated, while the upper one is kept in motion by a constant force. The slip length is obtained from the steady-state velocity profile. The simulations are performed at 300 K and 1 atm. For the level of detail needed here, a precise description of the head group (the one towards the solid wall), of the specific binding mechanism and of the solid crystal is not necessary. The relevant information is the hexagonal structure of the SAM and its surface density. Hence, the sample is modelled as a (111) Lennard-Jones (LJ) solid surface parallel to the O*xy* plane. Previous studies showed that, after equilibration, the molecules of the coating spontaneously achieve a tilt angle *θ*_{tilt}≃30^{°} [10]. This fact is exploited here to reproduce the inclination of the *ab initio* simulation for the 18-alkyl chain. The case with 0^{°} inclination is obtained by fixing the chain carbon atoms. Sensitivity tests to the water model were also performed, using Transferable Intermolecular Potential 3 Point (TIP3P), Transferable Intermolecular Potential 4 Point (TIP4P) [18] and TIP4P/2005 [19] models for water, see Chinappi & Casciola [10] for further details.

The mean velocity profile is shown in figure 2*b* for the TIP4P/2005 water model with 30^{°} tilt angle of the alkyl chains. For each water model, the slope is consistent with the respective viscosity, with TIP4P/2005 (*μ*=0.855 m Pa s^{−1}) getting closest to the experimental value of *μ*=0.896 m Pa s^{−1} [20]. Defining the wall surface as the average position of the methyl group carbon atom of the chains, *L*_{s} is obtained from the position where the velocity profile extrapolates to zero. The results are summarized in table 1. No difference is observed by changing forcing direction and intensity (data not shown). Figure 2*c* shows the number density profile within the fluid (*ρ*_{f}, *z*>0) and the solid (*ρ*_{s}, *z*<0), formed either by the LJ crystal (uncoated case, solid line) or by the SAM (coated surface, *θ*_{tilt}=30, dashed line). The hydrophobic behaviour of the LJ crystal appears from the first strong peak in the water density (*z*≃3.5 Å), which is instead much weaker for the hydrophilic SAM coating. Denoting the bulk values by a superscript ‘b’, this behaviour is quantified by the depletion layer thickness (last line of table 1),

## 4. Discussion

The tilt angle of coating molecule chains presents a wide variability in the dependence of the specific wall-binding mechanism, §2. Actually, molecules, such as OTS and 18-alkyl, differing only by the binding group, attain largely different tilt angles when deposited on a hydrogenated Si (111) surface. Their ground-state structures have been obtained by both fully relaxed BFGS optimization and constrained static calculations where the tetrahedral coordination of the sylane group (OTS) and of the entire carbon chain (18-alkyl) have been kept fixed. Concerning OTS, three –OH groups were supposed to replace the chlorine atoms originally present in the trichlorosylane group. The substitution reaction occurs in water solution where the OTS molecules can easily adhere to the Si surface through an oxygen bridge by eliminating a hydrogen molecule. The 18-alkyl ground-state structure has other equivalent orientations with respect to the (111) Si surface. The relevant transition energy barriers should then be measured to evaluate the temperature stability of the tilt angle. In addition, in the OTS case, the –OH group orientation with respect to the chain axis and the H:Si surface may affect the tilt angle. A deeper analysis of such structural problems is under way, and will be reported in the near future.

Despite the significant differences between the single molecule, zero-temperature simulations and the conditions found in actual coatings, i.e. finite temperature and interaction with neighbouring chains, these results suggest that the tilt angle may change from case to case, potentially affecting the slippage. In fact, our all-atom classical MD results show no relevant difference in slip length between cases with normal and tilted chain coatings (compare case C and F (TIP3P) and E and G (TIP4P/2005) in table 1). This behaviour is confirmed by addressing the depletion layer thickness, which globally quantifies the hydrophobicity of the surface, which, within statistical accuracy, is the same for the two tilt angles. A statistically significant, although slight, variation is instead observed in the classical simulations when changing the model for the water molecule (simulations C, D and E). In particular, the recently developed TIP4P/2005 model gives a slip length *L*_{s}≃13 Å, significantly larger than the TIP3P value *L*_{s}≃6 Åand a slightly larger depletion layer thickness. As a general conclusion drawn from the present and from previous [10,21] MD simulations, the intrinsic slip for liquid water over a smooth surface coated with methyl-terminating SAMs is in the range 0.5–2 *nm*. This result is robust under changes in both the SAM structure (tilt angle) and the water model. It is one order of magnitude smaller than the most credited experimental value [22]. In our opinion, this discrepancy could hardly be ascribed to inaccuracies of the classical MD force field or to details of the specific simulation procedure.

This turns our attention to other possible explanations, like the presence of gas nanobubbles trapped in nanoscale defects. Several studies on the structure of OTS coatings suggest that perfectly smooth coatings such as those addressed in ideal MD simulations do not exist in practice. Cottin-Bizonne *et al.* [22] and Joseph & Tabeling [23] provide estimates of a few nanometres for the height of typical asperities in their coatings. Additionally, spectroscopic variable angle ellipsometry, neutron reflection and atomic force microscopy carried out by Styrkas *et al.* [24] suggest that OTS films often consist of nearby multi-layer domains of different thickness, ranging from one to three and even four times the monolayer thickness.

In conclusion, the small intrinsic slip on hydrophobic homogeneous monolayer coatings provided by classical MD simulations does not emerge as an artefact of poor/improper modelling. The discrepancy with experimental values is most probably owing to imperfections in the deposited layer, a conjecture still to be confirmed by addressing realistic nanoscale defects and possible trapping of nanobubbles.

## Acknowledgements

Help with some MD simulations by D. Gentili and L. Mammucari is acknowledged. Computing resources were made available by CASPUR under HPC Grant 2010.

## Footnotes

One contribution of 25 to a Theme Issue ‘Discrete simulation of fluid dynamics: applications’.

- This journal is © 2011 The Royal Society