Nonequilibrium molecular dynamics simulation of diffusion at the liquid-liquid interface

Molecular Dynamics simulations are performed to study the dynamical properties of molecules in the presence of a liquid-liquid (L/L) interface. In the vicinity of the interface the movement of the particles, coupled with the thermal fluctuations of the interface, can make the evaluation of properties such as the self-diffusion coefficient, particularly difficult. We explore the use of the Evans-Searles Fluctuation Theorem [D. Evans and D. Searles, Phys. Rev. E 50, 1645 (1994)] to obtain dynamical information of molecules in distinct regions of a model L/L system. We demonstrate that it is possible to analyse the effect of the interface on the mobility of molecules using a nonequilibrium approach. This information may provide a valuable insight into the understanding of dynamics of interphase mass transfer.


I. INTRODUCTION
3][4] By the same token, their theoretical analysis is usually limited to relatively simple models whose direct comparison with experiments is not straightforward.The ability to model both idealised and more realistic systems using molecular simulations allows one to bring closer theory to experimental results. 5he statistical mechanical description of an interface in a thermodynamic system is necessarily position dependent.Configurational properties, e.g., density and the diagonal components of the pressure tensor, will be a function of the system coordinates. 6,7 is-a-vis, the movement of the particles, coupled with the thermal fluctuations of the interface, can make the evaluation of dynamical properties, e.g., self-diffusion, in a given region of a L/L system difficult to obtain. 8he diffusion properties in interfacial systems have been computed previously, with the majority taking the approach of dividing the system into regions of a given width normal to the interface.For each region, the diffusion coefficients have been computed using either mean square displacements [9][10][11][12] or the velocity autocorrelation functions. 13,14 n a second method, which was only applied to a liquid-vapour system, a criteria for the outmost liquid molecules is defined and their dynamics a) ccorreia@imperial.ac.uk b) a.galindo@imperial.ac.uk c) Author to whom correspondence should be addressed.e.muller@ imperial.ac.uk are followed for a restricted period of time. 15Another method proposed by Liu et al., 16 also applied to a liquid-vapour system, creates a set of virtual boundary conditions and computes the survival probability of a molecule in a region of interest.
An alternative approach was proposed by Duque et al., 17 in which the intrinsic sampling method is applied to analyse the dynamical process at the liquid-vapour interface by tracking the displacements of the surface molecules.All these approaches are based on an equilibrium description of the system.Their computation may be hindered by poor sampling statistics due to either the reduced sample size or time.Some of these limitations are inevitable due to the thermal displacements of molecules in the system leaving and entering zones with different dynamical properties.Others may be related to the thermodynamic nature of the L/L system, where the efficient measurement of diffusion properties using equilibrium methods is restricted to molecules inside their own phase.In this sense, the use of non-equilibrium methods, where an artificial external field is applied on the system to drive the dynamics, might prove fruitful.
0][21][22] These theorems take into account distributions of work along thermodynamically non-equilibrium paths.For systems that are in a steady state not far from equilibrium, the Evans-Searles Fluctuation Theorem can be used to derive Green-Kubo relations for linear transport coefficients. 19hen valid, these expressions are consistent with those resulting from linear response theory. 23he purpose of our work here is to explore the relation between the Fluctuation Theorem and Green-Kubo relations to compute the dynamical properties of a given molecule in the presence of a L/L interface.Taking the diffusion as a representative measurement of the mobility of the molecules at different points in the L/L system, we are able to obtain the dynamical information about the molecules.Our main conclusion is that, using a non-equilibrium approach, not only are we able to obtain valuable information about the response of molecules to a weak field close and far away from the interface but we may also overcome some of the limitations of most equilibrium methods.
The paper is organised as follows.In Sec.II, we describe the methodology and the equilibrium properties of the L/L interface used in this work.Section III is divided in Subsections III A and III B illustrating the equilibrium and nonequilibrium approaches to analyse the dynamical properties of molecules in the presence of a L/L interface.In Sec.IV, we discuss the results obtained from the application of the fluctuation theorem and finally we end with a short summary and conclusion.

