0 Copyright 1994 American Chemical Society
The ACS Joumal of
Surfaces and Colloids FEBRUARY 1994 VOLUME 10, NUMBER 2
Letters Monte Carlo Simulation of Curvature-Elastic Interfaces P. Pieruschka' and S. MarEelja Department of Applied Mathematics, Australian National University, Canberra 2601, Australia Received June 14, 1993. In Final Form: August 19,199P We introduce a new Monte Carlo simulationtechnique for the equilibration of curvature-elasticinterfacial systemssuch as surfactant f i b s dissolved in bulk media. The method is based on a continuousrepresentation of the interfaces and can accurately evaluate all relevant surface averages of the curvature Hamiltonian. We apply the method to sponge-like surfactant systems with the standard Hamiltonian, where the only parameters are the bending and saddle-splay moduli, K and ii. Random bicontinuous surface states are found to be stable for low bending rigidity and small negative saddle-splay modulus, justifying the use of Gaussian ensemblesas approximationsto real interfacialensembles. Topologicalchangesof the random sponge states as a function of the ratio iilK are analyzed; the results over a wide range of K suggest that only this universal ratio determineswhether the final states resemble disordered minimal surfaces,disordered lamellar, or connected parabolic geometries. 1. Introduction
Self-assembly in systems of oil and/or water in the presence of amphiphilic surfactant often leads to mesoscopic structures of bulk media separated by a thin film of surfactant. In the case of oil and water the film is monomolecular, while for a single solvent the surfactant forms a bilayer. In both cases the phase diagrams are veryrich. Possible geometries of the surfactant film range from simple dispersions of vesicles and stacked lamellar structures to topologically complex disordered bicontinuous spongegeometries or crystalline interfaces resembling periodic minimal surfaces. The mechanisms which lead to the formation of a certain geometry are only partly understood, but it is widely accepted that the relevant Hamiltonian is the curvature-elastic energy of the surfactant film and that interfacial tension can often be neglected.' There is a large variety of models for microemulsions. Descriptions as microscopic spin systems2 where two or *Abstract published in Advance ACS Abstracts, December 15,
1993. (1) Row, D.; Codon, C.; Cates, M. E. J. Phys. Chem. 1992,96,4174. (2) Gompper,G.; Schick,M. Phys.Rev. B 1990,41,9148,and references therein.
three spin values are assigned to mark oil, water, and surfactant molecules are very general in that they make no assumptions about the film geometry. However, the representation of the bending interaction in terms of spin interactions is difficult. Effective interface models which work directly with ensembles of interacting interfaces use the curvature-elastic Hamiltonian, which should be good as long as the monomer concentration is low.1~~ For the case of random sponge structures, an effective interface model has recently been linked to the BerlinKac mean-spherical approximation4 for two-state spin systems s(P) = fl. The continuous Gaussian spin field s ( i ) can be related to an ensemble of surfaces which are located at the interfaces of the positive ( s ( i ) > 0, e g . water) and negative ( s ( i ) < 0, e.g. oil) regions of the continuous spin field.5-8 (3) Cates, M. E.; Row, D.; Andelman, D.; Milner, S. T.;Safran, S. A. Europhys. Lett. 1988,5,733. (4) Berlin, T. H.; Kac, M. Phys. Reu. 1952,86,821. (5) Berk, N. F. Phys. Rev. Lett. 1987,58, 2718. (6) Teubner, M. Europhys. Lett. 1991,14,403. (7) Pieruschka,P.: MarEelja, S. J. Phys. ZI 1992,2,235. In eq 2 of this papera fador of 112was inadvertentlyomkted reaultingin the o v e d h a t e of the bending modulus by a factor of 2. The corrected values are in much better agreement with experiment. (8)Pieruschka, P.; Safran, S. A. Europhys. Lett. 1993, 22, 625.
0743-7463/94/ 2410-0345$04.50/0 0 1994 American Chemical Society
Letters
346 Langmuir, Vol. 10, No. 2, 1994 This theory derives the free energy directly from the curvature Hamiltonian for balanced microemulsions and bilayer systems (with zero spontaneous curvature) 7fC= 2K
ss(+)
c,+c
2
d S + ii
ss ss
clc2dS =
H = -(uV2s + VSVU) and aii
K = Z
;(i-&J
aii
(4)
with u = l/m and ii = uVs. The values of H and K can be calculated analytically from the readily known first 2~ H2 dS + ii KdS (1) and second derivatives of eq 3. The remaining problem is to determine the position of the surface, which changes where c1 and c2 are the local principal curvatures, and H with every update, and to measure its area. This is done and K are the mean and saddle-splay curvatures of the by dividing the sample volume V into subcubes small surface with total area S. K and ii are the bending and enough to allow for linear approximation within each saddle-splay moduli. It results in the free energy subcube and following the field gradient in several iterations until the surface is reached with good accuracy. (Alternatively,to measure the area one could also distribute small subcubes randomly within V, but for convenience we prefer simple subdivision. In principle, however, the where ( ) denotes ensemble averaging over the Gaussian method is completely lattice independent i.e. a lattice is field s$tes, k, is of the order of the inverse molecular size, not needed to define the surface. All information is given and v(k) is the spectral density of the Gaussian distribution by eq 3.) The approximate surface area S(2i) on site i is of Fourier amplitudes. From the free energy, all relevant evaluated by polygon tiling as in ref 11. Numericalresults thermodynamic quantities can be calculated and the have been compared with analytical expressionsfor regular results indicate the utility of Gaussian random interfaces surfaces and analytical results for ensemble averages of as model systems for soft elastic interfacial systems."** Gaussian random surfaces. (Relative errors for grid size Nevertheless it is useful to examine the foundations of the 32 per unit cell for the P-surface (PS)and Gaussian random Gaussian model in more detail and understand the limits surface (GRS)with v(k) = 6(k - ko) for surface to volume of its applicability to physical systems. ratio, Gaussian, and mean square curvatures are in The Gaussian reference ensemble (which is in real space percent: PS 0.16,0.26, -; GRS 0.27,0.60,1.0. The values equivalent to a simple ensemble given by two-body for the GRS were averaged over 25 independent configinteraction JJ s(F) V(i_-7)s(7) dFd7) is specified by its urations.) Our results were in good agreement with the spectral distribution v(k) which is in turn fully determined expressions in ref 6 and resolved the discrepancy between by the specific form of the curvature energy. For any the results of refs 5 and 6 in favor of ref 6. one-dimensi2naldistribution, here the distribution of an The total surface area conservation constraint is widely amplitude s(k), Gaussian form has the maximal e n t r ~ p y . ~ accepted to be a crucial physical requirement' and is Hence Gaussian surfaces are maximally random and can enforced in a natural way with the help of stretching be expected to be a good approximation to interfacial energyI2 AE, ' 1 2 K~ S(AS/Si)2,where S and Si are the ensembles with a high degree of thermal disorder. On the actual and initial areas and A S is their difference. AE, other hand the lamellar state is the ground state for stiff enters the Boltzmann factor after each MC step. A large membranes (for - 2 ~< ii < 0). The situation between those value of the stretch modulus KS keeps the fluctuation of boundary scenarios,where the balance of curvature energy, the surface area negligibly small (Figure la), i.e. K~ is used entropy, and the topological "chemical" potential ii under here as a simple implementation of the surface area the constraint of interfacial area conservation determines constraint. Simulation of stretchable surfaces with rethe structure, is not well understood. In particular the alistic K~ is possible but would break the scale-invariance question of the influence of the saddle-splay rigidity R is of the system which we want to avoid in the context of this a point of great interest. The formalism which is based paper. Without constraint the interfacial system is seen on Gaussian states cannothelp in resolving these questions. to gradually drift into trivial states with different surface A general and straightforward method-based only on area. This behavior is a useful check and should be the Hamiltonian %-is required to investigatethe stability characteristic for systems based on an interfacial Hamilof interfacial ensembles for given values of K and 8. tonian when the interfacial area is uncontrolled. The outlined MC algorithm is of very general nature 2. Methods and its more extensive application is only restricted by its large CPU demands. Obvious candidates for simulation Simulation of Continuous Surfaces. The interfaces are all Hamiltonians which comprise surface invariants. used by our MC simulation are given in the implicit In surfactant science, cases of interest are, among others, representation Hamiltonians which include spontaneous curvature, Hamiltonians which contain higher order bending terms (P, s(i) = 0 (3) F K , ...), or systems away from bulk symmetry. where s(F) is expressed as an expansion in an appropriate Application to Disordered Surfactant Interfaces base function space. A MC step is effected by a small We describe disordered surfactant interfaces as a Fourier change in the amplitudes of the base functions. After expansion of the implicit surface representation s(F). The each change of the interface the new state s' is again given expansion should be limited to the physical range of in analytical form. The curvature energy eq 1is evaluated k-values, up to some cut-off value .k, The number of on N surface points 2i and the step is accepted according degrees of freedom in the finite system under investigation to the usual criteria in the Metropolis scheme. The mean is then the volume of the corresponding phase space. and the saddle-splay curvatures are given bylo As short-wavelength fluctuations of the surface do not contribute to topological changes, they are usually taken
ss
(9) Shannon, C. E.;Weaver, W. The Mathematical Theory of Communication; University of Illinois Press: Urbana, IL,1949. (10) Weatherburn, C. E. Differential Geometry of Three Dimensions; Cambridge University Press: Cambridge, 1965.
(11)Gilat, G.; Raubenheimer, L. J. Phys. Rev. 1966, 144, 390. Raubenheimer, L. J.; Gilat, G. Phys. Rev. 1967,167,686. (12) Helfrich, W. 2. Naturforsch. 1973, BC, 693.
Letters loZ
1.0
Langmuir, Vol. 10, No. 2, 1994 347
I
I
1
.................................................................-.............
c
r
A
0.8
Ll
0.6 T=O.O
0.4
L
1.2 1.0
0
MC step Figure 1. (a) Relaxation of the internal energy (in units of the initial energies) of Gaussian random interfaces for K = 1, for different values of the ratio 7 = +/K. For 7 = 2 there is no equilibration even after an additional 2.0 X 106 MC steps, but a steady drift with nearly constant slope. Dashed lines show unstable runs for 7 = -0.1 (left dashed line) and for 7 = 3.0 (right dashed line). The dotted line showsthe constancyof the surface to volume ratio S/V during the run with 7 = 0.4. (b) Relaxation of Gaussian curvatures for several values of 7 in units of the initial Gaussian curvatures.
into account by the renormalization of the elastic moduli K and Z.13 In the resulting scale-invariant description,14 the upper limit of the k-range is proportional to the characteristic wave vector KO. The phase space volume is greatly reduced and the simplified, scale-invariant configurations are independent of the dilution. Use of the renormalized theory allows us to focus on the influence of the elastic moduli on the topological structure. However, in this simplified picture it is not possible to study the effect of short-wavelength fluctuation on the elastic constants. Moreover, the symmetric-asymmetric transition in L3 phases cannot be studied in a scale-invariant ensemble. We chose the size of the simulation box as seven units of the structural length scale = 27r/ko, the grid size 64 X 64 X 64, with periodic boundary conditions. Finite size effects were checked by rerunning samples of twice the side-lengthand a grid density of 128,without visible effect on the results. Finite-size effects were only observed for side-lengths smaller than 4x0. If the upper limit in Cz-space is taken as a small multiple of ko,the number of degrees of freedom is of the order of 5 X lo6and the simulation becomes untractable. We were therefore forced to study systems where all k-values fall into a limited range Ak around the characteristic wave vector KO. From our earlier work7v8we know that a narrow spectral range is physically realistic for high values of the surface elastic modulus.
Figure 2. Typical initial Gaussian sponge configuration used for the MC simulation. The sample volume used has a side length of 7x0, and Gaussian spectral distribution a 6(k - ko). For this configurationthe initial surfaceparameters are (S/v>i= 2/(d3?r) 0.367, ( K ) i = -'/& (H2)i = '/m in units of ke6 Relative errors in the measurement of the surface parameters for a 64 X 64 X 64 grid are = 1-2 % for S/ Vand ( K )and a systematicoverestimate for ( W )of =lo%.
All runs were started with a superposition of random waves with a characteristic wave vector ko,corresponding to Berk's model of disordered ~tructure.~ After a typical run of some 2 x lo5 MC steps, using 4 X lo4 degrees of freedom (3 X 104wave vector components and lo4phases) the internal energy reached an equilibrium state (Figure la) with the width Ak = O.lk0. With a side-length 7x0,the number of degrees of freedom of the system should be N = 4.1rko2 Ak-V = lo5 and the states are reasonably well described by our basis set. With a larger basis and much longer running times, one should expect to arrive at states with a broader range of 12-values. In runs close to the local stability boundary, ii = - 2 ~ of , the disordered phase, the observed first stage of the equilibration, where one finds a relatively rapid configurational relaxation, was followed by a steady downward drift, possibly toward a phase transition. Unfortunately. it takes a much larger computational effort to study this process, and this is not done in the present work. The algorithm is most efficiently implemented on parallel computers. We used a Connection Machine I1 with 8192 processors. 106 steps used about 20 h of CPU, making a large number of runs untractable. However, on a more recent model, CM V, the performance is already improved by a factor of 2. 3. Results Monte Carlo runs were started with Gaussian random surfaces with a spectral function v(k) = 6(k - ko), Figure 2. The final equilibrium structure depends on both elastic moduli, K and ii. It is convenient to discuss the results separately in terms of their dependence on the absolute values of K and ii and on the temperature independent, dimensionless topological ratio which we denote i?
(13) GoluboviC, L.; Lubensky, T. C. Phys. Rev. B 1989,39,12110,and references therein. (14) Porte, G.; Delsanti, M.; Billard, I.; et al. J. Phys. ZZ 1990, I , 1101.
7=-K
As outlined above, it is natural to assume that Gaussian
348 Langmuir, VoZ. 10, No.2, 1994
. 0.6
%*-
2.2
0.5
1.0
Letters
0.2 0.0 -0.2
-0.4
-0.6 -0.5
0.0
1.5
! 2.5
2.0
z
Figure3. Regionsof differenttopologicalexcitationsdependent on the topological parameter T = - i i / ~ . The curves depict change in the internal energies and curvatures upon equilibration, expressed relative to the values in the initial state: ( 0 ) internal energy; (A) mean square curvature; ( 0 )Gaussiancurvature. The curvature Hamiltonian is unstable in the shaded regions. The verticle dotted lines mark the approximate positions of the characteristic topological ratios ~~1and 7 ~ 2 . We use the labels (ii), (iii!, v d (iv) from the text to distinguish topological tendencies in the respective regions of 7.
states are stable for “low” values of the bending rigidity. A typical experimentalvalue for the stiffness of disordered systems is K of the order of 1. Even though the equilibration runsfor K = 1 in Figure 1 show different relaxation behavior, the distribution of amplitudes of the final states is never far removed from the initial Gaussian. The deviation is only visible for the boundary cases r = 0.0 and r = 2.0. This indicates that Gaussian states are indeed good model ensembles for realistic interfacial sponge states, as long as the rigidity K is small. A quantitative measure of the suitability of random Gaussian state description is obtained from the change in average internal energy upon MC equilibration at different values of K. For a sequence of values K = l/3,1,2, and 3, we found at 7 = 0.3 that the was =14, final value of the mean square curvature (P)f 23,41, and 55% below its initial value. In all those cases the final topology is approximately the same; the final values (K)f are 2-3% above the initial value. In Figure 1 we saw that equilibrationsof internal energies and saddle-splay curvatures differ markedly for different values of r, We can study the r-dependence in more detail in Figure 3 where the curvatures of the final states (P)f and (K)f are compared to their initial values, indicating characteristic topological changes. The 7-dependencecan be classified as follows: (i) T < 0 and T > 2. No stable runs were found in this region. This is easily understood from the well-known local stability criterion for the bending energy - 2 ~< f < 0 and is a good check of the numerical algorithm. However, the stability limits may change when higher order curvature terms are included. Because we do not consider higher order terms which would lead to physically reasonable states of small vesicles (but would break the scaleinvariance),we did not follow the evolution in the unstable regime. (ii) 0 5 T S 0.3. The final state has a lower mean square curvature (P)fthan the initial state and higher mean
saddle-splay curvature (K)f. From the Gauss-Bonnet theorem J clc2 dS = 47r(n, - nh) = 27r XE (where n, and nh are the numbers of components and handles and XE is the Euler characteristic), we can conclude that equilibrium structures in this region are rich in tubular handles and low in mean square curvature energy. This is reminiscent of periodic minimal surfaces, although fiial structures here and in all other regimes are isotropic. Resemblance to cubic phases can be conveniently measured15-17 using the Euler characteristic suitably normalized according to the scaling properties of the systems. Like the free energy, the Euler characteristic of a surfactant system is an extensive quantity invariant to the change of scale. Therefore, under dilations the invariant scaled Euler characteristic is XEO V2/A3, where XEO is the Euler characteristic per unit building block. In the notation of ref 17, the areaper building block A0 can be written in the form A0 ALo2 where LO is the cell dimension and the constant A depends on the particular topological structure. Hecce the scale-invariant Euler characteristic is R E O = XE’/A~and is related to the saddlesplay curvature by j i ~ ’ = (K)/(27r(S/W2). For periodic minimal surfaces170.28 < bo < 0.37. For the isotropic disordered phases, we calculated as a = 0.24, function of the topological ratio (at K = 1) with bo 0.2, and 0.1 for r = 0, 0.4, and 2.0, respectively. An exceptionally long run (6X 105 MC steps) was performed for the extreme case K = 50. There (P)f/(P)i = 0.1 with = 0.27. It is worthwhile mentioning that high final the existence of nonperiodic, isotropic minimal surfaces is an unresolved problem.18 For K = 1.0 the difference between the initial and final state which is shown in the calculated average curvatures is not visually discernible in three-dimensional real space images of the interfaces. In order to illustrate the equilibration process in visual images, we enhance structural changes by lowering the temperatures at constant r. In Figure 4 we show a state which evolved for a stiff system r = 0.0, K = 10.0. In comparison with Figure 2 the interfacial shape is visibly smoother and has an increased number of handles. (iii) 0.3 S T 5 1.0. Here the equilibrated states are closest to random bicontinuous states (Figure 3). In particular, for 0.3 S r S 0.4 the topologies of the initial and final states are nearly unchanged. For 0.4 S 7 S 1.0 the final states are characterized by a decrease in both saddle-splay curvature and mean square curvature. The structure is hence less connected (and not evolvingtoward a minimal surface) and has at the same time flattened. This is consistent with the picture of connected, molten lamellae. We present a real-space structure for 7 = 0.4 in Figure 5 which shows an emerging locally lamellar structure. (iv) 1.0 S 7 I2.0. In this regime the final state of the system has less connectivity and higher ( H 2 ) f than the initial state. Our result supports the idea of an equilibrium structure which is on average less connected and more parabolic, Le. connected “wormlike” shapedg or clusters of deformed closed objects. Figure 6 shows a state for the case of r = 2.0. The loss of connectivity is obvious. The (15) Hyde, S. T. Colloq. Phys. 1990, C7,209.
(16) Anderson, D.; Wennerstriim, H.; Olsaon, U. J. Phys. Chem. 1989, 50, 1335. (17) Gompper, G.; Kraus, M. Phys. Reu. E, in p m , and referenma therein. The authorsdefine an order parameter field $(?) on lattice sitee and update those. Their simulationuses a Ginzburg-Landau free energy, Le. not a curvature Hamiltonian. (18) Hoffman, D. A. Colloq. Phys. 1990, C7,197. (19) Safran,S. A.; Turkevich,L. A.; Pincus, P. J. Phys.,Lett. 1984,45, 69.
Langmuir, VoZ. 10, No.2, 1994 349
Letters
Figure 4. Three-dimensionalimage of the interfacial structure, K = 10.0,7= 0 (regime ii) after 3 X 105MC updates. It has ratios of final and initial curvatures of (K)f/(K)i* 1.25 and (W)f/ (H2)i 0.25 and a final value of = 0.26.
-
'1
Figure 6. Three-dimensional image of an interfacial structure, K = 1.0,7= 2 (regimeiv) after 6.5 X 106MC updates. It has ratios of final and initial curvatures of (K)f/(K)i* 0.39 and (W)d (W)i* 1.33.
relevant. It should be mentioned that we have done runs with broader, variationally determined8 spectra to make sure that the result is not an artifact caused by the localized spectral distribution. Although we are far short of equilibrating broadened spectra, we were able to confirm that the boundaries rcl and rc2 (Figure 3) are still valid, with a deviation of a few percent. We can therefore distinguish two characteristic values of the topological parameter r where the change in saddlesplay and mean-square curvature equals zero rcl= 0.3-0.4
rl
J
Figure 5. Three-dimensional image of the interfacial structure, K = 5.0,7 = 0.4 (regime iii) after 9 X 105 MC updates. It has ratios of final and initial curvatures of (K)f/(K)i * 0.66 and (W)r/ (W)i 0.43.
tendency to form lamellar structure is here supressed by an unfavorable contribution from the saddle-splayenergy. Nevertheless in this regime-as in any other regimelamellar states would eventually form for K >> 1 when entropic contributions are negligible. It is now interestingto ask how universalthe dependence of the geometry on the topological parameter r is. This has been done by successively cooling the samples from K = 1to K = 5,lO. Although our computational power is barely sufficient to overcome the long relaxation times at K = 5, we won for all cases conclusive evidence that the general nature of the change in topology from initial to finalstateis unchanged. This indicatesthat the parameter r determines universally the topology throughout the examined range of K values which we deem physically
and
rcz= 1.0
(6)
The second value ~~2 has been interpreted by Porte et ~ 1 on the grounds of a local energy argument as the border value between a bicontinuous structure and vesicular shapes. Although we cannot confirm an instability toward vesicles as proposed by Porte et al.,we can conclude from our results that the structure at rc2 turns into more parabolic shapes. The first characteristic value at rcl = 0.3-0.4 is new and separates a region in which disordered minimal surfacelike or disordered lamellar structures are preferred. This is not self-evident from purely energetic argumentsas there is an energetic penalty for connectivity for any r > 0. It can only be explained on entropic grounds. rcl is a clear manifestation of a topological contribution to the entropy of the interfacial system. 4. Conclusion
In this paper we have introduced a new Monte Carlo technique3 which can simulate the configurational evolution of a large class of interfacial geometries. (Previous surface MC simulations known to the authors are dealing either with different physical systems, such as tethered and polymerized surfaces or with similar systems in other than curvature-elastic theoretical de~cripti0n.l~)It is based on a genuinely continuous surface representation, which allows accurateevaluationof surfacecurvature terms (20) Porte, G.;Appell, J.; Bassereau, P.;et al. J. Phys. (Park) 1989, 50, 1335.
.
~
350 Langmuir, Vol. 10,No. 2, 1994
used to describe the Hamiltonian of surfactant film systems. Thus configurational equilibration of the film structure is directly feasible under the physically crucial constraint of film area conservation. The only input parameters, the bending and the saddle-spay moduli, have a clear physical meaning. Moreover the method is easily extendable to more complex Hamiltonians, eg. including higher order curvature terms and spontaneous curvature. We applied the technique to the Helfrich Hamiltonian to study the influence of the bending moduli K and on random sponge geometries and our results are in agreement with several earlier works16120*21 which emphasized the importance of the saddle splay modulus R in determining the topological nature of the disordered surfactant phases: an entropic contribution to the free energy, coming from the fluctuations in the local topology, favors departures from the lamellar state; hyperbolic tubular (21) Wennerstram, H.; Ohon, U. Langmuir 1993,9,365.
Letters contacts become energetically less unfavorable for higher values of E. On the basis of MC simulations, we can add two quantitative remarks to the above discussion. The topological characteristics predominantly depend on the ratio of RIK, rather than the separate values of the two elastic moduli. As the magnitude of this ratio goes through 1 and ~~2 = 1the nature of two specific values ~ ~ 0.34.4 thermal excitations of the surfaces changes respectively from predominantly hyperbolic, highly connected structures via flat connected sheeta to loosely connected parabolic domains. Acknowledgment. The authors are grateful to M. Robbins, I. S. Barnes, S. T. Hyde, M. Teubner, and S. A. Safran for useful discussions and G. Gompper for the preprint of ref 17. P.P. was in part supported by a grant of the von Hoesslinschen Foundation of the City of Augsburg, Germany.
