Quasistatic Computer Simulation Study of the Shear Behavior of Bi

Nov 20, 2007 - The films become more fluid when the walls are moved out of alignment ... moving the confining substrates and analyzing the shear respo...
0 downloads 0 Views 164KB Size
Download PDF
4750

Langmuir 2008, 24, 4750-4755

Quasistatic Computer Simulation Study of the Shear Behavior of Bi- and Trilayer Water Films Confined between Model Hydrophilic Surfaces Alexander Pertsin*,†,‡ and Michael Grunze† Angewandte Physikalische Chemie, UniVersita¨t Heidelberg, Im Neuenheimer Feld 253, 69120 Heidelberg, Germany, and Institute of Organo-Element Compounds, Russian Academy of Sciences, 28 VaViloV Str., Moscow 119991, Russia ReceiVed NoVember 20, 2007. In Final Form: January 26, 2008 In this paper, our previous simulations of the shear behavior of confined water monolayers (Pertsin, A.; Grunze, M. Langmuir 2008, 24, 135) are extended to water films two and three monolayers thick. The shear response of the films is studied in the quasistatic regime corresponding to the infinitely low shear rate. In certain ranges of wall-to-wall separations, bilayer films are found to be capable of sustaining shear stress, as is characteristic of solids, while remaining fluidlike in respect of the lateral order and molecular mobility. The relation between the solidlike and fluidlike properties of the films is dependent on the relative alignment of the walls and on the period of the wall lattice. The films become more fluid when the walls are moved out of alignment and when the wall lattice is uniformly compressed or stretched with respect to the “optimum” lattice that favors crystal-like packing. Trilayer films do not sustain shear stress in the whole range of wall-to-wall separations where these films are formed.

1. Introduction Since the late 1980s, it has repeatedly been reported that liquid films confined to a thickness of a few molecular diameters show, when subjected to shear, a dramatic increase in effective viscosity, which is usually associated with “solidification” of the confined liquid.1,2 An analysis of the frequency-dependent dynamic oscillatory shear spectra of such films in terms of the effective storage and loss moduli (G′ and G′′, respectively) reveals that the shear behavior of the films is determined by the inherent relaxation time of the liquid, τ, as defined by the intersection point of the frequency dependences G′(ν) and G′′(ν).3 At frequencies ν < τ-1, the shear response is mainly viscous (G′′ > G′), whereas, at ν > τ-1, it is dominated by the elastic component (G′ > G′′). In the latter case, the liquid film may experience the stick-slip transition, provided that the shear displacement is large enough. In general, molecularly thin aqueous films follow the above outlined viscoelastic behavior, with the relaxation time τ being on the order of 0.1-1 s, that is, much longer than that in bulk water. This has been demonstrated by Zhu and Granick4 in their dynamic oscillatory shear study of aqueous electrolyte films confined between mica surfaces in a surface force apparatus (SFA). In the same study,4 the authors formulated a physically sound criterion of film solidity: G′ . G′′ at ν f 0. For the aqueous films studied by Zhu and Granick, G′ was less than G′′ at all accessible ν < τ-1; that is, the films failed to solidify. † ‡

Universita¨t Heidelberg. Russian Academy of Sciences.