II. METHODOLOGY
All simulations were conducted on a system of atoms interacting through a truncated and shifted Lennard-Jones (LJ) 12-6 potential φ(r) with a cutoff radius r cut = 2.5 σ , where r is the scalar distance between a pair of interacting atoms of components i and j, respectively, and φ 0 is the value of the potential at r cut .The corresponding energy and length scales are defined by ij and σ ij , respectively.Unless otherwise indicated, reduced units are used throughout this work, i.e., σ ii = 1.0 and ii = 1.0.In order to obtain a minimally stable L/L interface, a binary interaction parameter δ 12 is introduced that decreases the potential well depth between atoms of unlike components 12 = δ 12 √ 11 22 .Depending on the ratio of the LJ diameters σ 11 /σ 22 and on the well depth 11 / 22 , one may expect the fluids to become immiscible at δ 12 ≈ 0.8. 13,24 n the present work, the L/L interfacial system is modelled using two symmetric fluids, i.e., σ 11 = σ 22 and 11 = 22 , and with a binary interaction parameter δ 12 = 0.2, thus ensuring an extreme immiscibility between the two components in the system.
The L/L interface was modelled using a system of N 1 = N 2 = 2048 particles inside a volume of dimensions L x = L y = 10.0 and L z = 51.2resulting in a global density of ρ = 0.8.The initial configuration was setup by placing the particles of component 1 inside a central volume of dimensions L x = L y = 10.0 and L z = 25.6 surrounded by the particles of component 2. The system was equilibrated for 10 6 steps using a timestep t = 0.001.The temperature was controlled to be T = 0.7 using a using a Nosé-Hoover 25 thermostat with mass m ξ = 50.0.The equations of motion were integrated using the algorithm proposed by Shinoda et al. 26 and implemented in the LAMMPS software package. 27After equilibration, thermodynamic properties were collected over production runs of 10 8 steps.The system is characterised by two immiscible liquids phases: one phase rich in particles of component 1 and the adjacent phase rich in particles of component 2. Due to its periodicity, there are two interfacial regions present in the system.Fig. 1  − z) for phases 1 and 2, respectively.It can be seen that there is no appreciable solubility of either component in the corresponding opposite phase; furthermore, a depletion zone is observed at each interface, i.e., the presence of the opposing phase "dries" the interface.Classical capillary wave theory [28][29][30] provides a means to analyse the average density profile in terms the fluctuations of a bare intrinsic profile, distinct from the equilibrium one. 29he theory describes these fluctuations by assuming the existence of a mathematical collective variable called the intrinsic surface, ξ (R) = ξq e iq•R that represents the instantaneous interfacial surface of a given phase at each transverse position R = (x, y).Different computational techniques have been introduced to enable the calculation of an unambiguous density profile, most notably the so called intrinsic sampling method. 31,32 n here ξ R, q u , is defined in terms of a set of surface atoms with coordinates R up to a wave vector cutoff q u that sets a limit on the resolution of the surface. 32f we take the simplest approximation of representing the intrinsic surface by the long wave limit of its Fourier representation, then we are effectively describing the surface as a sharp interface in the form of a Heaviside step function. 2,33,34 Te density of component 1 on the left hand side can be given by ρ(z) = ρ 1 θ [z − h 1 (x, y)], where h 1 (x, y) is defined here as the instantaneous position of the planar interface of a given configuration of phase 1 and θ (z) is a Heaviside step function.The time derivative of the average density over the (x, y) plane gives the equation for rate of change of the average density profile ∂ρ(z)/∂ z = ρ 1 δ(z − h 1 (x, y)) .If the fluctuations in h 1 (x, y) are small and if we invoke the central limit theorem, 2,34 then the delta functions may be approximated by a corresponding Gaussian function of the same width δ(z )) and can be integrated to give where h 1 and w 1 are the average position and width of the interface of component 1, respectively.Similar expressions may be derived for each one of the interfaces illustrated in  I gives a measure of the uncertainties in the calculations, as the system is symmetric with respect to components 1 and 2.

