- Email: [email protected]
Molecular dynamics study of self-diffusion in a liquid-liquid interface
Volume 147, number 2,3
CHEMICAL PHYSICS LETTERS
3 June 1988
MOLECULAR DYNAMICS STUDY OF SELF-DIFFUSION IN A LIQUID-LIQUID INTERFACE Marc HAYOUN CEN Saclay,D. Tech, 91191 Gif-w-YvetteCedex, France
Madeleine MEYER LPM-CNRS,I placeA. Briand. 92195 Meudon Cedex. France
and Pierre TURQ LaboratoireElectrochimie,UniversitP Pierreet Marie Curie, 8 rue Cuvier, 75005 Paris,France Received 11 March 1988; in final form 29 March 1988
The self-diffusion coefficients have been calculated in a model of a planar interface between two almost immiscible liquids. The simulations have been performed by the molecular dynamics method, using Lennard-Jones potentials. The diffusion process is isotropic in the bulk regions of the simulated system, whereas it exhibits anisotropic behaviour in the interfacial region: the transverse diffusion coefficient (parallel to the interface) is higher than the normal one. A qualitative explanation of this behaviour is suggested by comparison with the pressure tensor.
1. Introduction Most results obtained on the liquid-liquid interface (LLI) deal with static properties investigated from the experimental [ l-31 or from the theoretical [ 1,4-6 ] point of view. Simulation techniques, which have been fruitful for the study of the liquid-vapor interface [ 7,8 1, have not been applied to LLI, with the exception of a Monte Carlo simulation of the water-benzene LLI [ 91 and a molecular dynamics modelling of a LLI [ 10 1. In a previous paper [ 111, we have been able to model an interface between two immiscible liquids, using the molecular dynamics (MD) simulation technique and Lennard-Jones (LJ) pair interactions, In that work, we have shown that the average position and thickness ( = 3 atomic diameters) of this .LLI were stable over the time scale of the simulation. To characterize this interface we have computed the following properties: density and pressure profiles and interfacial tension (20 & 5 erg/cm’). The anal-
ysis of these profiles shows the existence of bulk phases extending over half of the overall volume of the model system. In this paper, using the same model [ 111, we focus attention on the mass transport properties. We have computed the self-diffusion coefficients as a function of the distance to the interface. The results exhibit a diffusional anisotropy in the interfacial region.
2. Model and computations The detailed characteristics of the model have been described [ 111, and the main features are recalled hereafter. The MD simulations have been performed in the (N, V, E) ensemble and the equations of motion were integrated using Verlet’s algorithm [ 121 with a time step 6t = 10 -I4 s. The simulated system contains 864 particles of liquid 1 (Ll ) and 864 particles of liquid 2 (L2). The computational cell is parallelipiped with edges L, L, and i, equal
0 009-2614/88/$03.50 0 Elsevier Science Publishers B.V. ( North-Holland Physics Publishing Division )
203
CHEMICAL PHYSICS LETTERS
Volume 147, number 2,3
3 June 1988
Table 2 MD results for liquid 1 (Ll) and liquid 2 (L2) equilibrated separately ‘) T (K)
P
( 10zocm-))
170 197
133.5 134.2
0.70 0.71
P Ll L2
(kbar)
a1The values for p and T have been averaged over 4000 time steps.
t
a Fig. 1. Schematic drawing of the MD simulation cell. Two adjacent boxes containing respectively Ll and L2 are presented in their initial configuration. The two interfaces a and b are indicated by the arrows; b has been generated by the periodic boundary conditions along Z. LX=&= 37 A, Lz= 69 A. The originof the referenceframeis at the centreofthe computational cell.
Table 1 Lennard-Jones parameters for the system used to simulate the LLI t,tlkr,” (K)
eal/G (K)
e,zlk, (K)
a (A)
100
200
75
3.591
with a IJ pair potential function already used in a study of unstable mixtures [ 131:
where LY- 1, 2 and fi= 1, 2 refer respectively to the atoms of L 1 and L2 (see table 1) . The potential cutoff radius r, is equal to 2.50 and the atomic masses are the same for the two liquids (m = 1.02747 x 1O-22 g). We have shown that this pair potential yields a stable interface between almost immiscible liquids
[111*
a) k, is the Boltzmann constant.
respectively to 37, 37 and 69 8, (fig. 1). In a first stage, we equilibrate separately the two liquids Ll and L2 at the same temperature and the same pressure. The two subsystems are then joined to build the computational cell. The initial position of the planar interface thus obtained is parallel to the XY plane as indicated in fig. 1. Periodic boundary conditions have been chosen. They generate a second interface in the 2 direction (fig. 1). The interactions between particles are calculated
The properties of the two subsystems (Ll, L2) equilibrated separately, are given in table 2. Table 3 contains results of the LLI simulation, including the estimated values of the solubilities between these two liquids.
3. Results and discussion The mass transport properties in the LLI have been investigated through the MD computation of the selfdiffusion coefficients. These quantities have been obtained from the usual least-squares fits of the meansquare displacements (MSD) [ 141. The selfdiffu-
Table 3 MD results for the LLI a) PI
P2
L2inLl (%)
T (IQ
P b’
( 1020cm-3)
Ll in L2 (%)
P
(1020cm-3)
(kbar)
(kbar)
17522
194f2
EO.3
zo.1
135.2
0.72
1.08
a) All these quantities have been averaged over 20000 time steps. p, nd pzreported here correspond to the plateau values of the density profiles. b, Without long-range correction.
204
3 June 1988
CHEMICAL PHYSICS LETTERS
Volume 147, number 2,3
Th
I
12
I1
a
b Fig. 2. Projection in the XZ plane (see fig. 1) of the regions used for the computation of the self-diffusion coefficients. 11:interfacial region in Ll. 12: interfacial region in L2. Bl: bulk region in Ll. B2: bulk region in L2. The four regions have the same total thickness AZ=%0 A.
sion coefficient along the X axis for instance, Dx, is calculated with the following relations: (sX2(r))=~~,~~~~,[Xi(7j+[)-X;(q)]2, T
(2)
(8X2(t))
(3)
=2D,t+const.,
where N is the number of particles taken into account and N, is the number of time origins. In all the computations performed, N, has been chosen equal to 250 and the maximum value of the time t is 5006t. The MSD, along the X, Y and Z axes, have been calculated in four different parts of the system, in order to obtain the diffusion coefficients characteristic of the interface and bulk regions. These four parts consist in slabs parallel to the average orientation of the interface. I1 and 12 are located on each side of the interface respectively in Ll and L2, whereas Bl and B2 correspond to the bulk of each liquid. Since the computational cell contains two stable LLI [ 111, the regions 11 and 12 are each composed of two equivalent slabs of thickness 4AZ, as indicated in fig. 2. The slabs Bl and B2, each of thickness AZ, are located in the middle of the bulk phases ( see fig. 2 ) . For 11, 12, B 1, and B2 the particle numbers were variable. Indeed, in these regions, the particles have been selected every lOOO&,over a run of 100006~. On average, the number of particles of type 1 used
,
0.0 -34.5
I
I
34.5 0.0 z (%I) Fig. 3. Upper part: percentage of L2 as a function of Z, averaged over 2000081 [ 111. Lower part: normal and transverse pressure profiles PN(Z) and PT( Z) (computed without long-range correction). At the second interface (b, see fig. 1) these profiles are extrapolated for Zz34.5 A (open symbols), using the image provided by the periodic boundary conditions. The verical dashed lines refer to the average positions of the interfaces. The contributions of the two liquids to the interfacial tension y are pointed out: Ll (dark areas), L2 (hatched areas). (ya
IV’,(Z)-J’,(Z)ldZ). is 74 for 11 and 119 for Bl, while for 12 and B2 the number of particles of type 2 are respectively 102 and 136. fAZ~2.5 8, has been chosen smaller than the interfacial thickness (Th= 10 A), deduced from the pressure and density profiles [ 111: (i) to take the significant part $Th of the interfacial region (e.g. for L2, between Z1 and Z2 in fig. 3), (ii) 4AZ must be smaller than jTh, the difference between JAZ and f Th being approximately given by (DlOOO&) ‘12, where D is a reasonable estimation of the diffusion coefficient. One can note that the four regions have the same total thickness AZ. We have also computed the MSD in the pure liquids Ll and L2, for comparison purposes. These calculations have been carried out on 864 particles. 205
-0
2
3 June I988
CHEMICAL PHYSICS LETTERS
Volume 147, number 2,3
6
4
8
._
12
10
(I 00 time steps) Fig. 4. Mean-square displacements (MSD) (8X2), (6Y*) and ( 6Z2), along X, Y and Z axes, for the regions II, B 1, Ii?,and B2 (see text and fig.2). Ll and L2 refer to pure liquids. The origins for the MSD along Y and Z axes are shifted respectively by 3006t and 60061, for a clearer graphic presentation. Least-squares fits have been performed in the range [ 100, SOO]& (continuous lines). lime
The computed values of the MSD vary linearly as a function of the time t, in the range [ 100, 500)&t. This corresponds to a diffusive motion. The straight lines, obtained by least-squares fits of the MSD, are plotted in fig. 4. The self-diffusion coefficients, given in table 4, have been computed according to relation (3 ). The diffusion coefficients for the bulk phases
(Bl, B2) and the pure liquids (Ll, L2) are isotropic. They have been calculated by averaging the three X, Y, Z components and are gathered in table 5. The main result is a significative anisotropy of the diffusion coefficients in the interfacial regions (11, 12). Indeed, the diffusion coefficients Dx and Dy are identical and differ from the one measured in the Z direction, Dz. We can thus define a transverse diffusion coefficient, DT, as 4(Dx+ Dy) and a normal one, D,.,, equal to DE These quantities, computed in 11 and 12 regions, are listed in table 5. In our LLI, this anisotropy always corresponds to DTgreater than DN, But the diffusion coefficients in the 11 region are both smaller than the corresponding bulk value (B 1 ), whereas for 12 & and DN are larger than in B2. Let us discuss now the differences between the bulk diffusion coefficients and those in the pure liquids (table 5). They cannot be attributed to the slight miscibilities of each liquid into the other one. These differences come from the respective values of the densities of bulk and pure liquids. This can be shown using the expression of Levesque and Verlet of the self-diffusion coef’licients as a function of the density [151. We consider now the comparison of the transverse diffusion coefficients with the normal ones, both measured in the LLI. Obtaining the highest values in planes parallel to the interface can be related to the decrease of the transverse component of the pressure in the interface (see fig. 3 ). Indeed, a decrease of the pressure increases the self-diffusion coefficient. But unfortunately, in spite of some qualitative explanations [ 16 1, a quantitative prediction is not possible in a simple manner. On the other hand, the ratio DT/DN is equal to 1.11 in the case of 11, whereas its value is 1.29 for 12. This dissymmetry of the interface is also observed in the transverse pressure profile, since its minimum does
Table 4 Diffusion coefficients in the X, Y and Z directions ( 10-j cm2/s) ”
DY DY DZ
11
Bl
Ll
12
B2
L2
3.3f0.4 3.3k0.2 3.0f0.2
3.5kO.4 3.720.4 3.8kO.4
4.0f0.2 4.1f0.2 3.9kO.2
2.2 f 0.3 2.2k0.3 1.7f0.2
1.3f0.3 1.3to.3 1.3kO.2
0.95+0.05 0.9 f0.1 0.92f 0.08
‘) Regions II, B 1,12, and B2 are defined in the text and in fig. 2, and Ll and L2 refer to pure liquids.
206
3 June 1988
CHEMICAL PHYSICS LETTERS
Volume 147, number 2,3 Table 5 Normal and transverse diffusion coefficients ‘) I1 DN ( lO-5 cm*/s) Dr (10-j cm’/s) Dbulr( 10m5cm*/s) (Zp) (lOZocm-‘)
Bl
12
Ll
L2
1.3fO.l 194 k2
0.92 k 0.03 197 ?2
1.7kO.2 2.2kO.2
3.OkO.2 3.3k0.2 176 +2
B2
3.7 r 0.2 175 +2
4.0fO.l 170 f2
184 f2
a) Regions 11, B 1,12, and B2 are defined in the text and in fig. 2, and Ll and 12refer to pure liquids.
not coincide with the average position of the interface (see fig. 3). This minimum is shifted towards the bulk of L2. The dissymmetry is illustrated in fig. 3 by the relative contributions of Ll and L2 to the interfacial tension. The comparison of the normal diffusion coefficient with the bulk one is more complex to analyse, but we can propose a qualitative explanation based on the properties of liquid mixtures [ 17 3. The diffusion in Bl is faster than in B2, and when moving from the bulk phases to the interface, the highest coefficient decreases whereas the lowest increases. This behaviour can be related to the partial interpenetration of the two liquids (see fig. 3). In the inter-facial region one cannot refer to a variation of composition, but nevertheless the neighbourhood of each particle contains a variable percentage of the other species. It is the reason why the faster diffusing species is slowed down while the slower is accelerated.
4. Conclusion MD computations of self-diffusion coefficients in a LLI, using LJ potential functions, allows for the description of mass transport properties in such interfaces. The most important result of this study 1s the anisotropy of the diffusion coefficients of the two liquids in the inter-facial region, the transverse component being greater than the normal one. This anisotropy can be related, at least qualitatively, to the decrease of the transverse component of the pressure tensor. Further efforts are required to quantify these effects.
Acknowledgement
It is a pleasure to thank J.P. Hansen and V. Pontikis for useful discussions and for their help in obtaining computation facilities. M. Mareschal is also acknowledged for his valuable comments. Thanks are also due to J. Mathie for her technical assistance.
References [ I] R. Aveyard and B. Vincent, Progr. Surface Sci. 8 ( 1977) 59. [2] D. Langevin and J. Meunier, J. Phys. (Paris) 44 (1983) CIO-155. [ 31 D. Beysens and M. Robert, J. Chem. Phys. 87 ( 1987) 3056. [ 41 P. Tarazona, M.M. Telo da Gama and R. Evans, Mol. Phys. 49 ( 1983) 283. [5] M. Robert, Phys. Rev. Letters 54 ( 1985) 444. [6] M.M. Telo da Gama and K.E. Gubbins, Mol. Phys. 59 (1986) 227. [ 71 J.S. Rowlinson and B. Widom, Molecular theory of capillarity (Clarendon Press, Oxford, 1982) p. 173. [8] S.M. Thompson, K.E. Gubbins, J.P.R.B. Walton, R.A.R. Chantry and J.S. Rowlinson, J. Chem. Phys. 81 (1984) 530. [9 P. Linse, J. Chem. Phys. 86 ( 1987) 4177. 110 M. Meyer, M. Mareschal and M. Hayoun, I. Chem. Phys., to be published. 111 M. Hayoun, M. Meyer, M. Mareschal, P. Turq and G. Ciccotti, in: Proceedings of Chemical Reactivity in Liquids, Fundamental Aspects, Paris (7-l 1 September 1987), eds. M. Moreau and P. Turq, to be published. 112 H.J.C. Berendsen and W.F. Gunsteren, in: Molecular dynamics simulation of statistical-mechanical systems, eds. G. Ciccotti and W.G. Hoover (North-Holland, Amsterdam, 1986). [ 131 M. Shoen and C. Hoheisel, Mol. Phys. 53 ( 1984) 1367. [ 141 A. Rahman, Phys. Rev. 136A (1964) 405. [ 151 D. Levesque and L. Verlet, Phys. Rev. A 2 ( 1970) 25 14. [ 161 J. Naghizadeh and S.A. Rice, J. Chem. Phys. 36 ( 1962) ‘2710. [ 171 G. Jacucci and I.R. McDonald, Physica A 80 (1975) 607.
207
CHEMICAL PHYSICS LETTERS
3 June 1988
MOLECULAR DYNAMICS STUDY OF SELF-DIFFUSION IN A LIQUID-LIQUID INTERFACE Marc HAYOUN CEN Saclay,D. Tech, 91191 Gif-w-YvetteCedex, France
Madeleine MEYER LPM-CNRS,I placeA. Briand. 92195 Meudon Cedex. France
and Pierre TURQ LaboratoireElectrochimie,UniversitP Pierreet Marie Curie, 8 rue Cuvier, 75005 Paris,France Received 11 March 1988; in final form 29 March 1988
The self-diffusion coefficients have been calculated in a model of a planar interface between two almost immiscible liquids. The simulations have been performed by the molecular dynamics method, using Lennard-Jones potentials. The diffusion process is isotropic in the bulk regions of the simulated system, whereas it exhibits anisotropic behaviour in the interfacial region: the transverse diffusion coefficient (parallel to the interface) is higher than the normal one. A qualitative explanation of this behaviour is suggested by comparison with the pressure tensor.
1. Introduction Most results obtained on the liquid-liquid interface (LLI) deal with static properties investigated from the experimental [ l-31 or from the theoretical [ 1,4-6 ] point of view. Simulation techniques, which have been fruitful for the study of the liquid-vapor interface [ 7,8 1, have not been applied to LLI, with the exception of a Monte Carlo simulation of the water-benzene LLI [ 91 and a molecular dynamics modelling of a LLI [ 10 1. In a previous paper [ 111, we have been able to model an interface between two immiscible liquids, using the molecular dynamics (MD) simulation technique and Lennard-Jones (LJ) pair interactions, In that work, we have shown that the average position and thickness ( = 3 atomic diameters) of this .LLI were stable over the time scale of the simulation. To characterize this interface we have computed the following properties: density and pressure profiles and interfacial tension (20 & 5 erg/cm’). The anal-
ysis of these profiles shows the existence of bulk phases extending over half of the overall volume of the model system. In this paper, using the same model [ 111, we focus attention on the mass transport properties. We have computed the self-diffusion coefficients as a function of the distance to the interface. The results exhibit a diffusional anisotropy in the interfacial region.
2. Model and computations The detailed characteristics of the model have been described [ 111, and the main features are recalled hereafter. The MD simulations have been performed in the (N, V, E) ensemble and the equations of motion were integrated using Verlet’s algorithm [ 121 with a time step 6t = 10 -I4 s. The simulated system contains 864 particles of liquid 1 (Ll ) and 864 particles of liquid 2 (L2). The computational cell is parallelipiped with edges L, L, and i, equal
0 009-2614/88/$03.50 0 Elsevier Science Publishers B.V. ( North-Holland Physics Publishing Division )
203
CHEMICAL PHYSICS LETTERS
Volume 147, number 2,3
3 June 1988
Table 2 MD results for liquid 1 (Ll) and liquid 2 (L2) equilibrated separately ‘) T (K)
P
( 10zocm-))
170 197
133.5 134.2
0.70 0.71
P Ll L2
(kbar)
a1The values for p and T have been averaged over 4000 time steps.
t
a Fig. 1. Schematic drawing of the MD simulation cell. Two adjacent boxes containing respectively Ll and L2 are presented in their initial configuration. The two interfaces a and b are indicated by the arrows; b has been generated by the periodic boundary conditions along Z. LX=&= 37 A, Lz= 69 A. The originof the referenceframeis at the centreofthe computational cell.
Table 1 Lennard-Jones parameters for the system used to simulate the LLI t,tlkr,” (K)
eal/G (K)
e,zlk, (K)
a (A)
100
200
75
3.591
with a IJ pair potential function already used in a study of unstable mixtures [ 131:
where LY- 1, 2 and fi= 1, 2 refer respectively to the atoms of L 1 and L2 (see table 1) . The potential cutoff radius r, is equal to 2.50 and the atomic masses are the same for the two liquids (m = 1.02747 x 1O-22 g). We have shown that this pair potential yields a stable interface between almost immiscible liquids
[111*
a) k, is the Boltzmann constant.
respectively to 37, 37 and 69 8, (fig. 1). In a first stage, we equilibrate separately the two liquids Ll and L2 at the same temperature and the same pressure. The two subsystems are then joined to build the computational cell. The initial position of the planar interface thus obtained is parallel to the XY plane as indicated in fig. 1. Periodic boundary conditions have been chosen. They generate a second interface in the 2 direction (fig. 1). The interactions between particles are calculated
The properties of the two subsystems (Ll, L2) equilibrated separately, are given in table 2. Table 3 contains results of the LLI simulation, including the estimated values of the solubilities between these two liquids.
3. Results and discussion The mass transport properties in the LLI have been investigated through the MD computation of the selfdiffusion coefficients. These quantities have been obtained from the usual least-squares fits of the meansquare displacements (MSD) [ 141. The selfdiffu-
Table 3 MD results for the LLI a) PI
P2
L2inLl (%)
T (IQ
P b’
( 1020cm-3)
Ll in L2 (%)
P
(1020cm-3)
(kbar)
(kbar)
17522
194f2
EO.3
zo.1
135.2
0.72
1.08
a) All these quantities have been averaged over 20000 time steps. p, nd pzreported here correspond to the plateau values of the density profiles. b, Without long-range correction.
204
3 June 1988
CHEMICAL PHYSICS LETTERS
Volume 147, number 2,3
Th
I
12
I1
a
b Fig. 2. Projection in the XZ plane (see fig. 1) of the regions used for the computation of the self-diffusion coefficients. 11:interfacial region in Ll. 12: interfacial region in L2. Bl: bulk region in Ll. B2: bulk region in L2. The four regions have the same total thickness AZ=%0 A.
sion coefficient along the X axis for instance, Dx, is calculated with the following relations: (sX2(r))=~~,~~~~,[Xi(7j+[)-X;(q)]2, T
(2)
(8X2(t))
(3)
=2D,t+const.,
where N is the number of particles taken into account and N, is the number of time origins. In all the computations performed, N, has been chosen equal to 250 and the maximum value of the time t is 5006t. The MSD, along the X, Y and Z axes, have been calculated in four different parts of the system, in order to obtain the diffusion coefficients characteristic of the interface and bulk regions. These four parts consist in slabs parallel to the average orientation of the interface. I1 and 12 are located on each side of the interface respectively in Ll and L2, whereas Bl and B2 correspond to the bulk of each liquid. Since the computational cell contains two stable LLI [ 111, the regions 11 and 12 are each composed of two equivalent slabs of thickness 4AZ, as indicated in fig. 2. The slabs Bl and B2, each of thickness AZ, are located in the middle of the bulk phases ( see fig. 2 ) . For 11, 12, B 1, and B2 the particle numbers were variable. Indeed, in these regions, the particles have been selected every lOOO&,over a run of 100006~. On average, the number of particles of type 1 used
,
0.0 -34.5
I
I
34.5 0.0 z (%I) Fig. 3. Upper part: percentage of L2 as a function of Z, averaged over 2000081 [ 111. Lower part: normal and transverse pressure profiles PN(Z) and PT( Z) (computed without long-range correction). At the second interface (b, see fig. 1) these profiles are extrapolated for Zz34.5 A (open symbols), using the image provided by the periodic boundary conditions. The verical dashed lines refer to the average positions of the interfaces. The contributions of the two liquids to the interfacial tension y are pointed out: Ll (dark areas), L2 (hatched areas). (ya
IV’,(Z)-J’,(Z)ldZ). is 74 for 11 and 119 for Bl, while for 12 and B2 the number of particles of type 2 are respectively 102 and 136. fAZ~2.5 8, has been chosen smaller than the interfacial thickness (Th= 10 A), deduced from the pressure and density profiles [ 111: (i) to take the significant part $Th of the interfacial region (e.g. for L2, between Z1 and Z2 in fig. 3), (ii) 4AZ must be smaller than jTh, the difference between JAZ and f Th being approximately given by (DlOOO&) ‘12, where D is a reasonable estimation of the diffusion coefficient. One can note that the four regions have the same total thickness AZ. We have also computed the MSD in the pure liquids Ll and L2, for comparison purposes. These calculations have been carried out on 864 particles. 205
-0
2
3 June I988
CHEMICAL PHYSICS LETTERS
Volume 147, number 2,3
6
4
8
._
12
10
(I 00 time steps) Fig. 4. Mean-square displacements (MSD) (8X2), (6Y*) and ( 6Z2), along X, Y and Z axes, for the regions II, B 1, Ii?,and B2 (see text and fig.2). Ll and L2 refer to pure liquids. The origins for the MSD along Y and Z axes are shifted respectively by 3006t and 60061, for a clearer graphic presentation. Least-squares fits have been performed in the range [ 100, SOO]& (continuous lines). lime
The computed values of the MSD vary linearly as a function of the time t, in the range [ 100, 500)&t. This corresponds to a diffusive motion. The straight lines, obtained by least-squares fits of the MSD, are plotted in fig. 4. The self-diffusion coefficients, given in table 4, have been computed according to relation (3 ). The diffusion coefficients for the bulk phases
(Bl, B2) and the pure liquids (Ll, L2) are isotropic. They have been calculated by averaging the three X, Y, Z components and are gathered in table 5. The main result is a significative anisotropy of the diffusion coefficients in the interfacial regions (11, 12). Indeed, the diffusion coefficients Dx and Dy are identical and differ from the one measured in the Z direction, Dz. We can thus define a transverse diffusion coefficient, DT, as 4(Dx+ Dy) and a normal one, D,.,, equal to DE These quantities, computed in 11 and 12 regions, are listed in table 5. In our LLI, this anisotropy always corresponds to DTgreater than DN, But the diffusion coefficients in the 11 region are both smaller than the corresponding bulk value (B 1 ), whereas for 12 & and DN are larger than in B2. Let us discuss now the differences between the bulk diffusion coefficients and those in the pure liquids (table 5). They cannot be attributed to the slight miscibilities of each liquid into the other one. These differences come from the respective values of the densities of bulk and pure liquids. This can be shown using the expression of Levesque and Verlet of the self-diffusion coef’licients as a function of the density [151. We consider now the comparison of the transverse diffusion coefficients with the normal ones, both measured in the LLI. Obtaining the highest values in planes parallel to the interface can be related to the decrease of the transverse component of the pressure in the interface (see fig. 3 ). Indeed, a decrease of the pressure increases the self-diffusion coefficient. But unfortunately, in spite of some qualitative explanations [ 16 1, a quantitative prediction is not possible in a simple manner. On the other hand, the ratio DT/DN is equal to 1.11 in the case of 11, whereas its value is 1.29 for 12. This dissymmetry of the interface is also observed in the transverse pressure profile, since its minimum does
Table 4 Diffusion coefficients in the X, Y and Z directions ( 10-j cm2/s) ”
DY DY DZ
11
Bl
Ll
12
B2
L2
3.3f0.4 3.3k0.2 3.0f0.2
3.5kO.4 3.720.4 3.8kO.4
4.0f0.2 4.1f0.2 3.9kO.2
2.2 f 0.3 2.2k0.3 1.7f0.2
1.3f0.3 1.3to.3 1.3kO.2
0.95+0.05 0.9 f0.1 0.92f 0.08
‘) Regions II, B 1,12, and B2 are defined in the text and in fig. 2, and Ll and L2 refer to pure liquids.
206
3 June 1988
CHEMICAL PHYSICS LETTERS
Volume 147, number 2,3 Table 5 Normal and transverse diffusion coefficients ‘) I1 DN ( lO-5 cm*/s) Dr (10-j cm’/s) Dbulr( 10m5cm*/s) (Zp) (lOZocm-‘)
Bl
12
Ll
L2
1.3fO.l 194 k2
0.92 k 0.03 197 ?2
1.7kO.2 2.2kO.2
3.OkO.2 3.3k0.2 176 +2
B2
3.7 r 0.2 175 +2
4.0fO.l 170 f2
184 f2
a) Regions 11, B 1,12, and B2 are defined in the text and in fig. 2, and Ll and 12refer to pure liquids.
not coincide with the average position of the interface (see fig. 3). This minimum is shifted towards the bulk of L2. The dissymmetry is illustrated in fig. 3 by the relative contributions of Ll and L2 to the interfacial tension. The comparison of the normal diffusion coefficient with the bulk one is more complex to analyse, but we can propose a qualitative explanation based on the properties of liquid mixtures [ 17 3. The diffusion in Bl is faster than in B2, and when moving from the bulk phases to the interface, the highest coefficient decreases whereas the lowest increases. This behaviour can be related to the partial interpenetration of the two liquids (see fig. 3). In the inter-facial region one cannot refer to a variation of composition, but nevertheless the neighbourhood of each particle contains a variable percentage of the other species. It is the reason why the faster diffusing species is slowed down while the slower is accelerated.
4. Conclusion MD computations of self-diffusion coefficients in a LLI, using LJ potential functions, allows for the description of mass transport properties in such interfaces. The most important result of this study 1s the anisotropy of the diffusion coefficients of the two liquids in the inter-facial region, the transverse component being greater than the normal one. This anisotropy can be related, at least qualitatively, to the decrease of the transverse component of the pressure tensor. Further efforts are required to quantify these effects.
Acknowledgement
It is a pleasure to thank J.P. Hansen and V. Pontikis for useful discussions and for their help in obtaining computation facilities. M. Mareschal is also acknowledged for his valuable comments. Thanks are also due to J. Mathie for her technical assistance.
References [ I] R. Aveyard and B. Vincent, Progr. Surface Sci. 8 ( 1977) 59. [2] D. Langevin and J. Meunier, J. Phys. (Paris) 44 (1983) CIO-155. [ 31 D. Beysens and M. Robert, J. Chem. Phys. 87 ( 1987) 3056. [ 41 P. Tarazona, M.M. Telo da Gama and R. Evans, Mol. Phys. 49 ( 1983) 283. [5] M. Robert, Phys. Rev. Letters 54 ( 1985) 444. [6] M.M. Telo da Gama and K.E. Gubbins, Mol. Phys. 59 (1986) 227. [ 71 J.S. Rowlinson and B. Widom, Molecular theory of capillarity (Clarendon Press, Oxford, 1982) p. 173. [8] S.M. Thompson, K.E. Gubbins, J.P.R.B. Walton, R.A.R. Chantry and J.S. Rowlinson, J. Chem. Phys. 81 (1984) 530. [9 P. Linse, J. Chem. Phys. 86 ( 1987) 4177. 110 M. Meyer, M. Mareschal and M. Hayoun, I. Chem. Phys., to be published. 111 M. Hayoun, M. Meyer, M. Mareschal, P. Turq and G. Ciccotti, in: Proceedings of Chemical Reactivity in Liquids, Fundamental Aspects, Paris (7-l 1 September 1987), eds. M. Moreau and P. Turq, to be published. 112 H.J.C. Berendsen and W.F. Gunsteren, in: Molecular dynamics simulation of statistical-mechanical systems, eds. G. Ciccotti and W.G. Hoover (North-Holland, Amsterdam, 1986). [ 131 M. Shoen and C. Hoheisel, Mol. Phys. 53 ( 1984) 1367. [ 141 A. Rahman, Phys. Rev. 136A (1964) 405. [ 151 D. Levesque and L. Verlet, Phys. Rev. A 2 ( 1970) 25 14. [ 161 J. Naghizadeh and S.A. Rice, J. Chem. Phys. 36 ( 1962) ‘2710. [ 171 G. Jacucci and I.R. McDonald, Physica A 80 (1975) 607.
207