(1) Yoshizawa, H.; Israelachvili, J. J. Phys. Chem. 1993, 97, 113000. Brushan, B.; Israelachvili, J. N. Nature 1995, 374, 607. Gee, M. L.; McGuinnan, P. M.; Israelachvili, J. N.; Homola, A. M. J. Chem. Phys. 1990, 93, 1895. Kumacheva, E.; Klein, J. J. Chem. Phys. 1998, 108, 7010. Raviv, U.; Laurat, P.; Klein, J. Nature 2001, 51, 413. Raviv, U.; Perkin, S.; Laurat, P.; Klein, J. Langmuir 2006, 22, 6142. Granick, S. Phys. Today 1999, 52, 26. Demirel, A. L.; Granick, S. Phys. ReV. Lett. 1996, 77, 2261. Demirel, A. L.; Granick, S. J. Chem. Phys. 2001, 115, 1498. Zhu, Y.; Granick, S. Langmuir 2003, 19, 8148. Sakuma, H.; Otsuki, K.; Kurihara, K. Phys. ReV. Lett. 2006, 96, 046104. (2) Raviv, U.; Klein, J. Science 2002, 297, 1540. Klein, J.; Raviv, U.; Perkin, S.; Kampf, N.; Chai, L.; Giasson, S. J. Phys.: Condens. Matter 2004, 16, S5437. (3) Dhinojwala, A.; Granick, J. Am. Chem. Soc. 1997, 119, 241. (4) Zhu, Y.; Granick, S. Phys. ReV. Lett. 2001, 87, 096104.

The solidity criterion suggested by Zhu and Granick4 attracts interest to the shear behavior of molecularly thin liquid films at extremely low shear rates. Unfortunately, the possibilities of dynamic shear SFA experiments with nanoconfined liquid films are rather limited, both in respect of the choice of the substrate material and in respect of the range of accessible shear rates. This imparts importance to the methods of computer simulation. A great advantage of these methods is the possibility of solving the problem in a step-by-step manner, starting from a simplified model system and then making the model more and more sophisticated by adding new features and details. In this way, the factors responsible for the solidification of confined liquid films can be successively analyzed and understood. The most popular and accessible computer simulation technique is molecular dynamics (MD). In its nonequilibrium version (NEMD), shear-force experiments are simulated by explicitly moving the confining substrates and analyzing the shear response of the confined fluid. Such an approach, however, suffers from a few shortcomings. These shortcomings have been discussed in detail by Bordarier et al.,5 and here we only note that in NEMD simulations one is forced to use unrealistically high shear rates to obtain a satisfactory signal-to-noise ratio for the quantities of interest. An alternative approach, usually referred to as quasistatic, is free from the problems involved in NEMD shear-force simulations. This approach has been suggested and developed by Schoen, Diestler, Cushman, and their coauthors in studies of the shear behavior of simple van der Waals liquids confined between various solid substrates.5,6 The quasistatic approach corresponds to the infinitely small shear rate and treats the shear deformation of a liquid film as a succession of equilibrium states corresponding (5) Bordarier, P.; Schoen, M.; Fuchs, A. H. Phys. ReV. E 1998, 57, 1621. (6) Schoen, M.; Cushman, J. H.; Diestler, D. J. Mol. Phys. 1994, 81, 475. Schoen, M.; Diestler, D. J.; Cushman, J. H. J. Chem. Phys. 1994, 100, 7707. Curry, J. E.; Zhang, F.; Cushman, J. H.; Schoen, M.; Diestler, D. J. J. Chem. Phys. 1994, 101, 10824. Curry, J. E.; Cushman, J. H. Mol. Phys. 1995, 85, 173. Schoen, M.; Hess, S.; Diestler, D. J. Phys. ReV. E 1995, 52, 2587. Bock, H.; Schoen, M. J. Phys.: Condens. Matter 2000, 12, 1545. Schoen, M.; Bock, H. J. Phys.: Condens. Matter 2000, 12, A333. Bock, H.; Schoen, M. J. Phys.: Condens. Matter 2000, 12, 1569. Bock, H.; Diestler, D. J.; Schoen, M. Phys. ReV. E 2001, 64, 046124.

10.1021/la7036313 CCC: $40.75 © 2008 American Chemical Society Published on Web 03/18/2008

Shear BehaVior of Bi- and Trilayer Water Films