III. DIFFUSION IN THE PRESENCE OF A LIQUID-LIQUID INTERFACE A. Equilibrium molecular dynamics
Previous studies of the molecular mobility in an inhomogeneous system employ the approach of dividing the system into regions of time-varying particle density ρ(z) and computing the velocity autocorrelation function of the particles in each region.For the system studied here, slabs of half width ε k are taken in the z direction and we consider velocities normal (v z ) and tangent (v x , v y ) to the interface.Because the number of particles on each slab is a fluctuating variable, these need to be accounted accordingly in the calculation of the correlation function where k is the slab index, z k is the slab position, and ε k is half the slab width.The summation is over the number of particles N a of component a, v i is the velocity of particle i in the corresponding Cartesian direction, and δ i, k (t) samples the number of particles of type i that are in slab k at time t.The product δ i, k (0)δ i, k (t) samples only the particles that are in slab k both initially and after some time t.The self-diffusion coefficient in each slab region is then estimated from the Green-Kubo integral of the corresponding velocity autocorrelation function where integration time is limited to the residence time of the particle in the corresponding region.In this sense, the above equation may be seen as an approximate measure of the diffusion coefficient.A total of eight regions are used to divide the system and the corresponding velocity autocorrelation function of each component is computed for each slab.The resulting self-diffusion coefficients in each slab are collected in Table II and illustrated in Fig. 2. Close to the interface, an anisotropy between the diffusion coefficients parallel  13,14 with disparate conclusions.Mayer et al. 10 and Buhn et al. 14 studied the L/L interface between two simple LJ systems and both noted an increase of the tangential diffusion coefficients when approaching the interface, as observed in the present work.Mayer et al. 10 argued that such enhancement is due to a decrease in the values of the components of the pressure tensor tangent to the interface (P xx , P yy ).In contrast to the latter, Benjamin 8 did not observe any anisotropy in the diffusion coefficients across the L/L water/1,2-dichlorothane system.This behaviour was attributed to the fact that this particular system has a more structured interface where the orientational dynamics are more restricted as a result of the long range interactions between the water molecules. 8n Fig. 3, the Cartesian components of the correlation functions of component 2 close to the interface are presented.In addition to the observed break up in isotropy (Fig. 3 top), a local maxima of the velocity autocorrelation function is present at t ≈ 0.37, as illustrated by the corresponding derivative (Fig. 3 bottom).This is unique to the correlation functions computed near the interface and is an indication that a particle experiences an acceleration not seen in the bulk phase.At the simulated temperature T = 0.7, the distance travelled by a particle over this period is approximately δz ≈ t × √ kT /2 = 0.22 which is commensurate with the half average interface width w α = 0.49.Although it is fair to assume that the interface is playing a role in the latter case, either through the presence of a depletion zone in the density of particles ρ(z) or through an anisotropy in the pressure components, this analysis may only yield a coarse view of the dynamics of the particles.The reasons are twofold: first, the complication associated with diffusing molecules entering and leaving the different domains means that Eq. ( 4) can only be evaluated for a finite amount of time; second, the presence of capillary wave fluctuations of the interface is inducing a blurring in the properties of slabs.

