J. Phys. Chem. C 2007, 111, 15493-15504
15493
Anisotropic Self-Diffusion in Nanofluidic Structures† Henry Bock* School of Engineering and Physical Sciences, Heriot-Watt UniVersity, Edinburgh EH14 4AS, United Kingdom
Keith E. Gubbins Department of Chemical and Biomolecular Engineering, 911 Partners Way, North Carolina State UniVersity, Raleigh, North Carolina 27695
Martin Schoen Stranski-Laboratorium fu¨r Physikalische und Theoretische Chemie, Sekr. C 7, Fakulta¨t fu¨r Mathematik und Naturwissenschaften, Technische UniVersita¨t Berlin, Straβe des 17. Juni 135, D-10623 Berlin, Germany ReceiVed: March 7, 2007
By means of equilibrium molecular dynamics simulations, we investigate self-diffusion in a “simple’’ fluid confined to nanoscopic slit-pores. The pore walls are decorated with wettable and nonwettable chemical “stripes’’ that alternate in the x direction and are assumed infinitely long in the y direction. We consider the impact of pore width as well as variations of the width of the wettable stripes dwet. Depending on these model parameters and the thermodynamic conditions, the confined fluid may exist as one of three morphologically distinct phases: a thin fluid film adsorbed by the wettable stripes, a fluid bridge spanning the gap between the (aligned) stripes on the two opposite substrates, or a nanostructured liquid where molecules occupy the entire space between the substrate surfaces. By analyzing mean square displacements, velocity autocorrelation functions, and their power spectra, a detailed picture of mass transport and its relation to substrate decoration emerges. In particular, we find that the axial symmetry of the self-diffusion tensor D with respect to the z axis is approximately preserved in liquidlike phases, whereas no such symmetry exists as far as both film and bridge phases are concerned.
I. Introduction If a fluid is confined to tiny volumes of nanometer dimensions, the properties of such a nanoconfined fluid are altered markedly compared with the bulk phase under otherwise identical thermodynamic conditions. This is because in nanoconfined fluids the length scale on which the fluid-fluid forces vanish competes with the dimensions characterizing the volume to which fluid molecules are restricted. However, in the past, the overwhelming number of studies was concerned with equilibrium properties such as the structure of the confined fluid and its phase behavior.1-4 Comparably few studies are concerned with transport properties of confined fluids where most previous workers focus on mass transport as the perhaps simplest transport phenomenon. Diffusion of fluids in confinement is important in a variety of practical applications. For instance, Skipper et al. investigated diffusion of water and simple organic molecules in 2:1 clays by means of elastic neutron scattering to understand the swelling of these materials which is important in the context of frost heaving.5 Neuman et al. employ the so-called fluorescence recovery after photobleaching (FRAP) technique to study the lateral diffusion of surfactant molecules between two solid surfaces as a function of relative humidity.6 Another important experimental technique to study diffusion in confined geometries †
Part of the “Keith E. Gubbins Festschrift”. * To whom correspondence should be
[email protected].
addressed.
E-mail:
is nuclear magnetic resonance (NMR) spectroscopy.7-9 Stepisnik and Callaghan10 and Duh et al.11 employed NMR to measure the velocity autocorrelation function (VACF) of tagged particles in confinement, a quantity of central interest also to the present work because it permits one to get a handle on the self-diffusion of such particles from a strictly microscopic perspective. Given this status of experimental work, it is perhaps not too surprising that quite a few theoretical studies of diffusion in confined geometries have also appeared in the literature to date. For example, Almeras et al. studied the relation between diffusion and wetting phenomena at fluid-solid interfaces.12 In their study, they focus on hydrodynamic boundary conditions at the interface following in spirit an earlier study from the same group.13 In a separate study, Bocquet and Barrat used mode coupling theory to calculate self-diffusion coefficients in confinement.14 Krishnan and Ayappa used so-called Γ distributions to develop a density-of-states model to describe self-diffusion in strongly inhomogeneous confined fluids.15,16 MacElroy et al. employed kinetic theory and molecular dynamics (MD) simulations to investigate transport properties of nanoconfined fluids.17 In most cases, however, theoretical approaches are solely based on computer simulation techniques such as MD simulations. Here the focus was mostly on “simple’’ hard-sphere or Lennard-Jones-types of fluids in which each molecule possesses only translational degrees of freedom.18-22 In a few cases, molecular fluids such as nitrogen,23 carbon dioxide,24 alkane melts,25 and water26 or polar model fluids of the Stockmayer type27 have also been considered. Besides pure fluids, binary
10.1021/jp071861y CCC: $37.00 © 2007 American Chemical Society Published on Web 07/26/2007
15494 J. Phys. Chem. C, Vol. 111, No. 43, 2007 mixtures of alkanes have been studied by Galliero et al., who are not concerned with self- but thermal diffusion and the Soret effect in confinement.28 However, in almost all of these earlier studies, fluid molecules are either restricted to cylindrical or slit-nanopores. In other words, the geometry of the spaces to which fluid molecules are constrained by the external fields representing the substrates were rather simple. On the contrary, Murad et al. have systematically investigated the role of the geometry of confinement on self-diffusion.29 In their study, these authors compare results obtained for “simple’’ pore geometries like nanoslits or -cylinders with those for more sophisticated geometries such as square and cubic pores or pores where the atomic roughness of the substrates is taken into account. This latter situation has also been considered in the works by Somers and Davis30 and earlier by Schoen et al.18 However, regardless of the particular geometry of confinement, the nature of the fluid molecules, or the degree of sophistication with which the solid substrates are modeled, some general consequences of the mere presence of such substrates should also be emphasized. Just like the scalar pressure of a homogeneous, isotropic bulk-phase needs to be replaced by a pressure (or stress) tensor4 to adequately describe the thermodynamic state of a confined fluid, the diffusion constant needs to be replaced by a diffusion tensor D. As was previously emphasized by Schoen et al.18 and Liu et al.,31 D ) D1 sufficiently far away from the confining substrates where D is the self-diffusion constant and 1 is the unit tensor. Sufficiently close to a homogeneous (i.e., unstructured) solid surface, however, D is axially symmetric where the diagonal elements D| and D⊥ of the second-rank matrix representing D are associated with self-diffusion constants parallel and perpendicular to the planar substrate surfaces of a nanoscopic slitpore.18,31 Under conditions of severe confinement the axial symmetry of D may no longer be preserved as was pointed out by Su et al.32 who focus on monolayer organic films confined between atomically structured surfaces. If these monolayer films are exposed to a shear strain by misaligning the solid substrates one may tune mass transport in terms of direction and magnitude. Under these conditions it was shown by Schoen et al.33 that self-diffusion in directions parallel with the planar solid surfaces may be anomalous in the sense that the mean-square displacement (MSD) is proportional to tβ where t denotes time and β < 1 is a dimensionless number. At sufficiently long times, the MSD eventually changes from this subdiffusive to ordinary diffusive behavior characterized by a linear time dependence of the MSD (β ) 1). In this work, we will focus on anisotropic self-diffusion characterized by a diffusion tensor that no longer exhibits axial symmetry. However, unlike in the work by Su et al.,32 we consider conditions of relatively weak confinement. More specifically, we employ planar but chemically heterogeneous solid surfaces confining fluid films of several molecular layers thickness in the z direction (i.e., perpendicular to the plane of the solid substrates). The heterogeneity of the substrates is effected by wettable and nonwettable chemical “stripes’’ of a certain width. These stripes alternate in the x direction and are assumed infinitely long in the y direction. Experimentally, the decoration of solid surfaces with such nanometer substructures has become feasible more or less routinely in recent years.34-36 At the same time, quite a few theoretical studies have been carried out to elucidate the morphology and the phase behavior of fluids confined by such
Bock et al. chemically decorated substrate surfaces.37-40 Because of this background of knowledge and because of the relevance of such structured surfaces as micro- or nanofluidic devices,41 it seems timely to devote this study to an investigation of mass transport in fluids confined by chemically heterogeneous, nanostructured solid surfaces. The remainder of this paper is organized as follows. In section II, we introduce our model system (see section II.A) and give a brief account of its equilibrium properties (see section II.B) before we focus on related dynamic properties (see section II.C). Section III is devoted to a presentation of our results where we begin with a quantitative discussion of anisotropic self-diffusion (see section III.A), discuss the impact of the chemical stripes with which our surfaces are decorated (see section III.B), and address confinement effects (see section III.C). We summarize our results and put them into perspective in the concluding section IV. II. Theoretical Background A. Model System. In what follows, we consider a simple fluid confined between two plane parallel solid surfaces separated by a distance sz along the z axis of the laboratory coordinate frame. The substrates are located at z ) (sz/2. In the x and y directions, the system is assumed infinite which is effected by applying periodic boundary conditions at the faces x ) (sx/2 and y ) (sy/2 of a rectangular computational cell of dimensions sx × sy × sz, respectively. The total configurational potential energy of our system may be cast as
U ) UFF + UFS
(2.1)
where the fluid-fluid contribution is given by
UFF )
1
N
N
∑ ∑ u(rij).
(2.2)
2 i)1 j)1*i
In eq 2.2, N is the number of fluid molecules, rij ) |ri - rj| is the distance between a pair of molecules located at ri and rj, respectively, and
u(r) )
{
ushift(r), r e rc r > rc 0,
(2.3)
is the fluid-fluid interaction potential. The so-called shiftedforce potential
ushift(r) ) uLJ(r) - uLJ(rc) +
|
duLJ(r) dr
r)rc
(rc - r) (2.4)
is related to the conventional Lennard-Jones (LJ) (12,6) potential
[(σr ) - (σr ) ]
uLJ(r) ) 4
12
6
(2.5)
by truncating the latter at a cutoff radius rc ) 2.5σ and by shifting it, such that both the potential and the force between a pair of molecules vanish continuously at rc. Here is the depth of the attractive well and σ may be viewed as the “diameter’’ of a spherically symmetric fluid molecule. Since ushift(r) vanishes by definition at the cutoff radius, u(r) is explicitly short-range which is advantageous in the context of computer simulations because then thermal averages need not be corrected for longrange interactions r > rc.4 For our present model, the fluid-substrate contribution to U in eq 2.1 is given by
Anisotropic Self-Diffusion in Nanofluidic Structures 2
UFS )
J. Phys. Chem. C, Vol. 111, No. 43, 2007 15495
N
∑∑Φ[k] (Ri) k)1 i)1
(2.6)
where superscript k refers to an interaction of fluid molecule i with the lower (k ) 1) and upper (k ) 2) substrate, respectively. In eq 2.6, Ri ≡ (xi, zi) ∈ R2 is a point in the x -z plane. For the present model with chemically striped solid surfaces, we split the interaction potential into repulsive (φrep) and attractive contributions (φatt). Specifically, we take
φ[k] rep
( (
4πnAσ2 σ (z) ) 5 z ( sz/2
2 φ[k] att (z) ) 2πnAσ
) )
σ z ( sz/2
10
(2.7a)
4
(2.7b)
and and σ have the same meaning (and values) as before in eq 2.5. In eqs 2.7, nA ) 2/l 2 is the areal density of the solid surface where l /σ ) 41/3 is the lattice constant assuming the (100) structure of the face-centered cubic (fcc) lattice. In these expressions, the “+’’ sign refers to the interaction between fluid molecule i and the lower (k ) 1) surface whereas the “-’’ sign is chosen if this molecule interacts with the upper (k ) 2) substrate. To account for a decoration of the substrates with wettable stripes of width dwet, we introduce a “switching’’ function
s(x; dwet) ≡
[ (
)
(
)]
dwet dwet 1 tanh x + - tanh x 2 2 2
(2.8)
which varies continuously in the x direction as one goes from regions controlled by nonwettable parts of the substrate to those in which the substrates are wettable. Plots in Figure 1 illustrate the variation of the switching function for the three cases dwet/σ ) 2, 4, and 8 studied in this work. The function is zero over those regions that are controlled by the nonwettable portion of the substrate; it changes to a value less or equal to one in spatial regions that are dominated by the wettable stripe with which the substrates are decorated. The narrower the stripe, the narrower s (x; dwet) is and the lower the maximum it attains. This accounts in a qualitatively realistic way for the reduced amount of wettable material as the width of the stripe shrinks. Consequently, s (x; dwet) serves to “switch on/off’’ attractive fluid-substrate interactions according to [k] Φ[k] (R) ) φ[k] rep(z) - φatt (z) s (x; dwet)
(2.9)
for fluid molecules entering or leaving the region controlled by the wettable stripe. However, it should be borne in mind that in principle decoration of a solid surface with chemical material may also have an impact on the range of the fluid-substrate interaction potential. This effect is related to the thickness of the material with which the substrate is decorated and is not accounted for by our present model for simplicity. For the system characterized by the potential energy represented by eqs 2.1, 2.3, and 2.6, we solve Newton’s equation of motion numerically by means of the standard velocity Verlet algorithm.42 B. Equilibrium States and Static Properties. From eqs 2.12.9, it is clear that each fluid molecule is subject to repulsive and attractive interactions with its neighbors and with the two solid substrates. Therefore, one expects phase transitions to occur in the confined fluid under suitable thermodynamic conditions. For a fluid confined between chemically striped surfaces, it is
Figure 1. Switching function s (x; dwet) as a function of position in the x direction (sx ) 15σ) and dwet ) 2σ (-), dwet ) 4σ (- - -), and dwet ) 8σ (- ‚ -). For the system with dwet ) 8σ, ds (x; dwet)/dx is also plotted as a function of x (b).
well-known37-40 that such transitions may arise between three morphologically distinct thermodynamic phases: thin liquid films adhering to the wettable lanes (of width dwet and infinite length in the y direction) on both substrates, nanostructured liquid (i.e., a higher-density phase filling the entire space between these substrates), and liquid bridges spanning the gap between the opposite wettable substrates regions. These three fluid morphologies are illustrated by representative plots of the local density
F(R) ≡
〈
N
δ(Ri - R) ∑ i)1
〉
)
〈N(R)〉 ∆xsy∆z
(2.10)
in Figure 2 where δ(Ri - R) denotes the Dirac δ-function and N(R) is the number of molecules located in a parallelepiped of side length ∆x ) ∆z ) 0.01σ and sy centered on R where sy varies typically between 30 and 120σ depending on the thermodynamic state (i.e., the average density of the confined fluid). In eq 2.10 and elsewhere in this work, angular brackets 〈...〉 signify a configurational average in the appropriate statistical-physical ensemble under consideration. For later comparison, it is interesting to compare the above morphologies with those forming between chemically homogeneous, wettable solid surfaces. The latter may be perceived as a special case of the present system realized by setting s (x; dwet) ) 1 in eq 2.9. In this case, fluid properties become translationally invariant in the x direction as well. This is illustrated by plots of the local density for a thin film and a typical liquidlike state of a fluid confined between such chemically homogeneous substrates (see Figure 3). In this reference system, bridge morphologies can obviously not arise because of the lateral homogeneity of the substrate surfaces. In this paper, we are interested in the equilibrium dynamics of all three morphologies off any phase transitions. We therefore need to know where pairs of the three phases (or all three of them) coexist in thermodynamic state space. For any fixed temperature T, these transitions can be located most conveniently in parallel Monte Carlo simulations in the grand canonical ensemble where equilibrium states are characterized by minima of the grand potential Ω. For the present system, it follows from the detailed discussion presented in ref 4 that Ω ) τyyV where τyy is a diagonal component of the stress tensor (in Cartesian coordinates) and V is the system volume. This “mechanical’’
15496 J. Phys. Chem. C, Vol. 111, No. 43, 2007
Bock et al.
Figure 2. Density profiles F (R) of representative film (a), bridge (b), and liquid (c) phases forming in a slit-pore with chemically striped substrate surfaces where sz ) 5.5σ and dwet ) 4.0σ. For clarity we show data for two periods 2sx ) 30σ in the x direction. Note that F (R) is translationally invariant in the y direction, that is normal to the x-z plane shown in the plots.
Figure 3. Density profiles F (R) of representative film (a) and liquid (b) phases forming in a slit-pore with chemically homogeneous substrate surfaces where sz ) 5.5σ. Note that F (R) is translationally invariant in the x and y directions.
expression for Ω is a direct consequence of the translational invariance of system properties in the y direction. A molecular expression for the latter can be derived in a straightforward fashion as explained in Appendix E.3.3 of ref 4. The final result is
τyy ) -
〈N〉kBT V
+
1 4V
〈
N
N
u′(rij)rij(rˆ ij ‚ eˆ y)2 ∑ ∑ i)1 j)1*i
〉
(2.11)
where u′(r) ≡ du/dr, rij ) rijrˆ ij is the distance vector between fluid molecules i and j, and eˆ y is a unit vector pointing along
the y axis of the Cartesian coordinate system. Equation 2.11 may be employed to compute Ω in GCEMC simulations. As explained in refs 4 and 38, plots of Ω/V give rise to different branches of negative slope where the magnitude of this slope is directly determined by the mean density of the confined fluid in each of its three morphologically distinct thermodynamic phases. Hence, on account of the different slopes, intersections between pairs of branches may exist indicating phase coexistence between that pair of phases (or all three of them at a single triple point) to which the intersecting branches pertain. Once these points of phase coexistence have been located at a given T, it is a simple matter to compute in GCEMC simulations the average number of fluid molecules 〈N〉 for representative state points in the single-phase region of each morphology. Hence, for the subsequent MD simulations one may take as input a fixed number of fluid molecules N ) 〈N〉 for each morphology. We note in passing that this approach is based upon the implicit hypothesis of equivalence of statisticalphysical ensembles or, in turn, the existence of a proper thermodynamic limit for our present system.4 In the context of dynamic properties to be analyzed in subsequent sections of this paper it is important to realize that (1) F(R) is obtained as an ensemble or time average and (2) plots in Figure 2, panels a and b, indicate that the fluid density is approximately zero outside those regions where most fluid molecules are located (i.e., where F(R) has values differing substantially from zero). Therefore, the overwhelming portion of molecules located in those regions stay there on average. This notion also prompts us to conclude that there is no substantial exchange between the neighboring higher-density regions [see Figure 2, panels a and b]. Notice, however, that for the film phase [see Figure 2a] the regions of high density located at the opposing walls are connected by an essentially unstructured region of low but nonzero density in the z direction. Moreover, as was shown earlier for a highly stratified liquid-like phase, molecules stay within their original fluid layer for surprisingly long times.18 This suggests that in the x direction almost all fluid molecules are restricted to spatial regions corresponding to the maxima of F(R), whereas the low but nonzero values of F(R) along the z direction reflect a certain degree of exchange of molecules between the opposite wettable stripes. We expect these structural features to manifest themselves in specific dynamical “modes’’ that we shall analyze below in greater detail. C. Dynamic Properties. At the core of this work, however, we wish to investigate the equilibrium dynamics of the three different morphologies of a confined fluid. To this end, it will prove convenient to consider the VACF defined as
CR(t) )
1 N
〈
N
〉
VRi(t0)VRi(t0 + t) t ∑ i)1
0
(2.12)
where VRi ≡ Wi‚eˆ R is the projection of the velocity (vector) of fluid molecule i onto the R axis of the Cartesian coordinate system. Because we are dealing with systems in a state of thermodynamic equilibrium, the thermal average in eq 2.12 is performed over a sufficiently large number of statistically independent time origins t0 along the trajectory of each of these molecules. To emphasize the difference between the time average such as in eq 2.12 and the ensemble averages introduced previously in section II.B, we attach the subscript “t0’’ to the angular brackets wherever appropriate. Our results below are typically based upon 105 such time origins.
Anisotropic Self-Diffusion in Nanofluidic Structures
J. Phys. Chem. C, Vol. 111, No. 43, 2007 15497
Another quantity related to the VACF is the mean square displacement (MSD) which is obtained from the definition of the time-dependent displacement of any tagged molecule
∆ri(t) ≡ ri(t) - ri(t0) )
∞
∫t t dt′ Wi(t′)
N
∑ tf∞ i)1
〈|∆rRi(t)| 〉t0
0
2
)
2Nt
∫0∞ dt′ CR(t′) ) DRR,
R ) x, y (2.14)
sz2 6
-
16sz2 π4
∞
∑ j)1
1
π4
(2j - 1)4
[
exp -
sz2
]
(2.15)
which defines the diffusion constant Dzz associated with mass transport in the z direction despite the fact that, strictly speaking, one anticipates only nondiffusive behavior in this direction for finite sz. By nondiffusive we are referring to the fact that for finite sz it follows from eq 2.15 that
lim tf∞
1
N
∑ 〈|∆rzi(t)| 〉t 2
0
N i)1
)
sz2 6
(2.16)
that is the MSD in the z direction assumes a plateau determined by the degree of confinement (i.e., by the distance sz between the substrate surfaces). For a diffusive process, on the other hand, one can easily verify from eq 2.14 that the corresponding N
〈|∆rRi(t)|2〉t ∑ tf∞ i)1
lim
0
)∞
(2.17)
Notice that mass transport in a direction perpendicular to the plane of the solid substrates becomes diffusive as sz f ∞ (i.e.,
)
(2j - 1)4
96
(2.18a)
(2.18b)
it follows from eq 2.15 that we recover ordinary diffusive behavior in the long-time, large-pore limit. In other words
lim lim
1
N
∞ ∑ 〈|∆rzi(t)|2〉t ) ∫0 ds Cz(s) ) Dzz 0
2Nt i)1
(2.19)
for mass transport in the direction perpendicular to the plane of the two solid substrates. An in-depth investigation of the dynamic behavior of a fluid confined between the present nanostructured solid substrates will benefit greatly if the above analysis in the time domain is amended by a parallel analysis in the frequency domain. To that end, we introduce the power spectrum of the VACF, that is its Fourier transform via ∞
C ˜ R(ω) )
∫
1 dt CR(t) exp(iωt) x2π -∞
(2.20)
According to a proof presented in the book of Hansen and McDonald43 both CR(t) and C ˜ R(ω) are real and even functions of their respective arguments. Because of eqs 2.14 and 2.20, it is then easy to verify that
C ˜ R(0) )
(2j - 1)2π2Dzzt
π2
1
∑ j)1
)
0
)
8
N
∑〈|∆rzi(t)|2〉t N i)1
1
(2j - 1)2
∞
tf∞ szf∞
where for the present anisotropic fluid the diffusion constant DRR is a diagonal element of the diffusion tensor D. Equation 2.14 is essentially based upon a solution of the diffusion equation with boundary conditions appropriate for an infinite system as discussed by Hansen and McDonald in section 7.2 of their book.43 In the z direction, fluid motion in the present system is restricted on account of the solid substrates. Hence, the diffusion equation needs to be solved with different boundary conditions that take into account the fact that the diffusive flux must vanish at the plane of the solid substrates, that is at z ) (sz/2.44 Based upon this modified solution of the diffusion equation, it was shown by Schoen et al.18 that
1
∑ j)1
(2.13)
Following the derivation sketched in Hansen and McDonald’s book,43 one may derive an expression linking the MSD to a time integral over the VACF defined in eq 2.12. A discussion of some finer points, which are central to this derivation but are not normally presented in the literature, is deferred to Appendix A. If the VACF vanishes sufficiently quickly with increasing time, we may take t f ∞ in eq A5 such that43
lim
in an infinitely large slit-pore). Expanding, for example, the exponential function on the right side of eq 2.15 in a MacLaurin series in terms of 1/sz2 and using the identities
xπ2D
RR,
R ) x, y
(2.21)
As we shall see below, an analysis of power spectra will permit us to identify characteristic dynamic modes that can be linked to the dimensions of specific nanofluidic structures forming at the wettable chemical pattern with which the substrates are endowed. III. Results In what follows, quantitative results for the properties introduced in section II.C will be given in the customary “reduced’’ (i.e., dimensionless) units. In particular, time t will be expressed in units of [mσ2/]1/2, lengths will be given in units of the fluid “diameter’’ σ, and temperature in units of /kB. Throughout this work, we shall focus on a single temperature T ) 0.7 which is sufficiently subcritical with respect to the bulk gas-liquid critical point as well as with respect to both the film-bridge and bridge-liquid critical points in the confined system.38 Energies are expressed in units of where we reemphasize that the same applies to both fluid-fluid and fluid-substrate interaction for simplicity [see eqs 2.5 and 2.7]. Other reduced units are derived as suitable combinations of the ones given above. A. Anisotropic Diffusion and Fluid Morphology. Because of the difference between eqs 2.14 and 2.19 for sufficiently small sz, one immediately anticipates anisotropic diffusion as far as the direction normal (z) to and parallel (x and y) with the solid surfaces are concerned. In other words, even if the substrates are chemically homogeneous [i.e., if s (x; dwet) ) const, see eq
15498 J. Phys. Chem. C, Vol. 111, No. 43, 2007
Bock et al.
TABLE 1: Chemical Potential µ and Corresponding Mean Density Gj for Characteristic Morphologies of Confined Fluid Phases as Functions of the Width of the Chemical Stripe dwet (T ) 0.7, sz ) 5.5)a Dxxa dwet
µ
Fj
2.0 2.0 2.0 4.0 4.0 4.0 8.0 8.0 8.0
-7.87 -7.65 -7.39 -8.30 -7.90 -7.40 -8.70 -8.42 -7.40
0.021 0.171 0.435 0.023 0.178 0.483 0.030 0.185 0.531
Dyya
Dzzc
phase MSD VACF MSD VACF D h xx/D h yy MSD b
film bridge liquid film bridge liquid film bridge liquid
2.523 0.238 0.106 1.396 0.167 0.088 0.615 0.092 0.069
2.521 0.238 0.105 1.397 0.167 0.088 0.615 0.092 0.069
3.285 0.330 0.100 2.463 0.256 0.086 1.813 0.148 0.074
3.288 0.330 0.100 2.467 0.257 0.086 1.812 0.148 0.073
0.768 0.721 1.050 0.567 0.651 1.023 0.339 0.622 0.939
0.142 0.112 0.068 0.073 0.089 0.054 0.040 0.064 0.042
a Also given are diagonal components of the diffusion tensor D. From the long-time slope of the MSD and the integral over the VACF [see h RR, R ) x and y corresponds to an average of DRR from eq 2.14]. b D MSD and VACF. c From a fit of eq 2.19 to the long-time behavior of the MSD.
2.8], D should be axially symmetric with respect to the z axis, that is
(
D| 0 0 D ) 0 D| 0 0 0 D⊥
)
(3.1)
where D| ) Dxx ) Dyy, D⊥ ) Dzz. Off-diagonal elements of the matrix representing D are zero because in an equilibrium system the motion of any fluid molecule in one spatial direction is independent on average from that in any other direction. If the substrates are chemically decorated, on the other hand, the symmetry of D may be reduced reflected by Dxx * Dyy. If and under what circumstances this in-plane anisotropy in selfdiffusion arises in our present system will be investigated in this section of the paper. Before attending to that discussion, it seems worthwhile pointing out that thermodynamic states have been selected such that the mean density Fj of film, bridge, and liquid phases are roughly maintained as dwet varies between 2.0 and 8.0. The corresponding values of the chemical potential µ shift downward because the confined fluid is subject to an increasing attraction by the solid substrate as dwet becomes larger (see Figure 1 and also Table 1 for a summary of thermodynamic states). Therefore, it is not feasible to fix Fj as dwet varies. A representative plot of the MSD for a typical liquid state is presented in Figure 4a for dwet ) 4.0. As one can see from the plot, the MSDs for diffusion in the x and y directions nearly coincide and are much larger than the MSD in the z direction. Moreover, one notices from the plot that the MSDs in the x and y directions are linear functions of time for t J 10 as they should be in the case of ordinary diffusive mass transport. On the contrary, the corresponding MSD in the z direction is characteristic of nondiffusive mass transport and can be well described by the expression given in eq 2.15 where we truncate the infinite sum on the right side of the equation at jmax ) 3. By fitting that expression to the curve plotted in Figure 4a, we determine the value of Dzz given in Table 1. Despite the fact that the MSD in the z direction seems to level off over the time range considered, one realizes that its plateau value of about 2.0 is much lower than the limiting value of about 5.0 that one would have expected from eq 2.16 for sz ) 5.5. However, on account of repulsive fluid-substrate interactions, there is an excluded volume in the immediate vicinity of each substrate inaccessible to fluid molecules. Its
Figure 4. Mean square displacement 〈|∆rRi|2〉 as a function of time for R ) x (s), R ) y (- - -), and R ) z (- ‚ -); parts a-c refer to nanostructured liquid, fluid bridge, and film morphologies, respectively, where only data for R ) x and y are presented in parts b and c. Notice the difference in scale on the ordinate.
thickness in the z direction is roughly 1.0 at either substrate as one can verify from a more detailed inspection of the local density (see, for example, Figure 2). Therefore, using an effective substrate separation seff z ≡ sz - 2 in eq 2.16 a limiting 2/6 = 2.0 is estimated which agrees nearly ) value of (seff z perfectly with the plateau value that one reads off Figure 4a in the limit of long times. The same value is obtained for different values of dwet as we have verified. As expected the above anisotropy also appears in the corresponding normalized VACFs defined as [cf., eq 2.12]
C ˆ R(t) ≡
CR(t) CR(0)
where
CR(0) )
)
1
〈
〈
m
B
〉
N
VRi2(t0) t ∑ N i)1
x
〉
N
∑VRi(t0) VRi(t0 + t) t , Nk T i)1
m
R ) x, y, or z (3.2)
) 0
∞
∫ dVR VR 2πk T B
2
-∞
0
( )
exp -
mVR2
kBT , 2kBT m R ) x, y, or z (3.3) )
Anisotropic Self-Diffusion in Nanofluidic Structures
J. Phys. Chem. C, Vol. 111, No. 43, 2007 15499
assuming a Maxwellian velocity distribution for a system at thermodynamic equilibrium. One may also think of eq 3.3 as a direct consequence of the equipartition theorem which holds regardless of whether we are dealing with a bulk or a confined fluid and irrespective of the detailed nature of the confining surfaces.45,46 In Figure 5, we illustrate the anisotropic mass transport plotting C ˆ R(t) as a function of time for all three spatial directions. However, unlike in the corresponding Figure 4, we present plots only for bridge [see Figure 5a] and film phases [see Figure 5b]. Plots of the normalized VACF for the liquid phase do not exhibit any significant anisotropy in the x and y components in agreement with the corresponding MSDs plotted in Figure 4a as we have verified. The reader should realize that anisotropy is much smaller in the VACFs than it is in the MSDs. The reason is that the value of the MSD at a particular time t is related to a time integral over the corresponding VACF up to that time as eq A5 shows. Hence, the effect visible in particular in Figure 4, panels b and c, at a specific time t must be regarded as the accumulated anisotropy visible in the parallel plots of the VACFs in Figure 5, panels a and b, up to that time. However, we have verified that the effects visible in Figure 5 are significant given the statistical accuracy of our data for which we estimate typical error bars of about 10-3 from the plots of Figure 5a. Anisotropy in mass transport manifests itself also in the respective entries for self-diffusion constants in Table 1 where the deviation between Dxx and Dyy can be quite substantial. The results plotted in Figures 4 and 5 suggest that the anisotropy of in-plane self-diffusion (i.e., self-diffusion in the x and y directions) is a consequence of the chemical nanostructure of the solid surfaces. This notion is corroborated by noticing that each fluid molecule is exposed to a force in the x direction exerted on it by the solid substrate which we obtain from eq 2.9 as
F[k] x (R) ≡ )
|
ds (x; dwet) ∂Φ[k](R) ) φ[k](z) ∂x z dx
[ (
)
(
)]
dwet dwet φ[k](z) tanh2 x - tanh2 x + 2 2 2
Figure 5. Normalized VACF’s C ˆ R(t) as functions of time t for R ) x (s), R ) y (- - -), and R ) z (- ‚ -); (a) fluid bridge, (b) film phase for dwet ) 4, sz ) 5.5.
(3.4)
For lines z ) const, this force is directly proportional to the derivative of the switching function with respect to the x coordinate of a fluid molecule. A plot of ds (x; dwet)/dx in Figure 1 shows that, along any path of constant z ) c, F[k] x (R) tends to prevent fluid molecules from leaving the region controlled by the attractive chemical stripe in the (x direction. Moreover, because of F[k] x (R) there is a tendency for fluid molecules originally “outside’’ the chemical stripe to get attracted by it and pulled “inside’’. In other words, if the morphology of the confined fluid matches the structure of the chemically decorated surface, this morphology will be stabilized. However, for liquidlike phases, where the fluid morphology does not match the structure of the decorated substrate [see Figure 2c], F[k] x (R) is subdominant to fluid-fluid forces on account of the liquid’s relatively high density. This, in turn, prevents the in-plane anisotropy of self-diffusion from arising. That the in-plane anisotropy of self-diffusion is intimately linked to the chemical structure of the substrates can be verified by comparing the results discussed so far with those obtained for the reference system for which we presented plots of the local density already in Figure 3. The microscopic structure of both an adsorbed thin film [see Figure 3a] as well as that of a liquid phase [see Figure 3b] are translationally invariant in both
Figure 6. MSD as a function of time for the film morphology in the reference system with chemically homogeneous walls [see Figure 3a]. Plots in panel a for the self-diffusion in the x and y directions are indistinguishable [cf., Figure 4c]. Anisotropy in mass transport between x and y directions on one hand and the z direction is preserved as a comparison between the plots in parts a and b of the figure illustrate. Notice the difference in scale of the ordinate in parts a and b.
the x and y directions. Consequently, plots of the MSD in Figure 6 for the thin film in the reference system do not exhibit any difference between self-diffusion in the x and y directions contrary to the system with chemically striped substrate surfaces presented in Figure 4c. B. Impact of Chemical Stripes. To analyze the in-plane anisotropy in mass transport in greater detail, we reemphasize that this anisotropy arises only for film and bridge morphologies
15500 J. Phys. Chem. C, Vol. 111, No. 43, 2007
Bock et al. C ˜ y(ω), which manifests itself primarily as a difference in the self-diffusion constants [see eq 2.21]
C ˜ x(0) ∝ Dxx * Dyy ∝ C ˜ y(0)
Figure 7. Power spectra C ˜ R(ω) as functions of frequency ω for film morphologies with dwet ) 2 (s), dwet ) 4 (- - -), and dwet ) 8 (- ‚ -); (a) R ) x, (b) R ) y, and (c) R ) z. Plots are presented for sz ) 5.5.
for reasons explained above in section III.A. Generally speaking, the observed anisotropy is always smaller for bridge compared to film morphologies (see Figures 4 and 5). The reason is that unlike the film morphologies, which are restricted in the z direction to the immediate vicinity of the solid substrates, bridge morphologies extend all the way from one solid surface to the other. Thus, a substantial fraction of bridge morphologies exists in a spatial region in which F[k] x (R) is vanishingly small because φ[k](z) decays rather quickly as one moves away from the substrate surfaces [see eqs 2.7b and 3.4]. Thus, the aforementioned confinement in the x direction due to a nonvanishing F[k] x (R) affects only a certain portion of the bridge morphology, whereas the entire film morphology is subjected to such a transverse confinement. Henceforth, we shall therefore focus exclusively on mass transport in film morphologies because the anisotropy is largest for these morphologies but qualitatively representative for the denser bridge morphologies as well. To analyze characteristic features of the motion of fluid molecules, we concentrate on power spectra of the VACF defined in eq 2.20. Plots of C ˜ R(ω) in Figure 7 illustrate the anisotropy in mass transport regardless of the width of the chemical stripe. For the system with dwet ) 2 the in-plane anisotropy in diffusion turns out to be smallest [see Figure 7, panels a and b]: apart from a quantitative difference between C ˜ x(ω) and
(3.5)
the functional dependence of both power spectra on the frequency ω is rather similar in that they decay to zero monotonically. However, closer inspection of C ˜ x(ω) reveals a very broad and weak shoulder around ω ) 0.171. As we shall see shortly, this feature reflects the fact that fluid molcules are “bound’’ to the region controlled by the chemical stripe. However, as the plot of the switching function in Figure 1 reveals, the restriction on motion in the x direction is rather weak for dwet ) 2 and particles may relatively easily “escape’’ from the region governed by the chemical stripe. As Figure 1 also suggests, the impact of the chemical stripe can be expected to be more pronounced with increasing dwet. Hence, we expect that increasing the width of the chemical stripe to dwet ) 4 increases the anisotropy of the in-plane dynamics. This notion is indeed confirmed by the plot of C ˜ x(ω) in Figure 7a which no longer exhibits the monotonic decay visible in the corresponding plot for dwet ) 2. Instead a maximum is visible at a fequency ω = 0.125. These features are indicative of an increasing restriction of mass transport in the x direction. However, fluid molecules can still move freely along the infinitely long wettable lanes in the y direction. Consequently, the corresponding plot of C ˜ y(ω) in Figure 7b still decays monotonically. It is then tempting to associate the maximum in C ˜ x(ω) with a characteristic time scale of motion of fluid molecules in the x direction between the junction at which the wettable and nonwettable portions of the substrate adjoin to one another. The plot of ds (x; dwet)/dx in Figure 1 shows that sufficiently close to either solid substrate a particle starting at x ) 0 and moving in the (x direction will eventually encounter a repulsive force [i.e., an extremum in ds (x; dwet)/dx] as it approaches the junction between the wettable stripe and its nonwettable sourroundings. If this interpretation is valid we would anticipate the maximum in C ˜ x(ω) to shift downward in frequency as dwet increases. This is because the distance between the repulsive force barriers at the boundary of the wettable stripe increases the larger dwet becomes. This follows because fluid molecules possess an average kinetic energy
〈Ekin〉 )
m 2N
〈 〉 N
Wi‚Wi ∑ i)1
3 ) kB T 2
(3.6)
which we obtain directly from eq 3.3. Under the conditions of our MD simulations 〈Ekin〉 = 1.05. Hence, as dwet increases molecules travelling in the x direction between the boundaries of the chemical stripe would need more time t ∝ dwet-1 until they are eventually reflected by the repulsive barrier. Because of the inverse relation between “time’’ and “frequency’’, this longer time translates into a shorter frequency characteristic of this mode of motion. Plots in Figure 7a confirm this notion. From these plots, we notice a substantial downward shift in frequency of the maximum in the power spectrum with increasing dwet. Compared with both C ˜ x(ω) and C ˜ y(ω) the power spectrum associated with mass transport in the z direction exhibits a more complex structure as dwet increases. For example, we notice from the plot of C ˜ z(ω) in Figure 7c that this power spectrum exhibits a pronounced and nearly stationary maximum at a frequency of about ω ≈ 0.2 regardless of dwet. Because of this observation
Anisotropic Self-Diffusion in Nanofluidic Structures
J. Phys. Chem. C, Vol. 111, No. 43, 2007 15501
we may associate the appearance of the first maximum in C ˜ z(ω) at relatively low frequency (corresponding to a relatively long time scale on which this motion takes place) with a (more or less) periodic motion of fluid molecules travelling in the space between the physisorbed films on both substrates [see Figure 2a]. Moreover, as dwet increases from 2 to 8, plots in Figure 7c reveal that a second maximum appears in plots of C ˜ z(ω) at a frequency ω > 1.0. The plot in Figure 7c shows that this second maximum shifts continuously to larger frequencies as dwet decreases. This immediately suggests that the second maximum in C ˜ z(ω) represents a mode of motion caused by the chemical decoration of the solid substrates. C. Confinement Effects. From the results displayed in the preceding section, it seems reasonable to hypothesize that characteristic features visible in the plots of C ˜ x(ω) and C ˜ z(ω) reflect periodic motion of fluid molecules in spatial regions controlled by the chemical pattern on the substrate or the substrates themselves. The time scale on which this motion takes place should then be related to the dimensions of the chemical pattern or the distance sz between the substrates. To test this hypothesis, we carried out MD simulations in which we fix dwet ) 4 and vary the substrate separation over the range 5.5 e sz e 9.5. Clearly, the space accessible to a fluid molecule that is not adsorbed by the wettable stripes varies as we vary sz. However, it is important to realize that the distance these nonadsorbed fluid molecules may travel is not sz but rather
seff h < sz z ≡ 2l
Figure 8. Plots of the contour C [1] enclosing the fluid film physisorbed at the lower substrate [see eq B1] (- - -). Also shown is the line lout(x) defined in eq B6 (s) for sz ) 5.5, dwet ) 4 (a), and dwet ) 8 (b), respectively. Hatched areas represent nonwettable and the crosshatched area wettable portions of the substrate.
(3.7)
The reduction of the space accessible to any molecule travelling freely between the substrates comes about because of the increase of repulsion such a molecule “feels’’ when it approaches either substrate. The latter is caused by either the presence of the nonwettable substrate portions or by those molecules forming the film physisorbed at the wettable stripes. The reduction of space may be associated with the lines lfilm(x) and lout(x) defined in eqs B3 and B6, respectively (see Figure 8). Hence, it seems sensible to define
hl ≡
[ ( )
( )]
sx sx 1 lout - + lfilm(0) + lout 3 2 2
(3.8)
on account of the variation of lfilm(x) and lout(x) illustrated by the plots in Figure 8. Employing then this operationally sensible but albeit somewhat arbitrary definition we deduce hl ) 1.658 from data plotted in Figure 8a for sz ) 5.5 and dwet ) 4. Because the physisorbed film at either substrate is not affected by the presence of the other one, hl should remain constant as sz increases as we have verified. If the first maximum in the plot of C ˜ z(ω) reflects mass transport in the z direction of the nonphysisorbed molecules we expect ωmax ∝ 1/seff z because ω itself is inversely proportional to the time required to traverse the distance seff z . A plot in Figure 9 confirms our expectation. of ωmax versus 1/seff z Moreover, our data points in Figure 9 can be represented by a straight line through the origin. In other words, the first maximum in the plot of C ˜ z(ω) vanishes as sz f ∞ (i.e., 1/seff z f 0) which it must because confinement effects become negligible in this limit and molecules travel freely in the z direction as in a corresponding bulk system [see eq 2.19]. These results are further corroborated by considering the average thermal energy associated with a fluid molecule. From this energy one may estimate an average thermal velocity with
Figure 9. Frequency ωmax at the first maximum of C ˜ z(ω) [see Figure 7c] as a function of the effectiVe inverse pore width 1/seff z (see text) for dwet ) 4.
which a fluid molecule moves in the z direction. We define the latter as [see eq 3.3]
Vj ≡
x
kBT m
(3.9)
In the dimensionless units defined at the beginning of section III and for the temperature T ) 0.7, we obtain Vj = 0.837. This number may be used to calculate a wavelength defined as
λz ≡
Vj ωmax
(3.10)
Employing in this expression ωmax = 0.191 for the system with sz ) 5.5 and dwet ) 4 as an example [see Figure 7c], we calculate λz/2 = 2.191. This value corresponds closely to λz/2 ≈ seff z = 2.2 estimated from the plot in Figure 8a and eqs 3.7 and 3.8.
15502 J. Phys. Chem. C, Vol. 111, No. 43, 2007 The factor of 1/2 arises because a full period corresponds to the motion of a fluid molecule from one substrate to the other and back. This somewhat crude but useful analysis is further supported by an inspection of the broad second peak in C ˜ z(ω) at ωmax = 1.23 [see Figure 7c]. Using again eqs 3.9 and 3.10, we obtain λz/2 = 0.34. To rationalize this peculiar value, let us define the attractive domain caused by the wettable stripes as the spatial region enclosed by the contour C [k] defined in eq B1. Because of this definition, the size of that domain in the z direction corresponds to the spatial region in which most of the physisorbed molecules are located. Thus, the value λz/2 = 0.34 can be read off directly from the plot in Figure 8a as the width of the domain G [1] in the z direction. It then seems sensible to interpret the frequency at the second maximum in C ˜ z(ω) as a reflection of molecular motion in the attractive well caused by the wettable stripes. Since the width of the attractive domain in the z direction remains constant as sz increases, the location of the second maximum in C ˜ z(ω) remains constant as well. However, as one may infer from the plot of the switching function in Figure 1, both depth and range of the attractive well vary with dwet. Consequently, the second maximum in C ˜ z(ω) is located at different frequencies as one varies dwet as plots in Figure 7c illustrate. For example, ω = 1.430 is deduced from the plot of C ˜ z(ω) for dwet ) 8. With this value, we calculate λz/2 = 0.293 which corresponds nicely to the size of G [1] in the z direction as one can verify by inspecting Figure 8b. Further evidence for the above interpretation is provided by considering the maximum in C ˜ x(ω) located at ω = 0.125 from which we estimate λx/2 = 3.348 through an expression analogous to eq 3.10. Again an inspection of the plot in Figure 8 reveals that this corresponds to the size of the attractive domain G [1] caused by the wettable stripes in the x direction. When dwet is increased to 8, this maximum shifts to a smaller value of ω = 0.059 from which we estimate λx/2 = 7.058. Again this wavelength resembles closely the size of the domain G [1] in the x direction as the plot in Figure 8b reveals. IV. Summary and Conclusions In this work, we focus on anisotropic diffusion in confinement where the confining solid substrates are composed of two different chemical materials: a nonwettable solid decorated with wettable stripes of variable width dwet. In MD simulations of a fluid at thermodynamic equilibrium, we calculate both the MSD and the VACF for molecular motion in each of the three spatial directions. From both quantities we obtain the selfdiffusion constants DRR (i.e., diagonal elements of the diffusion tensor D). GCEMC simulations of the present system reveal the existence of three different fluid morphologies, namely a thin film physisorbed by the wettable stripes, a fluid bridge in which fluid material spans the gap between the wettable stripes located on both substrates, and a nanostructured liquid-phase filling the entire volume between the plane parallel substrate surfaces that is accessible to the fluid molecules. When MSDs and VACFs for all three morphologies are compared, it turns out that mass transport is anisotropic as far as film and bridge morphologies are concerned. On account of the mere presence of solid substrates, mass transport in the z direction (i.e., perpendicular to the plane of the solid surfaces) is nondiffusive indicated by a plateau approached by the MSD in the limit t f ∞ (sz finite). Mass transport is diffusive, on the other hand, in the x and y directions
Bock et al. where the chemical decoration of the solid surfaces causes an anisotropy reflected by, for example, the self-diffusion constants Dxx * Dyy where Dxx < Dyy. This inequality is preserved as far as film and bridge morphologies are concerned. For liquid morphologies, on the other hand, the inequality is not valid within error bars of our simulations. In this latter case, D exhibits approximately axial symmetry characteristic of cases in which the confining substrates are chemically homogeneous, that is if the external field representing those substrates depends only on the normal distance between it and the center of mass of a fluid molecule. Thus, we conclude that anistropic diffusion parallel to the solid substrates is a direct consequence of the chemical structure with which the substrates are decorated. This notion is corroborated by a parallel analysis of selfdiffusion for a film morphology in a reference system in which no such chemical substructure is imprinted on the substrates. In this latter case, the axial symmetry of D is preserved even for film morphologies. Mere confinement by the substrates, regardless of whether they are chemically decorated or homogeneous, always causes self-diffusion in the x and y directions to be distinctly different from the nondiffusive mass transport in the z direction. From a slightly different perspective, fluid molecules in film and bridge morphologies may be viewed as being exposed to two different types of confinement. One of these is due to the aforementioned presence of the solid substrates. The second type of confinement is caused by the chemical stripes that give rise to a repulsive force acting on molecules that move in the x direction from the center of the region controlled by the stripes toward the junction between wettable and nonwettable portions of the substrate. Since the stripes are infinitely long in the y direction no such force acts on molecules moving into that direction. It is this type of confinement that destroys the axial symmetry of D and causes anisotropic self-diffusion in directions parallel with the substrate plane. An analysis of the power spectra reveals a dependence of the position of their maxima ωmax on characteristic length scales present in our model. For example, we observe a dependence ˜ x(ω) is concerned. Such a of ωmax on 1/dwet as far as C dependence is to be expected if fluid molecules with a certain kinetic energy move essentially without suffering collisions with other molecules for a certain period of time. That this is so is corroborated further by the fact that in regions, where F(R) differs appreciably from zero for the film phase [see Figure 2a], only very few molecules are located. This makes plausible the assumption that for intermediate times each of these molecules experiences only a small number of collisions with its neighbors so that it moves essentially ballistically. Moreover, the observed dependence of ωmax on 1/dwet in plots of C ˜ x(ω) suggests that molecules “oscillate’’ in the potential well generated by the wettable portion of the substrate. This corroborates our expectation that molecules rarely leave the high-density regions.18 Notice, however, that as t f ∞ molecules will eventually undergo a sufficient number of collisions such that mass transport in the x direction eventually becomes diffusive in the limit of sufficiently long times. Compared with C ˜ x(ω) the structure of the corresponding power spectrum C ˜ z(ω) turns out to be more complex because it may exhibit two peaks rather than just a single one. Our analysis shows that the one appearing at smaller frequencies is independent of dwet but scales linearly with 1/sz, whereas the one at higher frequencies depends on dwet but is independent of sz. This suggests that the lower-frequency peak (corresponding to motion
Anisotropic Self-Diffusion in Nanofluidic Structures
J. Phys. Chem. C, Vol. 111, No. 43, 2007 15503
on a longer time scale) may be associated with a more or less collisionless exchange of molecules between the wettable stripes along the z direction. On the contrary, the maximum in C ˜ z(ω) showing up at larger frequencies may be interpreted as a reflection of motion within the regions of highest density at either substrate. In other words, in the z direction, molecules tend to “oscillate’’ in spatial regions controlled by the attractive well of the fluid-substrate potential for a certain period of time. However, in the z direction, the attractive well of the fluidsubstrate potential is narrower than in the x direction. This makes it easier for fluid molecules to “escape’’ the attractive well in the z direction compared with the x direction. Hence, these molecules may travel freely between the substrates once they have overcome the attractive well “binding’’ them to the wettable stripe. The repulsive barrier caused by the chemical stripes diminishes as one moves away from the solid surfaces because of the decay of φ[k] att (z) with increasing z [see eq 2.7b]. As one goes from film to bridge to liquid morphologies, a smaller and smaller fraction of the molecules is therefore exposed to the repulsive forces acting upon them at the junction between wettable and nonwettable portions of the substrate surfaces. This explains why the anisotropy of self-diffusion in the x and y directions decreases in the direction “film’’ to “liquid’’ (i.e., as the mean density of the confined fluid increases). On the basis of the assumption of ballistic motion over a sufficiently long period of time, we are also able to relate peaks of C ˜ R(ω) at characteristic frequencies ω to dimensions of nanofluidic structures. These nanofluidic structures are imprinted on the fluid by the chemically patterned solid substrates. In the present case, these patterns are represented by wettable lanes (i.e., a chemical stripe that is infinitely long in the y direction). Hence, we expect a similar mapping between characteristic frequencies of C ˜ x(ω) and nanofluidic structures if the chemical nanopattern on the substrate has a different shape such as a circle48 or an ellipse49 that have been considered in earlier studies of the phase behavior of confined fluids between such patterned substrates. Our results may also have important repercussions for the fabrication of nanodevices in which fluids need to be transported across a substrate through chemical “lanes’’ from one location to another. This is the typical situation one would have to deal with in a chemical “laboratory’’ on a chip where fluids need to be transported to a special region on the solid surface. There they may be brought to a reaction and subsequently analyzed at yet another location. From the present results, one would conclude that these chemical lanes need to be constructed in an optimal way if the fluid material is supposed to stay on them. This may be effected best if the attraction of fluid molecules by the wettable portions of the substrates is strongest without causing the fluid material to solidify. Under these conditions, the repulsive force acting at the junction between wettable and nonwettable portions of the substrate will give rise to a high anisotropy in mass transport.
authors. It made possible several trips by K.E.G and H.B. to Berlin since 2004.
Acknowledgment. One of us (M.S.) is grateful to the Department of Chemical and Biomolecular Engineering at North Carolina State University and in particular to his host and coauthor of this work, Professor Keith E. Gubbins, for his hospitality during a sabbatical stay in the fall of 2004 during which this work was initiated. This appreciation is extended to the entire Gubbins research group. Financial support from the Sonderforschungsbereich 448 “Mesoskopisch strukturierte Verbundsysteme’’ is also gratefully acknowledged by all three
Because G
Appendix A. Details of the Derivation of eqs 2.14 and 2.15. In this Appendix, we present details in the derivation of eqs 2.14 and 2.15 that are usually not presented anywhere in the literature. Following Hansen and McDonald,43 we realize from eq 2.13 that (t0 ) const) N
1
〈|∆rRi(t)|2〉t ∑ 2N i)1
) 0
∫t t dt′ ∫t t′ dt′′ CR(t′ - t′′) 0
s≡t′-t′′+t0 )
0
∫t t dt′ ∫t t′ ds CR(s - t0) 0
0
(A1) which holds for R ) x, y, or z. Defining now a function
fR(t′) ≡
∫t t′ ds CR(s - t0)
(A2)
0
we may carry out the integration over t′ by parts in eq A1 to obtain
∫t tdt′ f (t′) ) t′f (t′)| tt - ∫t tdt′ t′f′ (t′) ) tf (t) - ∫t tdt′ t′f′ (t′) 0
0
0
0
(A3)
because f(t0) vanishes by definition [see eq A2]. Employing Leibniz’s rule for the differentiation of a parameter integral,47 we obtain
f′(t′) )
d dt′
∫t t′ ds CR(s - t0) ) CR(t′ - t0)
(A4)
0
because the integrand does not depend on t′ and the lower limit of integration, t0, is a constant. Inserting now eq A4 into eq A3, changing variables according to t′ f s in the integral on the right side of eq A3, we eventually obtain
1
N
t 〈|∆rRi(t)|2〉t ) ∫t ds ∑ 2Nt i)1 0
0
( ) 1-
∫0t - t
0
s t
CR(s - t0) )
(
ds 1 -
)
s + t0 t
CR(s) (A5)
where we reemphasize that eq A5 is valid for R ) x, y, or z. To obtain the far right side of eq A5, we also employed the transformation s f s˜ ≡ s - t0 and then renamed the integration variable according to s˜ f s in the resulting expression. B. Definition of hl in eq 3.7. As a measure of extent of a physisorbed film we take the width of a domain G [k] film enclosed by the contour (see Figure 8)
C
[k] film [k] film
≡ {(x, z)|Φ[k] (R) ) 〈U〉/N}, k ) 1, 2
(B1)
has been defined such that the inequality
Φ[k] (R) e
〈U〉 , ∀R ∈ G N
[k] film,
k ) 1, 2
(B2)
is satisfied, C [k] film encloses that region in the x-z plane where most molecules of the physisorbed film are located. Moreover, it turns out that along the contour C [k] film the local density F (R) is constant. Hence ) Ffilm C
15504 J. Phys. Chem. C, Vol. 111, No. 43, 2007
{|
l film(x) ≡ z F (R) ) Ffilm C ∧
|
∂F(x, |z|) >0∧ze0 ∂|z| x
Bock et al.
}
(B3)
is that part of C [1] film (lower substrate) where a molecule travelling in the -z direction with |x| e dwet/2 would begin to experience repulsion by the molecules of the physisorbed film (see Figure 8). Outside the region of the physisorbed film, that is in regions characterized by dwet/2 e |x| e sx/2 and |z| e sz/2, F (R) is a monotonically increasing function of |z| almost everywhere (except for a vanishingly small region in the immediate vicinity of the junction between wettable and nonwettable substrate portions). Hence, solutions z˜ of the constitutive equation
∂2F (R) ∂z2
|
|x|)sx/2
! ) 0, z e 0
(B4)
correspond to two points of maximum change of F (R) at maximum distance from the wettable stripe in the x direction near the lower substrate. Taking
Fout C )
[ ( ) ( )]
sx sx 1 F - , z˜ + F , z˜ 2 2 2
(B5)
we may define a second line analogous to the one introduced in eq B3 through the expression
lout (x) ≡ {z|F (R) ) Fout C ∧ z e 0}
(B6)
along which a fluid molecule is exposed to a rapidly increasing repulsion by the nonwettable substrate portions. References and Notes (1) Gelb, L. D.; Gubbins, K. E.; Radhakrishnan, R.; SliwinskaBartkowiak, M. Rep. Prog. Phys. 1999, 62, 1573. (2) Schoen, M. Computational Methods in Surface and Colloid Science; Boro´wko, M., Ed.; Marcel Dekker: New York, 2000; p 1. (3) Alba-Simionescu, C.; Coasne, B.; Dosseh, G.; Dudziak, G.; Gubbins, K. E.; Radhakrishnan, R.; Sliwinska-Bartkowiak, M. J. Phys.: Condens. Matter 2006, 18, R15. (4) Schoen, M.; Klapp, S. H. L. ReViews in Computational Chemistry; Lipkowitz, K. B., Larter, R., Cundari, T. R., Eds.; John Wiley & Sons: New York, 2007; Vol. 24. (5) Skipper, N. T.; Lock, P. A.; Titiloye, J. O.; Swenson, J.; Mirza, Z. A.; Howells, W. S.; Fernandez-Alonso, F. Chem. Geol. 2006, 230, 182. (6) Neuman, R. D.; Park, S.; Shah, P. J. Phys. Chem. 1994, 98, 12474. (7) Shenderovich, I. G.; Buntkowsky, G.; Schreiber, A.; Gedat, E.; Shariv, S.; Albrecht, J.; Golubev, N.; Findenegg, G. H.; Limbach, H. H. J. Phys. Chem. B. 2003, 107, 11924. (8) Mitzithras, A.; Strange, J. H. Magn. Reson. Imaging 1994, 12, 261. (9) Kimmich, R.; Stapf, S.; Maklakov, A. I.; Skirda, V. D.; Khozina, E. V. Magn. Reson. Imaging 1996, 14, 793.
(10) Stepisnik, J.; Callaghan, P. T. Physica B 2000, 292, 296. (11) Duh, A.; Mohoric, A.; Stepisnik, J. J. Magn. Reson. 2001, 148, 257. (12) Almeras, Y.; Barrat, J.-L.; Bocquet, L. J. Phys. 2000, 10, 27. (13) Bocquet, L.; Barrat, J. L. J. Phys.: Condens. Matter 1996, 8, 9297. (14) Bocquet, L.; Barrat, J. L. Europhys. Lett. 1995, 31, 455. (15) Krishnan, S. H.; Ayappa, K. G. J. Phys. Chem. B 2005, 109, 23237. (16) Krishnan, S. H.; Ayappa, K. G. J. Chem. Phys. 2006, 124, 144503. (17) MacElroy, J. M. D.; Pozhar, L. A.; Suh, S. H. Colloids Surf. A 2001, 187, 493. (18) Schoen, M.; Cushman, J. H.; Diestler, D. J.; Rhykerd, C. L., Jr. J. Chem. Phys. 1988, 88, 1394. (19) Suh, S. H.; Rho, S. B.; Kim, S. C. J. Chem. Eng. Jpn. 1993, 26, 431. (20) Kim, H.; Cho, C. H.; Lee, E. K. J. Theor. Comput. Chem. 2005, 4, 305. (21) Liu, Y. C.; Wang, Q.; Liu, X. F. J. Chem. Phys. 2005, 122, 044714. (22) Mittal, J.; Errington, J. R.; Truskett, T. M. Phys. ReV. Lett. 2006, 96, 177804. (23) Ravi, P.; Murad, S. Mol. Simul. 1993, 11, 93. (24) Yang. X. N.; Zhang, C. J. Chem. Phys. Lett. 2005, 407, 427. (25) Winkler, R. G.; Schmidt, R. H.; Gerstmaier, A.; Reineker, P. J. Chem. Phys. 1996, 104, 8103. (26) Gallo, E.; Rovere, M.; Spohr, E. Phys. ReV. Lett. 2000, 85, 4317. (27) Froltsov, V. A.; Klapp, S. H. L. J. Chem. Phys. 2006, 124, 134701. (28) Galliero, G.; Colombani, J.; Bopp, P. A.; Duquay, B.; Caltagirone, J. P.; Montel, F. Physica A 2006, 361, 494. (29) Murad, S.; Ravi, P.; Powles, J. G. J. Chem. Phys. 1993, 98, 9771. (30) Somers, S. A.; Davis, H. T. J. Chem. Phys. 1992, 96, 5389. (31) Liu, P.; Harder, E.; Berne, B. J. J. Phys. Chem. B 2004, 108, 6595. (32) Su, Z.; Cushman, J. H.; Curry, J. E. J. Chem. Phys. 2003, 118, 1417. (33) Schoen, M.; Cushman, J. H.; Diestler, D. J. Mol. Phys. 1994, 81, 475. (34) O ¨ ner. D.; McCarthy, T. J. Langmuir 2000, 16, 7777. (35) Zeppenfeld, P.; Diercks, V.; David, R.; Picaud, F.; Ramseyer, C.; Girardet, C. Phys. ReV. B 2002, 66, 085414. (36) Harkema, S.; Scha¨ffer, E.; Morariu, M. D.; Steiner, U. Langmuir 2003, 19, 9714. (37) Ro¨cken, P.; Tarazona, P. J. Chem. Phys. 1996, 105, 2034. (38) Bock, H.; Schoen, M. Phys. ReV. E 1999, 59, 4122. (39) Bock, H.; Diestler, D. J.; Schoen, M. J. Phys.: Condens. Matter. 2001, 13, 4697. (40) Hemming, C. J.; Patey, G. N. J. Chem. Phys. 2004, 121, 6508. (41) Pihl, J.; Karlsson, M.; Chiu, D. T. Drug DiscoVery Today 2005, 10, 1377. (42) Frenkel, D.; Smit, B. Understanding Molecular Simulation; Academic Press: London, 2002. (43) Hansen, J. P.; McDonald, I. R. Theory of Simple Liquids, Academic Press: London, 2006. (44) Hall, P. L.; Ross, D. K. Mol. Phys. 1978, 36, 1549. (45) Valleau, J. P.; Diestler, D. J.; Cushman, J. H.; Schoen, M.; Hertzner, A. W. J. Chem. Phys. 1991, 95, 6194. (46) Diestler, D. J.; Schoen, M.; Hertzner, A. W.; Cushman, J. H. J. Chem. Phys. 1991, 95, 5432 (47) Arfken, G. Mathematical methods for physicists; Academic Press: London, 1985. (48) Sacquin, S.; Schoen, M.; Fuchs, A. H. Mol. Phys. 2002, 100, 2971. (49) Sacquin-Mora, S.; Fuchs, A. H.; Schoen, M. J. Chem. Phys. 2004, 121, 9077.