to different lateral alignments of the confining substrates. Inasmuch as G′′ f 0 at ν f 0, plotting the simulated shear stress as a function of shear strain allows determination of G′ from the slope of the initial segment of the stress-strain curve. It should be noted that the quasistatic approach is not a universal tool for exploring the shear behavior of confined liquid films. This approach, where the viscous component of the shear stress is disregarded, is only useful in the limiting case of extremely low shear rates, which are of interest in the present work. For aqueous films, the lowest values of ν accessible in dynamic oscillatory shear experiments are more than 1 order of magnitude less than the respective τ-1.4 In this frequency range, the shear process is expected to be slow enough for the film to relax to a structure close to the equilibrium one at each particular lateral alignment. So, quasistatic simulation results should in principle be comparable with dynamic shear experiments at the lowest frequencies. For further details of the quasistatic approach in comparison with NEMD, the reader is referred to the above cited paper by Bordarier et al.5 In our recent paper,7 the quasistatic approach was applied to monolayer water films confined between hydrophilic model walls. The equilibrium states of the system at different lateral alignments were simulated by the grand canonical Monte Carlo (GCMC) technique. Contrary to the data reported by Zhu and Granick,4 the simulation resulted in a nonzero G′, so that the solidity condition was fulfilled. In our view, the disagreement between the simulation7 and experiment4 originated from the difference between the experimental and simulated systems. On one hand, the dynamic shear measurements4 were concerned with aqueous electrolyte solutions confined between mica sheets. The latter, when in contact with water, are known to acquire a negative surface charge because the surface K+ ions of mica desorb and dissolve in water.8 At intersheet separations of one to two water diameters, the neutralization of this surface charge requires an increase in the local concentration of counterions between the mica sheets to as high as 15 M (about three water molecules per ion).4 On the other hand, our simulated model system represented pure ionless water confined between electrically neutral hydrophilic walls. Such a system was deliberately chosen for simulations to see whether the experimentally observed fluidity of aqueous electrolyte solutions between mica surfaces may be associated with the presence of ions, as suggested by Raviv, Klein, and co-workers.2 The shear stress-strain curves simulated for water monolayers each showed an extremum that defined the yield point and yield stress, Sc.7 In an imaginary shear force experiment, the simulated system would thus be capable of experiencing the stick-slip transition: With increasing applied shear stress, sliding of one confining wall relative to the other would not be possible until the applied stress exceeded Sc. In our simulations,7 the highest solidity (i.e., the highest Sc and G′) was observed when the monolayer film formed a two-dimensional lattice commensurate with that of the confining walls. However, the film solidity still persisted when the wall’s lattice period was altered, the commensurability was lost, and the films became fluidlike, as judged from their pair distribution functions and the mean-square displacements (MSDs) of their constituent molecules. Beyond the yield point, the initially crystal-like monolayer film melted but its structure and molecular motion differed in many respects from those typical of normal fluids. In particular, the dependence of MSD on the GCMC chain length showed chaotic jumps, which (7) Pertsin, A.; Grunze, M. Langmuir 2008, 24, 135. (8) McGuiggan, P. M.; Israelachvili, J. N. J. Mater. Res. 1990, 5, 2232.

Langmuir, Vol. 24, No. 9, 2008 4751

Figure 1. Top view of a fragment of the confining walls. The unit cell of the wall lattice is shown by the hatching. The open and filled circles are in the lower (z ) 0) and upper (z ) h) walls, respectively. The designations are explained in the text.