B. Non-equilibrium molecular dynamics
The Green-Kubo relations express the transport coefficients as integrals over corresponding time correlation functions. 36,37 nder appropriate conditions, these coefficients may also be evaluated from an analysis of the relevant fluctuating property through Einstein's relation.The connection between Green-Kubo and Einstein's relation in a nonequilibrium steady state has recently been established. 23,38 is relation makes use of the fluctuation theorem which relates the probability of observing trajectories of some duration t to the corresponding conjugate ones through a characteristic dissipation function t . 22,39 he dissipation function is de-fined as a dimensionless dissipated energy that is integrated along the system's trajectory and it represents a generalised form of entropy production. 38,39 nsider a dynamical system described by Eq. ( 5) where a colour charge s i is assigned to a subset of particles, coupling them to an external field of strength F e .The coupling coefficient is termed color to emphasise the point that it does not enter into the inter-particle interactions, but only defines how the particle interacts with the external field. 40The method used here resembles an experimental technique where the mobility of a probe molecule is studied using a controlled external magnetic field. 41The time evolution of the N-particle system under periodic boundary conditions is described by where q i and p i are the position and momentum of the ith particle with mass m i , s i = (−1) i − 1 for i ≤ 2 and zero otherwise, F e is the external field strength, ξ is the thermostat degree of freedom with mass m ξ = 50.0and β = 1/kT, with T = 0.7.The thermostat is only applied to the remaining (N − 2) particles in the system, i.e., it only acts on the atoms that are not coupled to the external field.These latter atoms constitute a heat bath to the two probe particles and ensure that these reach a non-equilibrium steady state.
For completeness, we summarise the main expressions of the Evans-Searles Fluctuation Theorem relevant to this work.The reader is referred to the original papers for an unabridged version. 22,23 or the system described by Eq. ( 5), the dissipation function t 23 takes the form where the dissipative flux is given by J ( ) = −(s 1 v z,1 + s 2 v z,2 ), the relative displacement between the probe particles is given by z 12 (t) = z 12 (t) − z 12 (0) with z 12 = z 1 − z 2 .For this system, the transient fluctuation theorem reads If the system reaches a steady state and if it is unique, 19,23,39 then when the dissipation function t is evaluated for deterministic steady state trajectories, the theorem may hold asymptotically, When the applied field is sufficiently weak, then after long enough times, at least compared with the characteristic Maxwell relaxation time of the dissipative flux, one may expect the contributions to the average of the dissipative flux J (t) = −v z,12 (t) to become statistically independent and to obey the central limit theorem, 19,23,38 i.e., the distribution would become a Gaussian near the mean As expounded by Sevick et al., 38 the time average dissipative flux J (t) = − v z,12 (t) increases with the field strength.
For weak enough fields where the mean is not too large compared with the standard deviation, Eq. ( 9) may be expressed as F e (11)   from which an Einstein relation may be obtained for the steady state velocity of the probe particle 1 relative to particle 2, In the zero-field limit, the variance of the distribution of v z,12 (t) is independent of the direction and magnitude of the field F e and to leading order in F e , it may be expressed as where the subscript • Fe=0 represents the field free average and it is assumed that the equilibrium average of the steady state velocity between the probe particles is zero v z,12 (t) F e=0 = 0.The equilibrium variance of the time averaged inter particle velocity is given by where ξ = t 1 − t 2 .In the long time limit, or at least at some time larger than the correlation time required for the time correlation function to decay to zero, 23 the second term on the right hand side of Eq. ( 14) may be neglected.Substitution of Eq. ( 14) into Eq.( 12) gives The maximum value of the field value F e is defined by the time it takes for the system to reach a steady state and the central limit theorem to be valid. 42If the integration time t in Eq. ( 15) is large enough compared with the characteristic relaxation time of the particles but small enough such that the local thermodynamic properties in the neighbourhood of that point do not change markedly, then Eq. ( 15) may provide a means to measure the diffusion properties at different points in the vicinity of the liquid-liquid interface.The objective of this work is therefore to explore the possibility of measuring the diffusion coefficient at different points in the L/L system using Eq. ( 15).

