Laterally driven interfaces in the three-dimensional Ising lattice gas
We study the steady state of a phase-separated driven Ising lattice gas in three dimensions using computer simulations with Kawasaki dynamics. An external force field acts in the direction parallel to the interface, creating a lateral order parameter current which varies with distance from the interface. Above the roughening temperature, our data for ‘shear-like’ linear variation of are in agreement with the picture wherein shear acts as effective confinement in this system, thus supressing the interfacial capillary-wave fluctuations. We find sharper magnetisation profiles and reduced interfacial width as compared to equilibrium. Pair correlations are more suppressed in the vorticity direction than in the driving direction; the opposite holds for the structure factor. Lateral transport of capillary waves occurs for those forms of for which the current is an odd function of , for example the shear-like drive, and a ‘step-like’ driving field. For a V-shaped driving force no such motion occurs, but capillary waves are suppressed more strongly than for the shear-like drive. These findings are in agreement with our previous simulation studies in two dimensions. Near and below the (equilibrium) roughening temperature the effective-confinement picture ceases to work, but the lateral motion of the interface persists.
pacs:05.40.-a, 05.50.+q, 68.05.Cf, 68.35.Rh
Dimensionality of space is a relevant parameter in condensed matter systems that exhibit large spatial fluctuations, e.g., thermal fluctuations of the order parameter near a second-order phase transition. Such critical fluctuations are correlated over large distances, which gives rise to singular behaviour characterized by a set of critical exponents whose values depend on . Moreover, the very existence of the phase transition depends on . Below the lower critical dimension fluctuations are strong enough to destroy the ordered phase, and hence there is no phase transition. In contrast, above the upper critical dimension fluctuations are no longer important, and the critical exponents become the same as in mean field theory.
An important dependence on dimensionality also occurs in systems where two distinct phases coexist. In such systems thermal fluctuations can be correlated over a long distance within the interfacial region. Then the interfacial correlation length (parallel to the interface) increases with increasing thickness of the interface, , where Fisher (), and the interface is termed rough. The value of the roughness exponent again depends on dimensionality . A statistical-mechanical description of interfacial degrees of freedom, capillary wave theory, was proposed by Buff, Lovett and Stillinger BLS (); BedeauxWeeks (), who suggested to model the interface as a sharp divide between the two phases, but one that would freely fluctuate. Capillary wave theory (CWT) predicts that for , i.e., fluctuations can destabilise the interface in . Then the thickness of the free interface is infinite in the thermodynamic limit. At the marginal dimension , the value corresponds to a logarithmic divergence. For , one has but is finite for all ; the interface is then said to be smooth.
These results originate from the choice of the simplest, Gaussian form of fluctuations, i.e., the probability for a local departure (height) of the interface from the reference plane is given by the Boltzmann weight , where is the thermal energy, and
where are coordinates parallel to the interface. is the interfacial stiffness, which for a continuum fluid is simply the tension of a free interface Fisher (); LipowskyFisher (); Jasnow (); Tarazona (). If the external potential is quadratic in , for example due to gravity, the equipartition theorem for quadratic degrees of freedom may be applied. The equilibrium interface pair correlation function for a translationally-invariant system is then found to be
where , and the limit of infinite interface dimension , has been taken. An upper cutoff on wave numbers is always assumed. is usually identified with some appropriate microscopic length in the interface region, e.g., the lattice spacing or the bulk correlation length Weeks77 (). The behaviour of is obviously -dependent – as is the interface width , defined by .
For certain microscopic models, such as the Ising model, which is equivalent to the lattice gas model of a fluid, exact results are available DBA_PTCP (); Jasnow (). A description starting from a microscopic Hamiltonian is particularly useful because it accounts for both interfacial and bulk fluctuations. In contrast bulk degrees of freedom of coexisting phases are absent in CWT, in which one considers only interfacial configurations. On the other hand, in the “classical” theories for the interface based on order parameter (or density) profiles, for example the van der Waals theory vdW (); RowlinsonandWidom (), interfacial fluctuations are absent and the structure of the interface is reduced to an inhomogeneity of the order parameter, i.e., it is of order of the bulk correlation length .
In , Ising interfaces are rough for all temperatures below the critical temperature , and , as revealed by exact results for the local magnetisation profile on a scale large compared to Douglas (), in agreement with CWT. However, in there is evidence Jasnow (); WGL (); Weeks80 () for the existence of a finite roughening temperature , below which the interface is smooth. This evidence is supported by Monte Carlo (MC) simulations HP () and by rigorous analysis of discrete random surface models, such as the solid-on-solid (SOS) and discrete Gaussian (DG) models Weeks80 (); Whiskers (). These models approximate the interface of the Ising model at low temperatures, for example the SOS model can be regarded as a certain anisotropic-coupling limit of the Ising model where the SOS interface configurations are selected from the Ising configurations by the requirement that there are no overhangs or bubbles DBA_PTCP ().
In the situations that we shall address, the interface is stabilised by the presence of two walls at spacing . This problem is interesting because the fluctuating interface will experience collisions with the constraining walls. The standard capillary wave model does not apply for this case; it has to be extended to take into account entropic repulsion from the walls. Confined Ising interfaces can be treated rigorously in . Results for the magnetisation profile Maciolek () indicate that in two dimensions, in spite of the entropic repulsion at its extremities, the interface meanders back and forth between the walls so that . This is markedly different from what is expected in three dimensions on the basis of the analysis of low-energy excitations in discrete random surface models Whiskers (). An energy-versus-entropy argument leads to the conclusion that interface configurations which result in interference with the boundary are needle-like. A rigorous analysis gives and . These results support conjectures from the phenomenological effective interface Hamiltonian Fisher (); FandF (); LipowskyFisher (), and for Ising interfaces they agree with MC simulation studies KKB ().
Dependence of the structure of equilibrium interfaces on the spatial dimensionality has repercussions also on relaxation dynamics of the fluctuating interface. Relaxation dynamics of capillary waves have been recently studied using molecular dynamics simulations of simple liquids in Tarazona (); Thakre (). For the liquid interface constrained between two walls at separation , a pronounced enhancement of the relaxation time of capillary waves was found; the most affected are relaxation times of waves with the wave number Thakre ().
These results for equilibrium suggest that a non-trivial dependence of the interfacial structure and dynamics on dimensionality may persist to non-equilibrium situations. In this paper, we study this problem for fluctuating interfaces that are driven into a steady state by the action of an external field parallel to the plane of the interface. The problems of roughness, spatial and temporal correlations of laterally driven interfaces were first addressed in the lattice gas driven by a spatially uniform force field via MC simulations LeungMon89 (); LeungZia93 (). These studies, which were of two-dimensional systems, found that the interface becomes less rough when drive is applied. The order-parameter (magnetisation) profiles become much sharper upon increasing the drive and the interfacial width is reduced as compared to the equilibrium value. By a suitable coarse graining of microscopic particle configurations, the local position (height) of the interface was defined, and the behaviour of the spatial interface height correlation function was studied. The results are consistent with a reduction of the interfacial correlation length . Moreover, the structure factor , defined as the the Fourier transform of the height-height correlation function displays deviations from the equilibrium capillary wave dependence as ; the data suggest a weaker singularity , which implies a reduction of the roughness exponent . In theoretical attempts to treat driven interfaces, one derives a dynamic equation for the interfacial degrees of freedom, starting from the time-dependent Landau-Ginzburg-type equation for the order parameter Leung88 (); ZiaLeung91 (). This approach leads to a non-local and non-linear equation for the interface height. A linear stability analysis of this equation for a spatially uniform drive parallel to the interface shows that temporal decay of fluctuations along the driving field is faster than that orthogonal to the driving field. However, predictions for the roughness of the interface do not agree with the simulation results in .
More recent interest in the theoretical challenges of laterally driven interfaces originates from an experiment on colloidal gas-liquid interfaces subjected to shear flow. In the experiment of Derks et al. DerksShear () a real-space visualisation of interfacial fluctuations revealed a reduced interfacial roughness when shear was applied. The width and correlation function of the height of the interface were analysed by fitting to the equilibrium CWT results. The fitting parameters in the analytic CWT expressions were the correlation length along the interface measured in the flow direction, and the surface tension. The authors concluded that was decreased but was increased. Recently, the problem of non-equilibrium fluctuations of a liquid-liquid interface under shear has been addressed theoretically ThieBick () within the framework of fluctuating hydrodynamics. This approach leads to a mode-coupling equation for the interface height which was solved using a perturbation theory. Results for the interfacial width are in agreement with the experiment of Ref. DerksShear (), but the results for the interfacial correlation length in the flow direction are not. The theoretical height-height correlation function and the structure factor imply a decrease of the correlation length in the direction of flow. Interestingly, in the direction perpendicular to the flow (vorticity direction), the correlation length seems to increase. Similar conclusions have been obtained from molecular dynamics simulations Thakre ().
Previously we have studied interfaces in the two-dimensional Ising strip driven by an external field that is applied parallel to the walls (and to the interface), and may vary in the direction perpendicular to the mean position of the interface, by using MC simulations with spin-exchange Kawasaki dynamics Shear2dPaper (); LM7Proc (); MovPaper (). These studies were partially motivated by the need to understand sheared fluid interfaces. Because our results were obtained in , and because of the simplified character of our model and its dynamics, we were not in the position to attempt a direct comparison with experimental data. However, our results were in partial qualitative agreement with Ref. DerksShear (). We found that the shear-like drive acts as an effective confinement on the system; a steady state is reached in which the magnetisation profile is the same as that in equilibrium, but with a rescaled length implying a reduction of the interfacial width. Pair correlation functions along the interface decay more rapidly with distance under drive than in equilibrium, and for cases of weak drive can be rescaled to the equilibrium result. Moreover, we find that interfacial transport can occur in an unexpected way parallel to the interface. The lateral flux of the order parameter at a planar interface induces lateral motion of the thermal capillary waves, provided that the flux is an odd function of distance from the interface.
In the present paper we study the same model system using the same approach, but in three spatial dimensions. We wish to investigate to what extent our findings from persist to higher dimensions. Moreover, dependence on dimensionality of various quantities characterising the structure and dynamics of the interface, as well as of transport properties in a driven state is interesting. One would like to know whether the different character of interfacial fluctuations, i.e., “wandering” in and “spikes” Whiskers () in has any repercussions. There are also new aspects that deserve to be studied: the behaviour of the system near the roughening transition, and the anisotropy effects introduced by the introduction of the vorticity direction. Last but not least, results in permit a more direct comparison to experiment.
The rest of this paper is organised as follows. In Section II we introduce the model and give details of the simulations. In Section III.1, we present and discuss results for the interfacial structure of the driven Ising system, via the magnetisation profile, interface width, and correlation functions in real and Fourier space. In Section III.2, we investigate the dynamics of the driven interface, showing results for the current and evidence for capillary wave transport. Finally we draw conclusions in Section IV.
Ii Model and simulation details
We consider a three-dimensional (3d) Ising model on a simple cubic lattice. On lattice sites sit spins , which may take values . In lattice gas language, the spin variables become particle occupation numbers , corresponding to the absence and presence of a particle at a site, respectively. The Hamiltonian for the system is
where indicates that the sum is over nearest neighbour sites and . The coupling constant , so that the interactions are ferromagnetic (attractive in lattice gas language).
The lattice has dimensions , with periodic boundary conditions applied in the and directions. (All lengths are expressed in units of the lattice constant ). Spin layers are located at heights , for a total of layers. The interface is induced by walls of fixed spins at the top () and at the bottom () planes of the lattice; these boundary conditions energetically favour parallel alignment of the interface with the plane, with the ‘+’ phase in the upper half of the volume, . We focus on slab-like lattice geometries, with and , so that the system is confined between the two walls, and the scaling length scale for the interfacial width is LipowskyFisher (); Whiskers ().
Time evolution of the system proceeds under Kawasaki Kawasaki () spin-exchange dynamics, which conserve the magnetisation, or equivalently the number of lattice-gas particles, locally. Elementary simulation moves consist of exchanging the values of two randomly chosen nearest-neighbour spins with probability ; in the lattice gas, this corresponds to a particle moving to a neighbouring empty lattice site. The key features of these dynamics are their conservation of order parameter, and their locality. These attributes are desirable from the point of view of simulating particle movement (albeit crudely) and the competition between diffusive motion and driven motion.
In the general case, an external force field acts on the system, driving in the plane. This field alters the Monte Carlo acceptance rates, to produce a modified Metropolis rate,
Here, is the inverse temperature (the Boltzmann constant is set to unity), and is the change in internal energy from the proposed exchange. is the work done by or against the external force field; for , the above rate reduces to the standard Metropolis one, which samples thermal equilibrium states. We are interested in the case of non-zero , when the system will reach a non-equilibrium steady state. The system is immersed in a heat bath at constant temperature , into which the work done is dissipated. The driving field is related to the work term by
where , , and the displacement vector between spins and is or . In the following we will consider the forms of the driving field which are perhaps the most relevant experimentally. Our main focus is the case of ‘shear-like’ linear variation of driving field with , and the field acting in the -direction only, such that the field components are , . Thus exchanges along are enhanced or suppressed, while exchanges in the and directions proceed with equilibrium rates (). We also study the case of a V-shaped spatial dependence, , , such that the drive acts in the same direction throughout the system. Some results for a spatially uniform driving field in the direction, , , and a step-like field, , are also included.
We have carried out extensive Monte Carlo (MC) simulations of the above model using single-spin and multi-spin NewmanBarkema (); GemmertBarkema () coding techniques, in the latter case generalizing the driven multi-spin algorithm used previously to 3 systems. The multi-spin method we have adopted allows simulation of 64 independent systems (and hence greatly enhanced statistics) on a 64-bit computer system, using efficient bitwise operations. The state of an Ising spin may be represented by one bit; thus the values of a particular site in 64 different systems can be stored in a 64-bit variable (the natural word size). Bitwise operations operate on all bits of a variable at once, and are computationally cheap instructions. By combining these operations appropriately, and generating random bits with the required probabilities, the desired acceptance rates can be produced.
Single-spin results provide a check as to the correctness of the multi-spin implementation. Parallelization via domain decomposition was employed to speed up simulations; the lattice was sub-divided along the direction, and appropriate synchronization used when exchanging spins on and near the boundary between two domains. Data processing was also parallelized, owing to the large quantity of data produced by the multi-spin algorithm. Results reported here are from total run lengths of MC sweeps ( trial moves), the slow evolution of Kawasaki dynamics to a steady state proving to be less of a problem in than in LM7Proc ().
The majority of the results shown here are for a system size and or , at a temperature , where () is the bulk critical temperature of the equilibrium 3 Ising system IsingTc3d (). We have checked that going to larger values of has a minor effect on the results only for the smallest considered wall separation , i.e., as long as the longitudinal correlation length , which grows exponentially with , is less than . We have also varied the temperature, firstly to investigate the effect of an increase to , and secondly to study the behaviour near and below the roughening transition. For an equilibrium system in the thermodynamic limit, the roughening temperature is HP (). In a finite system, the (pseudo-)transition will occur at a shifted value of temperature, which will be governed by the system size BarberFSS (). As temperature is increased in the smooth regime, either the interfacial width will reach the scale of , or the lateral correlation length will reach or – in either case, the interface reaches the rough regime. We have thus covered a range of temperatures around the roughening temperature in the simulations, from to . The roughening transition belongs to the universality class of the Kosterlitz-Thouless transition KT-trans (). The renormalization group method of Kosterlitz Kosterlitz () showed that the correlation length diverges very rapidly as from below:
where and are non-universal parameters. The numerical values for these constants obtained from MC simulation studies of the Ising interface in are and HP (). The shift of the pseudo-roughening temperature can be estimated from the condition , which yields . This is a very weak dependence on and for our system size it gives . At the same time the width of the interface diverges upon approaching the bulk roughening temperature, as . The condition yields , which is a much stronger dependence on the finite dimension than that obtained from the condition involving . For , , therefore we conclude that the shift of the roughening transition is governed by . Using that estimate of gives , so the range of simulation temperatures should include the equilibrium pseudo-transition temperature.
We first investigate the interfacial structure of the driven 3 Ising model in order to see whether also in three dimensions the effective action of drive is to increase the confinement of the interface.
iii.1.1 Magnetisation profiles
The magnetisation profile along the axis is calculated as
where the angles denote an average in the steady state. In a phase separated system changes sign across the interface, and attains values close to near the upper and lower walls respectively. In Fig. 2 we plot the magnetisation as a function of the scaled variable . Upon applying either shear-like or V-shaped drive, the magnetisation profile becomes ‘sharper’: changes sign more rapidly in the interfacial region, and there is a more extended flat region near either wall. The size of this effect increases with increasing . These trends are the same as in two dimensions Shear2dPaper (), but for given driving strength, we find that the magnitude of the effect is smaller in 3.
It is possible to rescale the driven profiles to collapse back onto the equilibrium result: see Fig. 2. We interpret this as the drive acting to reduce the effective distance between the walls from to , and thus to effectively increase the confinement of the system. This increase may be quantified by introducing the scaling , where is the ratio of the effective and actual wall separations. For shear-like drive with , we find , whereas for two dimensions Shear2dPaper (), for the same value of ( is named there), . One may more formally express this in terms of a finite size scaling function for the profile. The scaling relation for an equilibrium fluid confined between two walls was proposed by Fisher and de Gennes FdG (). In the absence of a bulk field it reads:
where is the bulk correlation length, and is the spontaneous magnetisation in bulk. and are finite-size scaling functions: is obtained from by changing the first scaling variable . Thus in equilibrium the shape of the scaling function can be varied by changing the wall separation at fixed or, equivalently, by changing the temperature at fixed . We find that driving changes the shape of the interfacial profile at fixed temperature and such that
where is a boundary correction that decays away from the walls on the scale of . Relation (8) admits (in the scaling regime) another interpretation of the result (9), namely, as the drive acting to reduce the effective temperature of the system at fixed actual distance between the walls. The estimation of the rescaling factor is obtained by rescaling the driven data to spline-interpolated equilibrium curves and minimizing the associated chi-squared statistic – the value of at the minimum is the optimal value. For the system, the 8 points in the centre of the system were included in the rescaling procedure – points further out are subject to stronger wall interaction effects. Both ways of driving yield very similar rescaling factors; more pronounced differences are observed for smaller values of and lower temperatures with the conclusion that effective confinement is stronger for the V-shaped drive. Rescaling fails for lower temperatures closer to and below the equilibrium bulk roughening transition temperature.
The magnetisation profile may also been used to study the behaviour of the interfacial width. We measure the width via the second moment of and study its variation with driving strength, wall separation , and temperature. Upon increasing driving strength or , the width reduces, as expected from the results for the full profile. For shear-like drive, we are also able to obtain data collapse for the behaviour of as a function of a scaling variable . Here is an adjustable exponent. The division of the width by corresponds to the expected equilibrium behaviour LipowskyFisher (); Whiskers (); KKB (), so that for , . The scaling behaviour of the width is shown in Fig. 3, for fixed temperature , and a variety of wall separations and drive gradients in the ranges , . Previously we obtained data collapse as a function of the same scaling variable in the 2 system Shear2dPaper () (there the width was scaled by the equilibrium behaviour, ). Remarkably, we obtain data collapse for the same value of exponent as for the system. From Fig. 3 we see that for small at the larger values of and 20, the data collapse is lost – we believe that this is because one moves out of the confined regime with for these parameters. For these wall separations, the longitudinal correlation length becomes comparable to the linear dimension of the interface, and the system crosses over to the regime where the dominant length scale is . This does not require a large increase in , because , where the transverse length scale for Gaussian interface fluctuations Fisher (); KKB (). Data collapse is regained for larger values of , because the effective wall separation is the controlling length scale (from the discussion of the magnetisation profile above), and is small enough for the system to be in the “confined regime”. The inset of Fig. 3 shows the variation of the width with drive gradient for shear-like and V-shaped drive. The trends are rather similar, with the width for given very slightly smaller for V-shaped drive – this is consistent with the previous conclusion that the confinement is stronger for this drive type.
iii.1.2 Spin-spin correlation functions
In order to characterize the behaviour of the sytem on the two-body level, we first consider the microscopic interface pair correlation function, i.e., the spatial spin-spin correlation function at the midplane, defined as
which depends on separations in both the and directions. Comparing the specific cases and provides information on the anisotropy in the and directions. In equilibrium, for , but for driven systems, the two functions may differ. Results in Fig. 4a show that shear-like drive causes both and to decay more quickly and for larger separations to saturate at larger asymptotic values than in equilibrium. In the direction, this finding is in agreement with hydrodynamics results ThieBick (); however in that study, the correlation length in the direction was found to increase under shear, contrary to the trend in our system. We defer further exploration of this difference to the discussion of the height correlations below, since the height variable provides a more direct point of comparison between the systems.
As in the 2 case, the spin correlation functions at intermediate separations may be transformed to the equilibrium result via a rescaling of the lateral coordinate, or : and Shear2dPaper (); see the inset of Fig. 4a for rescaling results for . The rescaling factors are obtained via the same method as for the magnetisation profile; in this case, very small values of or are cut-off in the procedure, as are the tails of the functions, so that the rescaling procedure is carried out over . The parameters may be interpreted as ratios of lateral interfacial correlation lengths in and out of equilibrium: . Rescaling the driven results produces and , so that the correlation length is reduced under drive – this tallies with the faster decay evident from Fig. 4a, which shows results for equilibrium and shear-like drive. The anisotropy may be measured by the ratio ; this is consistently slightly smaller than unity, leading to the surprising conclusion that the correlations are slightly more suppressed in the (vorticity) direction. As with the magnetisation profile, the effect of drive is much weaker in three dimensions than in two – for example, in 2, for , while in 3, for : a ten times larger field gradient is required to produce a comparable confinement. As a result, the rescaling procedure works for much larger values of than in 2 – a stronger drive is required to push the system into a regime which is too far from equilibrium for rescaling to be possible. The V-shaped drive has a similar effect on the interfacial correlations: for example, with , , corresponding to a slightly stronger confinement effect than with shear, consistent with the findings for the magnetisation profile.
Fig. 4b shows the effect of varying the temperature on the interfacial correlations. Lowering the temperature to (below the roughening transition in a bulk equilibrium system) in equilibrium results in correlations that decay only to for the largest separations. Since on the lattice the average interface position lies between two lattice points (for zero overall magnetisation), one measures the correlations just either side of the interface (which side does not matter, due to symmetry). Thus at zero temperature, for all , since the interface is perfectly flat at . This explains the observed increase of asymptotic values of for low . Moreover, the width of smooth interfaces is of order of the bulk correlation length, which at low temperatures is 2-3 lattice spacings. Therefore, for low temperatures essentially measures correlations in the bulk-like phase. We also see from Fig. 4b that driving the system enhances the asymptotic value further for , which indicates that at fixed temperature the bulk-like phase is more ordered under drive than in equilibrium. As for the magnetisation profile, rescaling does not work for the low temperatures – the drive affects the asymptotic value more than the decay rate for these temperatures.
iii.1.3 Height-height correlation functions
We now turn to the interfacial height-height pair correlation function. The interface height is obtained via a coarse-graining procedure, which produces a single-valued height function from the real microscopic configuration. The latter contains ‘bubbles’ of one phase in the other, and ‘overhangs’ at the interface, meaning a height cannot be defined from it directly. To coarse-grain, we use a simple method which we found to be successful in the 2 system, where it gives results equivalent Shear2dPaper () to a more complicated coarse-graining method DeVirgiliis (). The height is simply a sum of the spins over a column. Thus when there are equal numbers of ‘+’ and ‘-’ spins in a column, , while when there is a majority of one species, . The general height-height correlation function depends on spatial separations , and on temporal displacement :
where the angles indicate an average over time. We first consider the equal-time correlations, with one of the spatial separations set to zero: Fig. 5 shows results for . In two dimensions, (also ) in equilibrium exhibits strong anti-correlated regions for medium-to-large separations, presumed to be finite size effects LM7Proc (). These are not present in for the system sizes considered – the functions decay to zero without becoming signifcantly negative, indicating less severe finite size effects; an explanation may be the following. With conservative dynamics, a positive-height ‘bump’ must be accompanied by one with negative height, since . In , these must lie on the same layer (the only one), so an anti-correlation is measured in . However in there are layers, so the pair may be located in different layers, meaning does not necessarily display anti-correlations. Turning to the driven cases, we see that applying shear-like drive leads to more rapid decay of , as well as a smaller initial value ) (a measure of the interfacial width). The magnitude of this effect increases with increasing . In this section, the results shown are for shear-like drive, but the conclusions also apply to V-shaped drive – as for other static quantities, the effect is similar to shear, with slightly greater correlation-suppression.
As for other quantities, rescaling to the equilibrium result is possible in the same manner as in . For the height correlations, the rescaling takes the form . The values of and are those obtained from the rescaling of the magnetisation profile and the spin-spin correlation function for given simulation parameters. In this procedure was motivated by Weeks scaling in equilibrium BedeauxWeeks (): , where is the scaling function, and since in , the correct rescaling involves . Although Weeks scaling does not hold in , this procedure works reasonably well for , except at zero separation, where the rescaled values become larger than the equilibrium width. As with the spin correlations, this range of validity is much greater than in two dimensions.
Furthermore, we are able to fit the results for and for small-to-intermediate separations to the equilibrium capillary wave result for the height correlation function in ; see Fig. 6 for results for shear-like drive. This procedure was used previously by Derks et al. to describe their experimental data, where an excellent fit was obtained DerksShear (). By integrating Eqn. (2), the (equilibrium) capillary wave result for in and in the limit is BedeauxWeeks ()
where is the modified Bessel function of the second kind. The upper wave-number cutoff has been sent to infinity in order to obtain an analytic result; in order to regularize the integral at , a shift is introduced. Bedeaux and Weeks BedeauxWeeks () give , where is the original wave-number cutoff in the integral (2). Combining (12) with the capillary-wave result for the interfacial width,
we are able to substitute for the (unknown) in terms of the width, and obtain a fitting form with two parameters: the correlation length and the pre-factor . In the above, we have specialised to separations in rather than the radial distance usual in CWT, since isotropy is broken in the non-equilibrium situation. For separations in , the form for is the same, but different values of the parameters are expected – i.e., the correlation length () will be different, as will the pre-factor. The intepretation of the latter quantity is difficult. Indeed, the interface tension is an equilibrium concept and cannot be carried over directly to non-equilibrium situations, so the meaning of the pre-factor is not initially clear – here we just note its anisotropy.
The equilibrium fit wanders off the data for larger separations; this may be due to a (less serious) manifestation of the finite size effects encountered in , which were mentioned above, and the conserved order parameter. (This is most obvious for the equilibrium data on the log scale in Fig. 6, where the data diverge as they approach zero and become negative). For the driven cases, as drive becomes stronger, the fit works for a smaller range of separations – the example of shear-like drive is given in Fig. 6. This trend is expected from the findings for the rescaling procedures applied above – initially the system is “close enough” to equilibrium for CWT to be approximately applicable, but as drive increases, this ceases to be true. From the fits we obtain the equilibrium and non-equilibrium correlation lengths in the and directions, and – see the inset of Fig. 6 for their variation with . The trend of decreasing correlation length with increasing drive strength mirrors the one found in , although there we were not able to obtain the correlation length reliably, due to the difficulty of fitting the correlation function data over a reasonable range. We also note that is consistently smaller than , in agreement with the earlier conclusions, based on the behaviour of the spin-spin correlation functions, that correlations are slightly more strongly suppressed in the vorticity direction than in the driving direction.
The suppression of correlations we find in the drive () direction is in agreement with the hydrodynamics work of Thiébaud and Bickel ThieBick (), who studied phase separated fluids between two walls under shear. This trend is also the same as in the 2 Ising system, and both microscopic (, spin-spin) and coarse-grained (, height) measures of correlations give the same conclusion. Both the theoretical and simulation findings disagree, however, with the experimental results of Derks et al. DerksShear (), who found an increase of correlation length in the flow direction when shear was applied to a phase-separated colloid-polymer mixture. The fact that we have used the same method of fitting the height correlation data to the equilibrium capillary wave form as Ref. DerksShear (), makes the method of comparison the same, at least.
Finally, we consider the pre-factor resulting from the fit of the height correlation data to the CWT form (12). In equilibrium, the pre-factor is proportional to , where is the surface stiffness, which for a continuum fluid is the interfacial tension – the free energy associated with the interface. Out of equilibrium, this quantity is not defined, so the meaning of the pre-factor resulting from the fit is not clear. Numerically, we find that the pre-factor from fitting is a decreasing function of drive gradient for shear-like drive. If one defines a “non-equilibrium surface tension” from the CWT fit, and further takes the temperature to be fixed (i.e., the value of the parameter in simulations), the conclusion is then that this “tension” increases as the system is more strongly driven. This procedure of defining an effective non-equilibrium surface tension via the CWT fit was the approach taken in the analysis of the experiments of Derks et al. DerksShear (), who also found this tension to be an increasing function of shear rate. The CWT fits of experimental data in Ref. DerksShear () show that an increase in a “non-equilibrium surface tension” is accompanied by an increase in correlation length, but our simulations show the opposite relationship – a decreasing correlation length as the effective interfacial tension increases. Our findings seem to be inconsistent with the CWT result (see Eqn. (2) and below) . However, in our system at equilibrium , so that the effective increase of the confinement due to driving (reduction of ) wins over the effective increase of .
iii.1.4 Structure factor
To further complicate the situation, our finding of a decrease of correlation length in the vorticity () direction is in disagreement with Ref. ThieBick (), where an increase was found. Intruigingly, Thiébaud and Bickel found the structure factor to be unaffected in the direction by the application of shear ThieBick (). The same was concluded for the uniformly driven system from the analytic approach based on the time-dependent Landau-Ginzburg functional by Leung Leung88 (). The static structure factor in our system is accessible via a two-dimensional spatial Fourier transform of the equal-time height correlations, :
where denotes a two-dimensional spatial Fourier transform, , and , so that and lie on the range . In Fig. 7 we plot for equilibrium and driven systems, along either the or direction, as a function of . From Eqn. (2), in equilibrium, one expects . In Fig. 7 we fit the equilibrium data to this form. The data shown are along the direction with , although we have checked that the equilibrium structure factor behaves the same along the direction, as expected. For , this behaviour is indeed observable in the simulations, but for , the data diverges from the CWT prediction, when other powers of presumably become important.
Turning to the non-equilbrium behaviour, we see that for shear-like drive, is affected (suppressed) in both the drive (, blue crosses) and vorticity (, red squares) directions, but the effect is smaller in the vorticity direction. These results are in disagreement with the hydrodynamics results ThieBick (); however, since the effect in the drive direction is stronger than in the vorticity direction, the latter effect could possibly be of higher order than was considered in Ref. ThieBick (). The data for shear-like drive along are also fit to the equilibrium CWT form in Fig. 7; we see that as for equilibrium, the fit is reasonable for . Additionally, the intercept at is greater, indicating a smaller lateral correlation length, as found in real space above. Indeed, one can compare the parameters resulting from the CWT fits in real space and Fourier space. We find that the qualitative trend for the correlation length is the same, but do not obtain quantitative agreement – the values obtained from the real space fits are consistently larger. These differences are expected – for equilibrium, they can be caused by the finite system size and lattice discretization effects. For non-zero drive, the effect of deviations from CWT can be different in real and Fourier space. Additionally, the fits in Fourier space are for small (long wavelengths), while the real-space fits are for small separations, so the length scales the fits are applicable to is not necessarily the same. For V-shaped drive, we find that for given drive gradient , the results are similar to those for shear-like drive, with slightly greater suppression of the structure factor at small .
iii.2 Capillary wave transport
We now consider the dynamics of the height-correlation function defined in (11). The primary limit of interest is which we find shows evidence of capillary wave transport along the driving direction, for suitable forms of the current profile. Previously (in the context of ) MovPaper (), we conjectured that the capillary wave fluctuations on an interface will be transported by an external driving field, provided that the lateral order parameter current has a component which is an odd function of distance from the interface. Thus only purely even current profiles are expected to give no transport.
iii.2.1 Order parameter current
In the system, the order parameter current profile is defined as , where is the current profile of the spin species. The components and of are the net number of spins of that species moving in positive and directions per unit time at perpendicular coordinate . We also note that the Ising symmetry means that the order parameter current can also be written as ; the first definition may be applied to systems lacking the Ising symmetry, i.e., liquid-gas or liquid-liquid interfaces.
As shown in Fig. 8, the shear-like drive , gives a purely odd order parameter current component . Since the component of the field is zero, . For the V-shaped drive, the current is an even function of , since . We thus expect transport along for the shear-like drive, but none for the V-shaped drive. For the same value of , the currents for the two drive types almost coincide in the region , where the driving fields are the same in magnitude and direction.
In Figs. 8 and 8, is shown for various drive gradients , for temperatures above and below the (equilibrium) roughening transition. Looking at the high temperature data in Fig. 8 we see that for small , has maxima at the walls; as is increased, plateaus develop with the maximum current shifted slightly from the wall. Eventually for strong drive the maxima become localized near the interface. This reflects the competing effects of local drive strength and current carrier availability ( pairs): for large , the drive strength is essentially saturated at the walls, so the greater carrier density at the interface eventually becomes more important. Below the bulk roughening temperature (, Fig. 8), the current also has maxima at the walls for small , and quickly develops maxima at the middle two layers as is increased. These maxima appear for much weaker drive (approximately six times smaller ) than they do for . They are also localised to the two middle layers either side of the interface, and are much more pronounced than at the higher temperature; this indicates that at low temperatures the interface region is very sharp, reduced to approximately two lattice spacings. For strong drive, the greater carrier density at the interface again ‘wins’, and these maxima become global. We also note that is roughly five times smaller than that at the higher temperature, since the carrier density is much smaller due to the increased bulk and interfacial order.
Finally, Fig. 8 also shows an example of mixed symmetry in the current profile. The driving field is of the ‘V’ type, but with different values of in the upper and lower halves of the system: in the lower half, in the upper. Thus the total driving field can be written in the form , with , , showing the even and odd components explicitly. The current profile reflects the asymmetry in the drive: in the lower half of the system, the current matches that for a (symmetric) V-shaped drive with , while in the upper half, it matches that for either V or shear with (see Fig. 8) – the crossover occurs over a single lattice spacing.
iii.2.2 Space-time correlations.
We investigate whether the conjecture for the occurrence of capillary wave motion holds for systems by measuring for different forms of driving field, which produce differing current profiles. For current profiles with odd symmetry in (or more generally, profiles with an odd component), we expect to see evidence of capillary wave transport in . Fig. 9 shows that this is indeed the case – see the main panel for results for for shear-like drive, and the inset for V-shaped drive. For time difference , the peak lies symmetrically around due to the translational invariance ensured by the periodic boundaries along . However, at time differences , the peak moves to negative values for shear-like drive, indicating that now the greatest correlations are between spatially-displaced points. We interpret this to mean that wave-like height fluctuations are being coherently transported along the interface by the drive. For the V-shaped drive, the peak remains at for all , showing the absence of wave motion. In both cases, correlations decay with increasing time difference, due to thermal noise. We note that the rate of decay of correlations is much faster for driven systems than it is for equilibrium systems with Kawasaki dynamics.
We have also investigated other forms of driving field, for example spatially uniform driving field in the direction, , and step-like drive: . The former produces an even order parameter current profile whereas the latter an odd one, and results for (not displayed) show that wave motion does not occur for a uniform drive but occurs for a step-like drive, consistent with the conjecture. For the case of the asymmetric V-like drive discussed in the previous section, which has mixed symmetry, we expect to see wave movement, since the current profile is not purely even, but like the driving field itself, can be written as a sum of even and odd components. Indeed we find this is the case, with the peak of moving with time. From these results we conclude that the criterea for capillary wave motion are the same in the and Ising systems.
Having established the occurence of wave motion, it is natural to investigate the dependence of the wave velocity on system parameters. We measure the speed of the peak of , , and vary the driving strength ( for shear-like drive, for step-like drive), temperature and wall separation . As shown in Fig. 10, shows linear variation with or for fixed temperature and system size, for . For shear-like drive, the gradient of is close to 2 in this range. We also see that varying has a rather small effect on the peak velocity – doubling from 10 to 20 reduces the gradient of by only a few percent. Changing the temperature from to also has small effect, in the other direction – the peak moves faster for the higher temperature at given , . For the step drive, also seems to be linear in driving strength for small , with a reduced gradient compared to shear-like drive. For both forms of drive, non-linearity appears to set in for . For the mixed symmetry case, where the driving field could be written as , we find that the velocity of motion is smaller than that for a purely odd field with – the velocity in the mixed case is linear in and approximately 80% that of the pure case for . For lower temperatures near and below the equilibrium roughening temperature, we find that the interfacial motion still occurs, with a much reduced velocity; correlations also decay much more quickly with time.
iii.2.3 Dispersion relation
In order to characterize the dynamics of capillary waves we consider the evolution of the spatial Fourier modes of the height function defined by:
where and are two-dimensional vectors and the division is defined as . (Recall that and are integer coordinates of the lattice.) The sum denotes summations over integer and from 1 to and , respectively. Because is a real function, the Fourier transform has the conjugate symmetry , i.e., there are independent terms in (15). Each complex Fourier component can be written in terms of its modulus and phase :
where . If , each mode would correspond to a travelling wave moving in the direction with a velocity and (with no dispersion).
In a steady state, the phase shift of each mode is a fluctuating quantity, with a distribution which when measured with increasing time intervals spreads and decays quickly to zero. However, at short times, we are able to measure its mean value in unit time to obtain the dispersion relation of the frequency as a function of the wave vector as:
Results for calculated for equal to of an MC sweep are shown in Fig. 11 as a function of for several values of . We observe similar behaviour for all considered temperatures, including below the roughening temperature (data not shown). The shape of is very similar to that obtained in the two-dimensional system (see Ref. MovPaper ()), which suggests to use the same analytical formula for describing the dispersion relation:
with small- expansion
Indeed, as can be seen in the Fig. 11, our data for fit well with . Waves with larger are less dispersive in the sense that the corresponding coefficients and are smaller. The fastest mode is the one with , or a wavelength of two lattice spacings; can also be fitted using (18) with . As we pointed out in Ref. MovPaper (), if the dynamics of the height function are modelled by the following linear transport operator (propagator)
then the Ansatz in the form of the travelling wave yields the first and the second terms in the fit function (18). In (20), time is treated as a continuous variable whereas spatial derivatives are discrete: , and a 5-point stencil is used for . The third term, which for small gives the quadratic dependence in (19), can be obtained if one allows an imaginary contribution (with 3- and 5- point stencils) to the linear operator . One may be able to understand the presence of complex coefficients in the transport equation for the height function , by recognising that the plane wave solution above neglects the dependence of the amplitude on the wave number. In fact the average modulus of each complex Fourier component varies significantly with , even in the absence of driving (see the plot for the structure factor Fig. 7). Taking this into account, it should be possible to derive an equation for the amplitude as well as the phase – this may aid in understanding the presence of real and imaginary parts in the transport operator. The limit of the solution of the amplitude equation should yield the static structure factor, which may be compared to simulation data. We leave the interesting and difficult problem of deriving the full transport equation from these considerations to future work.
Iv Conclusions and outlook
We have presented evidence from Monte Carlo simulations that the main characteristics of interfacial structure and dynamics in the laterally driven stochastic lattice gas model are generic and persist to three dimensions. Far above the equilibrium roughening temperature the structure of the interface confined between two walls is affected by lateral driving in a way similar to that resulting from an increase of confinement (by reducing the distance between the walls) of the equilibrium system. However, the effect of drive is much weaker in three dimensions than in two, plausible due to the different nature of low-energy interfacial fluctuations, i.e., the “spike”-like excitations rather then interface “wandering”. In the case of the shear-like drive, our findings for the decay of the height-height correlation functions and for the structure factor are in partial agreement with recent results from fluctuating hydrodynamics ThieBick (). The discrepancy concerns the behaviour in the vorticity direction. We have found a decrease of correlation length in this direction whereas hydrodynamic calculations predict an increase. Also, our data show that the structure factor is suppressed both in the drive and the vorticity directions. In Ref. ThieBick () it is concluded that is unaffected in the vorticity direction. Moreover, our results for the interfacial width are in agreement with the experiment of Ref. DerksShear (), but the results for the interfacial correlation length in the flow direction are not. It would certainly be desirable to carry out more studies, experimental, theoretical and simulation-based, to clarify these discrepancies. Nearer (and also below) the bulk roughening temperature, the confinement/ordering effect of the drive is reduced, since there are no large fluctuations to “smoothen out”. The picture of an equilibrium system under a greater effective confinement no longer seems to apply for these low temperatures.
The conjecture made in Ref. MovPaper () for the occurence of lateral transport of interfacial fluctuations is supported by our results in three dimensions. The transport also occurs for low temperatures below the equlibrium roughening transition. However, the dynamics of the lateral propagation of capillary waves which gives rise to the observed dispersion relation is still not understood. A full treatment should include non-linear effects, and treat the time and wave-vector dependence of the Fourier amplitude. Further insight, especially into the effect of coupling between bulk and interfacial degrees of freedom, could be gained by considering theoretical models based on the order parameter bray ().
- (1) See, for example, M. E. Fisher, in Fluctuations, Interactions and Related Transitions, Jerusalem Winter School for Theoretical Physics, edited by D. Nelson, T. Piran, and S. Weinberg (World Scientific, Singapore, 1989), Vol. 5.
- (2) F. P. Buff, R. A. Lovett and F. H. Stillinger, Phys. Rev. Lett. 15, 621 (1965).
- (3) D. Bedeaux and J.D. Weeks, J. Chem. Phys. 82, 972 (1985).
- (4) R. Lipowsky and M. E. Fisher, Phys. Rev. B 36, 2126 (1987).
- (5) D. Jasnow, Rep. Prog. Phys. 47, 1059 (1984) and references therein.
- (6) R. Delgado-Buscalioni, E. Chacon, and P. Tarazona, Phys. Rev. Lett. 101, 106102 (2008); J. Phys.: Condens. Matter 20, 494229 (2008).
- (7) J.D. Weeks, J. Chem. Phys. 67, 3106 (1977).
- (8) D. B. Abraham, Structure and Phase Transitions in Surfaces, in Phase Transitions and Critical Phenomena, Vol. 10, ed. by C. Domb and J. Lebowitz (Academic, London, 1986), and references therein.
- (9) J. D. van der Waals, Z. Phys. Chem. 13, 657 (1894).
- (10) J. S. Rowlinson and B. Widom, Molecular Theory of Capillarity, (Oxford: Oxford University Press, 1982).
- (11) D. B. Abraham and P. Reed, Phys. Rev. Lett. 33, 377 (1974).
- (12) J. D. Weeks, G. H. Glimer, and H. J. Leamy, Phys. Rev. Lett. 31, 549 (1973).
- (13) J. D. Weeks in Ordering in Strongly Fluctuating Condensed Matter Systems, edited by T. Riste (Plenum Press, New York, 1980) p.293.
- (14) M. Hasenbusch and K. Pinn, J. Phys. A: Math. Gen. 30, 63 (1997) and references therein.
- (15) J. Bricmont, A. El Mellouki and J. Frohlich, J. Stat. Phys. 42, 743 (1986).
- (16) J. Stecki, A. Maciołek and K. Olaussen, Phys. Rev. B 49, 1092 (1994).
- (17) M. E. Fisher and D. S. Fisher, Phys. Rev. B 25, 3192 (1982).
- (18) T. Kerle, J. Klein, and K. Binder, Phys. Rev. Lett. 77, 1318 (1996).
- (19) Amol K. Thakre, J. T. Padding, W. K. den Otter, and W. J. Briels, J. Chem. Phys. 129, 044701 (2008).
- (20) K.-t. Leung, K. K. Mon, J. L. Vallés and R. K. P. Zia, Phys. Rev. Lett. 61, 1744 (1988); Phys. Rev. B 39, 9312 (1989).
- (21) K.-t. Leung and R. K. P. Zia, J. Phys. A: Math. Gen. 26, L737 (1993).
- (22) K.-t. Leung, J. Stat. Phys. 50, 55 (1988).
- (23) R. K. P. Zia and K.-t. Leung, J. Phys. A: Math. Gen. 24, L1399 (1991).
- (24) D. Derks, D. G. A. L. Aarts, D. Bonn, H. N. W. Lekkerkerker and A. Imhof, Phys. Rev. Lett. 97, 038301 (2006).
- (25) M. Thiébaud and T. Bickel, Phys. Rev. E 81, 031602 (2010).
- (26) T. H. R. Smith, O. Vasilyev, D. B. Abraham, A. Maciołek, and M. Schmidt, Phys. Rev. Lett. 101, 067203 (2008).
- (27) T. H. R. Smith, O. Vasilyev, D. B. Abraham, A. Maciołek, and M. Schmidt, J. Phys. Condens. Matter 20, 494237 (2008).
- (28) T. H. R. Smith, O. Vasilyev, A. Maciołek and M. Schmidt, EPL 89, 10006 (2010).
- (29) K. Kawasaki, Phys. Rev. 145, 145 (1966).
- (30) M. E. J. Newman and G. T. Barkema, Monte Carlo Methods in Statistical Physics (Oxford University Press, 1999).
- (31) S. van Gemmert, G. T. Barkema and S. Puri, Phys. Rev. E 72, 046131 (2005).
- (32) A. L. Talapov and H. W. J. Blöte, J. Phys. A 29, 5727 (1996).
- (33) M. N. Barber, in Phase Transitions and Critical Phenomena, edited by C. Domb and J. L. Lebowitz (Academic, London, 1983), Vol. 8, p. 145.
- (34) J. M. Kosterlitz and D. J. Thouless, J. Phys. C 6, 1180 (1973).
- (35) J. M. Kosterlitz, J. Phys. C 7, 1046 (1974).
- (36) M. E. Fisher and P. G. de Gennes, C. R. Acad. Sci. Paris Ser. B 287, 207 (1978).
- (37) A. De Virgiliis, E. V. Albano, M. Müller and K. Binder, Physica A 352, 477 (2005).
- (38) A. J. Bray, A. Cavagna and R. D. M. Travasso, Phys. Rev. E 64, 012102 (2001); ibid. 65, 016104 (2001).