were presumably associated with collective migrations of water molecules along the shear direction. In this paper, our previous simulations7 of confined water monolayers are extended to a range of larger wall-to-wall separations, capable of accommodating bi- and trilayer water films. 2. Method The simulation model and method have been described in detail in our previous article.7 In brief, the configuration of the computer experiment was similar to that of SFA measurements: The water film was confined in a slit between two parallel walls and was allowed to exchange water molecules with a (fictitious) bulk water reservoir at fixed temperature T, volume V, and chemical potential µ. Fluctuations of the number of confined water molecules, N, were simulated by means of the GCMC technique. To improve the acceptance of insertion attempts, which is usually too low in liquids, we used the excluded volume mapping,9 the Swendsen-Wang filtering,10 and a rotational bias procedure.11 The constraining walls were each represented by an array of force sites (or “atoms”) arranged on a two-dimensional hexagonal lattice of period l. The unit cell of the lattice was taken to be rectangular, with periods a ) x3l, b ) l (Figure 1). Most simulations were carried out with l ) 3 Å. This value corresponds to the formation of crystal-like water monolayers, which showed the highest Sc in our previous study.7 Some comparative simulations were also performed with l ) 2.76 and 3.2 Å. The lateral dimensions of the simulation box were 7a × 11b (36.4 × 33 Å at l ) 3 Å), whereas the box height was determined by the wall-to-wall separation (i.e., the slit width), h, which was incrementally varied in the range between 6.5 and 16 Å. The increase of the lateral dimensions to 11a × 18b (57.2 × 54 Å at l ) 3 Å) affected the simulation results only slightly. This suggests that the simulation box used was large enough to neglect the effects of artificial periodicity imposed on the confined water film by the periodic boundary conditions. The shear strain was applied to the upper wall in the direction of the y-axis, so that the positions of the respective force sites in the upper and lower walls were related by the following equations: x′ ) x, y′ ) y + βb, and z′ ) z + h, where the prime refers to the upper wall (Figure 1). Note that at β ) (k the walls are in registry, while at β ) (k + 0.5 they are out of registry, with k being an integer. The shear stress Tzy was defined as the force per unit area that must be applied in the y-direction to keep the upper wall fixed at a specified registry.12 With this definition, the film resists (assists) shear if Tzy > 0 (Tzy < 0). Considering that the shear stress in the y-direction, (9) Stapleton, M. R.; Panagiotopoulos, A. J. Chem. Phys. 1990, 92, 1285. (10) Swendsen, R. H.; Wang, J.-S. Phys. ReV. Lett. 1987, 58, 86. (11) Shelley, J. C.; Patey, G. N. J. Chem. Phys. 1995, 102, 7656. (12) In our previous article,7 Tzy was defined with the opposite sign in terms of the force exerted by the film on the upper wall.

4752 Langmuir, Vol. 24, No. 9, 2008

Pertsin and Grunze

Tzy(β), is antisymmetric with respect to β ) 0.5, that is, Tzy(β) ) -Tzy(1 - β), we restricted an analysis of the shear stress-strain relations by the symmetrically independent range of registries, 0 e β e 0.5. This range was explored by a series of independent GCMC runs at some selected values of β within the range (usually at 0.1 intervals). The lengths of the equilibration and production phases of the simulation at a given β were equal to 106 GCMC cycles, with each comprising N0 moves, where N0 is the initial number of water molecules in the given cycle. In doubtful cases, two or three independent GCMC runs using different starting configurations and different sequences of random numbers were carried out. The set of calculated quantities included the normal pressure, PN, the y-component of the shear stress, Tzy, the water density profile, F(z), and also the in-plane pair distribution function g(r; z). In calculations of the two latter quantities, the simulation box was divided into slices 0.1 thick, lying parallel to the walls. During the GCMC run, F(z) and g(r; z) were averaged within each individual slice and then referred to the z-coordinate of the slice center. The mobility of water molecules was assessed in terms of MSD, d, and also its lateral and normal components, dxy and dz, respectively, as defined by the equations

Figure 2. Dependence of normal pressure PN on the slit width h at l ) 3 Å and β ) 0. (The meaning of points A-D is explained in the text).

d ) dx + dy + dz dxy ) dx + dy dχ ) dχ(n) )

1

(1)

N

∑ (χ N

i 0

- χin)2,

(χ ) x, y, z)

i)1

where χi0 and χin are, respectively, the χth coordinates of the oxygen atom of molecule i in the beginning and the end of the Monte Carlo run of length n. The interaction of water molecules between themselves was described with the four-site TIP4P model,13 whereas the wallwater interaction potential was designed so as to mimic the interaction of the water molecule with a hydrophilic proton-acceptor surface. Each wall atom A interacted with water’s O and H atoms through inverse power (m - n) potentials: φ(r) ) [n(r0/r)n - m(r0/r)m]/(m - n),