IV. RESULTS AND DISCUSSION
A series of simulations are set up with probe particles of component 2 starting at representative points along the L/L system described in Fig. 1 and Table I.Given that the average position of each L/L interface is in the vicinity of h ≈ ±13.0, four different sets of simulations were performed with probe particles of component 2 starting at fixed points along the z axis of the L/L system, namely, z(0) = {18.0,16.0, 14.0, 12.0}.For the first three coordinates the probe particles start at points belonging to their own phase (see Fig. 1) but for z(0) = 12.0, the probe particles are initially in the phase belonging to component 1.We take advantage of the symmetry of the system where for each particle of component 2 considered, a mirror image particle of the same type is also pulled in the opposite direction, i.e., z i (0) = (−1) i |z(0)| and s i = (−1) i − 1 hence both particles start from symmetric positions and are driven in directions opposite to each other.The set of field strengths used are F e = {−2.0,−1.5, −1.0, 1.0, 1.5, 2.0}.For each field value starting from a given coordinate point z 0 a set of 2 × 10 4 independent trajectories are simulated.We use Eq. ( 5) from which averages were taken over periods of t = 10.0 using a timestep t = 0.001.From each set of trajectories the distributions functions of the inter-particle distance p( z 12 (t)) were computed.
The computed distributions of trajectories as a function of field strength for particles starting at z(0) = 18.0 and z(0) = 14.0 are depicted in Fig. 4. Positive field strengths will push the particles in the direction of the interface and negative ones will pull the particles away from it.In a bulk region away from the interface, z(0) = 18.0, we can observe that the computed distributions are similar, with a Gaussian shape, differing only in the average value.
From the symmetry of the responses to the field direction, we may conclude that there is no discernible influence of the interface on the particles in the region z 0 ≈ 18 ± 2.0.This result is consistent with the equilibrium results illustrated in Figs. 1 and 2 where only for values close to the interface anisotropy in the diagonal components of the pressure tensor is observed.In other words, from Fig. 4 one may observe that a molecule in the bulk region has responses which are similar regardless of the applied field and which differ only in the sign of the average distance travelled according to the direction of the driving force.Different force strengths only affect the average distance travelled.Similar results (not shown) can also be observed when the departing position is z(0) = 16.0.Closer to the interface, when the probe particles start from z(0) = 14.0, there is a clear difference between the steady state responses to positive and negative field directions.Moving away from the interface, the observed distributions have similar variances.However, when the field direction is positive and drives the particles towards the other phase, their mobility is reduced as a result of the potential energy barrier associated with the interfacial region.
The average velocity at time t defined in Eq. ( 15) is calculated from the first moment of the distribution where p( z(t)) is the normalised distribution function of the particle displacement at time t for a given field strength.The field dependent mobilities may then be obtained from the application of Eq. ( 15) when the systems reaches steady state.In Figs. 5 and 6, the results obtained from application of Eq. ( 8) are presented for non-equilibrium trajectories starting at z(0) = 18.0 and z(0) = 14.0, respectively.The absolute value of A is a measure of the displacement of the particles.Away from the interface, at z(0) = 18.0, the fluctua- tion theorem is quickly satisfied for relatively short times, t = 2.0, and no significant change can be observed at longer times, t = 10.0.However, when the trajectories start from z(0) = 14.0 (cf.Fig. 6, F e = −1.5 and F e = −1.0),although their response reaches a steady state at short times consistent with the local environment near the interface, it can be observed that their response becomes transitive at later times.This behaviour seems to be a blueprint of a particle probing the heterogeneity of the system by sampling two regions with distinct thermodynamic properties.A physically compelling explanation stems from the observation that at the end of the simulation, the probe particles have travelled an average distance of z 12 (t = 10.0)/2 ≈ 1.0 (see Fig. 4).This distance is, on average, long enough for the probe particle to reach the bulk region (see Fig. 1).The new environment will affect the dissipation rate of the system as the probe particles leave the region near the interface.
In Table III, the resulting diffusion coefficients obtained by application of Eq. ( 15) at different times and for different initial positions are presented.In steady state and in the linear regime, the diffusion coefficient should be independent of the applied field.This is the case of trajectories starting at z(0) = 18.0 where there is an asymptotic convergence of the slopes of the curves as time increases.Closer to the interface, for trajectories starting at z(0) = 14.0, this linear dependence breaks down, with the initial diffusion rates being significantly larger than the ones obtained at longer times.Comparison with the equilibrium results depicted in Fig. 2 This article is copyrighted as indicated in the article.Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 155.198.12.107On: Fri, 17 Oct 2014 13:56:07 shows that, within the statistical error, the mobility values obtained from non-equilibrium simulations are consistent with the corresponding equilibrium ones for different points in the L/L interface.It can be observed that the diffusion coefficients computed from equilibrium simulations are systematically smaller.This difference may be ascribed to the fact that the integration time of the Green-Kubo integral Eq. ( 4) is limited to the residence time of the particles inside the region of interest.Increasing the region size will allow longer time intervals but it would also lead to a stronger smearing effect, particularly in regions closer to the interface.In contrast to the trajectories where the particles started from points within their own phase, the case where z(0) = 12.0 merits a separate discussion.When the trajectories start from this position, the probe particles are, at least on average, outside their own phase.The probe particles are in an energetically unfavored state and experience a highly anisotropic force pulling  them towards their own phase.The influence of the interface overwhelms the applied external field with the end effect that the average displacement of the probe particles is negative towards the interface.
Einstein's relation Eq. ( 15) may not be applicable as is, because in this case there is an additional field F s derived from the presence of the interface that acts on the particles.If we assume that the interfacial field strength to be constant within the space travelled by the molecules, then these may still reach a steady state and the analysis leading to Eq. ( 15) may be substituted by where F e is the applied external field strength and F s is the surface force experienced by the particle.In Fig. 7, the resulting mean time displacements of the inter-particle distance as a function of the applied force are depicted.Although the mean displacements are systematically negative as shown by the corresponding distribution functions in Fig. 8, the relation between v z,12 (t) and the applied field is still linear.From the slope of these curves it is still possible to obtain the diffusion properties using Eq. ( 17) albeit with a non-zero coordinate value.The asymptotic value of the diffusion coefficient is then given D(t = 10.0)= 0.035 in good agreement with the equilibrium results depicted in Fig. 2.This article is copyrighted as indicated in the article.Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 155.198.12.107On: Fri, 17 Oct 2014 13:56:07

V. CONCLUSIONS
We have examined the diffusion properties of particles within a model liquid-liquid system with particular emphasis on the interfacial region.In agreement with previous work, results based on equilibrium simulations expose an anisotropy in the self-diffusion coefficients normal and tangential to the interfaces.
We have investigated the application of the fluctuation theorem and its relation with the Green-Kubo relations for the analysis of the diffusion properties in a model liquid-liquid system.It was shown that analysis using the steady state fluctuation theorem provides detailed information on the dynamical properties of the molecule at different points of the liquidliquid system.In addition, the fluctuation theorem allows the characterisation of the diffusion properties of molecules outside their own phase, on the other side of the interface, a property that is difficult to compute using equilibrium methods, where the asymmetry in forces experienced by the molecules avoids the accurate calculation of the diffusion properties.
The fluctuation theorem, through an analysis of the dissipative work necessary to drive the probe molecule, gives the possibility of describing the dynamics of the molecule at, and across, the interface through its free-energy profile.Although outside the scope of the present work, an in-depth description of the details of the interfacial region would be necessary for the non-equilibrium analysis of the interface dynamics and the correct measurement of the potential of mean force.Such approach would give a valuable insight into processes involving more realistic models of complex fluids, where the more ordered nature of the liquid phase and its interplay with the orientational entropy of the probe molecules are of considerable interest.

− 30 −FIG. 2 . 2 FIG. 3 .
FIG. 2. Tangent D x , D y and normal D z self-diffusion coefficients plotted as a function of the slab height z.The diffusion coefficients of components 1 and 2 are represented by crosses and circles, respectively.For each component, the symbols for D x and D y overlap.The density profiles of each component are overlaid as a guide.This article is copyrighted as indicated in the article.Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 155.198.12.107On: Fri, 17 Oct 2014 13:56:07

FIG. 4 .
FIG. 4. Normalised distribution functions of the response of the inter-particle distance z 12 (t) to different field strengths F e at t = 10.0 for trajectories starting at z(0) = 18.0 (top figure) and at z(0) = 14.0 (bottom figure).

FIG. 5 .
FIG. 5. Steady state fluctuation theorem plots for trajectories starting at z(0) = 18.0 using different field strengths.The solid line is the prediction of the steady state fluctuation theorem in Eq. (8).Top figure depicts the results at t = 2.0 and bottom figure depicts the results at t = 10.0.

FIG. 6 .
FIG.6.Steady state fluctuation theorem plots for trajectories starting at z(0) = 14.0 and using different field strengths.The solid line is the prediction of the steady state fluctuation theorem in Eq. (8).Top figure depicts the results at t = 2.0 and bottom figure depicts the results at t = 10.0.

− 10 − 8 − 6 − 4 −FIG. 8 .
FIG.8.Normalised distribution functions of the response of the inter-particle distance z 12 (t) to the external field of different strengths F e at t = 10.0 when the particles start at z(0) = 12.0.
depicts the average density profiles ρ(z) = δ(z i

TABLE I .
Structure of the liquid-liquid interfaces depicted in Fig.1.ha,kis the average position of interface k with reference to component a and w a,k is the interface width.In TableI, the resulting interface width and average position of each component, respectively, are collected.The departure from symmetry in the results collected in Table This article is copyrighted as indicated in the article.Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 155.198.12.107On: Fri, 17 Oct 2014 13:56:07 154101-3 Braga, Galindo, and Müller J. Chem.Phys.141, 154101 (2014)

TABLE II .
9,10,35ffusion coefficients D for each component in the directions tangent (x, y) and normal (z) to the interface as a function of the position z, ε k is half the slab width.The interfaces are located approximately at z ≈ −12.55 and z ≈ 13.05, respectively.,Dy) and normal (D z ) to the interface are observed.This anisotropy has been studied previously by Meyer et al.,9,10,35by Benjamin 8 and more recently byBuhn et al.