(m < n)

Figure 3. Water density profiles, F(z), at l ) 3 Å and β ) 0 for slit widths h indicated by the arrows in Figure 2.

(2)

where  and r0 are the depth and position, respectively, of the potential well. The potential parameters in eq 2 were calibrated based on the following conditions: (i) The binding energy of a water molecule to the wall with a given lattice period l is the same as that for the lowest energy TIP4P water dimer (-6.24 kcal mol-1); that is, a water molecule sees the wall as favorable for binding as another water molecule. (ii) The equilibrium separation of the water oxygen from the wall is in the range 3.2-3.4 Å. (iii) The most favorable configuration of the A‚‚‚O-H bond is close to linear, typical of hydrogen bonds. The values of potential parameters satisfying these conditions are presented in the preceding paper.7 Throughout the simulations, the direct interaction between the walls was ignored. This was justified by the fact that the structure and shear behavior of confined water at a fixed h are independent of the wall-wall interaction potential.

3. Results and Discussion Normal Pressure and Water Density Profile. We begin with the behavior of PN and its relation to changes in the water density profile F(z) on variation of h at β ) 0 (the walls in registry). In our previous simulations,7 the focus was on water monolayers formed at compressive PN observed at h e 6.4 Å (Figure 2). Now, we turn to larger values of h. As h is increased beyond 6.4 Å, PN becomes negative (stretching). The initially “monolayer” form of F(z) remains unimodal up to the minimum of PN at h (13) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926.

Figure 4. Typical configuration of the system at h ) 8.7 Å, l ) 3 Å, and β ) 0 (point B in Figure 2).

) 7 Å, and then, at h ) 7.2 Å, it changes to bimodal (point A in Figure 2 and the relevant F(z) in Figure 3). The next important change occurs at h ) 8.7 Å (point B in Figure 2), just after the maximum of PN is reached: The density profiles of the individual monolayers in the bilayer film become “isolated” in the sense that the water density at the midplane (z ) h/2) vanishes (Figure 3). A snapshot of the system at point B is shown in Figure 4. The monolayers remain isolated up to the minimum of PN at h ) 9.6 Å. At h g 9.8 Å (point C in Figure 2 and the corresponding profile in Figure 3), the profiles are characterized by a nonvanishing density at the midplane. The fourth characteristic point of the system is at h ) 10.6 Å, where a third water density maximum at the midplane appears, indicative of the formation of a trilayer film (point D in Figure 2 and the corresponding curve in Figure 3). At all larger values of h tried, up to 16 Å, PN remained compressive, while F(z) did not vanish at all z, except for the excluded volume regions near the walls. As with water monolayers,7 the magnitude of PN in bilayer films is dependent on registry β. This is illustrated in Figure 5

Shear BehaVior of Bi- and Trilayer Water Films

Figure 5. Dependence of normal pressure PN on registry β at h ) 8.9 Å and l ) 3 Å.

for the film formed at h ) 8.9 Å. Unlike the monolayer films, however, the bilayers show a decrease in PN with increasing β. With increasing h, the dependence of PN on β decayed and practically disappeared when PN changed its sign at h ) 9.2 Å. The implication of the effect of registry on PN is that, in SFA shear-force experiments under a fixed normal pressure, the slit width h may vary with shear strain. The magnitude of this variation, ∆h, can be roughly estimated as ∆h ) ∆PN(∂PN/ ∂h)-1, where ∆PN is the observed variation in PN. In the case illustrated in Figure 5, ∆PN ) 1.8 kbar, the derivative (∂PN/∂h) evaluated from Figure 2 as the slope of PN(h) at h ) 8.9 Å is about 60 kbar Å-1, so that the expected ∆h is as small as 0.03 Å. More significant variations in h can be expected at separations close to the maximum of PN(h) at h ) 8.6 Å. Regions of Nonzero Shear Stress. Before proceeding to a detailed analysis of the shear stress-strain relations in confined water films, we studied the behavior of Tzy as a function of h at a fixed registry, β ) 0.25. The aim was to locate the regions of slit widths where Tzy is distinct from zero. The choice β ) 0.25 was based on our previous experience with water monolayers,7 which exhibited the yield point, βc, close to the center of the symmetrically independent region of registry values (0 e β e 0.5). Our subsequent simulations of the dependence Tzy(β), which will be discussed below, have shown that the values of βc in the thicker water films are also close to 0.25, so that Tzy(0.25) can serve as a good estimate for the yield stress Sc. The dependence of Tzy(0.25) on the wall-to-wall separation is presented in Figure 6. The figure does not include the simulation results for trilayer films at h g 10.6 Å, where the magnitude of Tzy(0.25) did not exceed the statistical uncertainty. The behavior of Tzy(0.25) for bilayer films can be qualitatively rationalized in terms of the ability of their constituent monolayers to maintain lateral coupling between themselves. On going from point A to point B (see Figures 2 and 3), the separation between the two maxima of F(z) (hereafter lm(z)) steadily increases. Nevertheless, it remains too small for the formation of bridge hydrogen bonds between the monolayers corresponding to these maxima. (At h ) 8 Å, for example, lm ) 1.8 Å.) At the same time, the water density between the maxima steadily decreases with increasing h. As a consequence, the lateral coupling between the monolayers weakens and Tzy(0.25) decreases. In the range 8.2 Å e h e 8.6 Å, Tzy(0.25) vanishes and the bilayer films lose their solidity. At point B, lm reaches 2.5 Å, fairly close to the O‚‚‚O distance in the linear TIP4P water dimer (lOO ) 2.75 Å).13 This results in the appearance of a noticeable shear stress. The largest Tzy(0.25) occurs at h ) 9.1 Å. In this case, lm ) 2.74 Å, which is practically

Langmuir, Vol. 24, No. 9, 2008 4753

Figure 6. Dependence of shear stress Tzy on the slit width h at β ) 0.25 and l ) 3 Å.

Figure 7. Shear stress Tzy as a function of registry β for different h at l ) 3 Å.

coincident with lOO. At h g 9.2 Å, lm becomes larger than lOO and the difference between these quantities increases with increasing h. As a consequence, Tzy(0.25) drops and eventually vanishes at h ) 9.8 Å, when lm ) 3.2 Å. Shear Stress-Strain Relations. Now, we turn to the behavior of Tzy as a function of β, as depicted in Figure 7. The curves in Figure 7 refer to bilayer films that are formed in the range 8.7 Å e h e 9.2 Å and show the highest Sc. (At larger h, the bilayer films with a nonzero Sc had a qualitatively similar shear response.) The solidity of the films is evidenced by a nonzero slope of the initial segments of the curves. The magnitudes of Sc are given by the heights of the maxima of Tzy(β). The yield points are either at β ) 0.2 or at β ) 0.3, which justifies our preliminary estimates of the shear resistance of the films in terms of Tzy(0.25). It can be seen that the symmetrically independent range of β includes both shear resisting and shear assisting regions. A similar situation was earlier observed by Curry14 in her GCMC simulation of a system of spherical particles modeling an octamethylcyclotetrasiloxane monolayer confined between mica surfaces. By contrast, water monolayer films simulated in our previous study7 resisted shear at all β between 0 and 0.5. The films assisted shear only in the second half of period b (0.5 e β e 1), where Tzy changed its sign by symmetry. (14) Curry, J. E. J. Chem. Phys. 2000, 113, 2400. Note that the author interpreted the transition from the shear resisting to shear assisting ranges as an atomic-scale stick-slip transition. Such an interpretation implied that the upper wall would slip during the shear assisting range, if it would be allowed to move.

4754 Langmuir, Vol. 24, No. 9, 2008

Figure 8. In-plane pair distribution functions g(r; z) at h ) 9.1 Å and l ) 3 Å: (a) at β ) 0 and z ) h/2 ( 1.8 Å (close to the edges of the water density profile F(z)); (b) at β ) 0 and z ) h/2 ( 1.4 Å (at the maxima of F(z)); and (c) at β ) 0.4 and z ) h/2 ( 1.4 Å.

Pair Distribution Functions and MSDs. As with water monolayers,7 the wall lattice with a period l ) 3 Å imprinted a long-range lateral order on confined water bilayers. This order can be appreciated from Figure 8 that shows three in-plane pair distribution functions g(r; z) at h ) 9.1 Å. Curves (a) and (b) both refer to walls in registry (β ) 0) but differ in the position of their associated planes with respect to the characteristic points of the water density profile F(z). Curve (a) describes g(r; z) at z ) h/2 ( 1.8 Å, which is close to the edges of F(z) (i.e., near the very interface with the wall), whereas curve (b) refers to the maxima of F(z) observed at z ) h/2 ( 1.4 Å. (Note that the average water density in the former plane is 2 orders of magnitude lower than that in the latter.) Finally, curve (c) also refers to the maxima of F(z) but at β ) 0.4. As seen from Figure 8, the most pronounced lateral order is characteristic of the “interfacial” water (curve (a)), whose g(r; z) is practically identical to that of crystal-like water monolayers (see Figure 5 in ref 7). On the other hand, the overwhelming majority of water molecules, as represented by curve (b), show substantially weaker lateral order. Up to βc ) 0.3, the shear had practically no effect on g(r; z). However, a further increase of β to 0.4 (curve (c)) led to weakening of the lateral order. Although the bilayer water films formed in the separation range 8.7 Å e h e 9.6 Å are solid from the viewpoint of the shear stress-strain relations, they behave as fluids in respect of the lateral mobility of their constituent molecules. This can be appreciated from Figure 9 that shows the dependence of the lateral component of MSD, dxy, on the number of GCMC cycles, n, at h ) 9.1 Å. Unlike the crystal-like monolayer films formed at l ) 3 Å,7 which showed up to βc only slight oscillations of dxy(n) around a constant value, the behavior of dxy(n) in the shearresistant bilayer films is well described by a linear dependence typical of classical fluids. In a comparison of the dxy(n) curves in Figure 9, the slope of the curves, D, can be regarded as a measure of molecular mobility, similar to the diffusion coefficient in Einstein’s diffusion equation. The values of D are listed in the legend to Figure 9. It can be seen that the increase of β from the in-registry point (β ) 0) to the yield point (βc ) 0.3) leads to a comparatively small increase in D (∼30%). Beyond βc, however, D shows a substantial rise, such that the overall increase in D is about 3-fold. This result together with the above-discussed behavior of g(r; z) provide evidence that the confined water film becomes more fluid when the shear strain exceeds βc.

Pertsin and Grunze

Figure 9. Lateral mean-square displacement of water molecules as a function of the length of the simulation run, n, for different β at h ) 9.1 Å and l ) 3 Å.

Figure 10. In-plane pair distribution functions g(r; z) at the maxima of F(z) at h ) 9.2 Å and β ) 0 for compressed (l ) 2.76 Å) and stretched (l ) 3.2 Å) wall lattices. Table 1. “Diffusion Coefficient” D of Water Molecules Confined between Structureless and Structured Walls at h ) 9.2 Å, in Comparison with That for Bulk Water (10-5 Å2) structured walls bulk water

structureless wall

l ) 3.0 Å

l ) 2.76 Å

l ) 3.2 Å

90

22

0.8

6.1

2.7

Stretching and Compression of the Wall Lattice. In the case of water monolayer films,7 the uniform stretching or compression of the wall lattice by changing the lattice period l deteriorated the commensurability of the wall and monolayer lattices observed at l ) 3 Å. As a consequence, the lateral wallwater coupling weakened and Sc decreased. Similar trends are observed with confined bilayer films. The weakening of the lateral coupling manifests itself in the behavior of g(r; z), which shows substantial damping of oscillations with r (cf. Figure 10 with Figure 8, curve (b)). Another manifestation is a dramatic increase in the “diffusion coefficient”, D, which can be appreciated from Table 1. Aside from the simulated values of D at l ) 3.0, 2.76, and 3.2 Å, Table 1 also presents, for comparison, the respective results for bulk water and for water confined between structureless walls. In the latter case, the wall-water interaction was described by an inverse power (9-3) potential dependent only on the distance of the water oxygen from the wall. The well depth and position of the (9-3) potential were taken to be the same as for the structured

Shear BehaVior of Bi- and Trilayer Water Films

Figure 11. Shear stress Tzy as a function of registry β at h ) 9.2 Å for compressed (l ) 2.76 Å) and stretched (l ) 3.2 Å) wall lattices.

walls. The substantial decrease in D observed on going from bulk water to the structureless wall demonstrates the effect of reduced dimensionality due to confinement of the system in the z-direction, whereas the further decrease occurring on going from the structureless wall to the structured ones gives an idea of the contribution of lateral wall-water coupling to the mobility of water molecules. The shear stress-strain relations in compressed and stretched bilayer films are depicted in Figure 11. Comparing these data with the respective curve at l ) 3 Å (Figure 7), one can see that the uniform stretching or compression of the wall lattice by only 7-8% leads to an approximately 3-fold drop in the yield stress. Nevertheless, the systems with compressed and stretched wall lattices remain solid from the viewpoint of their shear behavior while being fluid in respect of their pair distribution functions and MSDs.

4. Conclusions In this article, we have extended our previous GCMC simulations7 of water monolayers confined between electrically neutral hydrophilic walls to water films two and three monolayers

Langmuir, Vol. 24, No. 9, 2008 4755

thick. In common with confined monolayers,7 bilayer films are found to be capable of solidifying when sheared in the quasistatic regime. Noticeable solidity of water bilayers (in the sense of the magnitudes of Sc and G′) is, however, observed in a rather narrow interval of wall-to-wall separations, only 0.9 Å in width (8.7 Å e h e 9.6 Å). It is in this interval where the separation between their constituent monolayers is favorable for the formation of bridge hydrogen bonds responsible for the lateral coupling between the monolayers. Small, albeit statistically significant values of Sc and G′ are also observed in the “transitional” range, 6.4 Å e h e 8.0 Å, where the water density profile F(z) is steadily transformed from the monolayer to bilayer form. An important feature of confined bilayer films is that their solidlike shear behavior is not associated with film crystallization. Whereas the films can sustain shear stress, the MSDs of their molecules are described by a linear function of the simulation chain length n, which is characteristic of fluids. Similar duality in the film properties was also observed in confined monolayers,7 but the latter were, in addition, capable of forming crystal-like commensurate systems with MSDs independent of n. In common with water monolayers, water bilayer films show a decrease in the lateral order and an increase in the lateral molecular mobility (i.e., become more fluid) when the wall lattice is uniformly compressed or stretched with respect to the “optimum lattice” (l ) 3 Å) that favors crystal-like packing. Similar changes in the structure and mobility occur in bilayer films when β is increased beyond βc. In this range, the behavior of bilayer films differ in two respects from that of monolayers. First, at β ) 0.4, water bilayers change from shear resisting to shear assisting, whereas water monolayers resist shear in the whole range of β between 0 and 0.5. Second, the curves dxy(n) in confined bilayer films show no giant oscillations and hops that were observed in confined monolayers7 and were presumably associated with collective migrations of molecules in the direction of shear. A likely reason is a much lower probability of such collective motions in the case of bilayer films. Acknowledgment. This research was supported by the Office for Naval Research and the Deutsche Forschungsgemeinschaft. We are grateful to Prof. S. Granick for helpful discussion. LA7036313