Kirstin Purdy
Advisors S.D. Mahanti, Aniket Bhattacharya
Department of Physics, Michigan State University
Research Experience For Undergraduates
August 14, 1998
In the past there have been theories about phase transitions in surfactant self-assembly. Phase transitions are defined usually by a peak in the specific heat of a substance, but a theory by Blankschtein, Thurston and Benedek(1) also states that in surfactant self-assembly a phase transition can also be identified through the second moment of the cluster distribution function, which is related to the fluctuations in the average size of the micelles formed in self-assembly. A simulation using a simple two dimensional lattice model and a Monte Carlo method of moving the surfactants was used to investigate this theory. The Monte Carlo movement consisted of a snake like reptation of the surfactant chains. Beginning with the results published by C.M. Care(2), which we believe exhibit phase transition properties, the simulation was initially carried out to reproduce his results for different strengths of head-solvent interactions. The properties of the surfactant self-assembly were explored in search of a phase transition. A maximum in both the specific heat and the average cluster size fluctuationswas observed . These peaks suggest that a phase transition is occurring, but the specific heat peak happens regularly at a temperature higher than that of the average cluster size fluctuations. The location of the peak in the cluster size fluctuation also seems to be dependent on the strength of the head-solvent interaction. Because of the instability of the average cluster size fluctuations the location of the phase transition is still uncertain.
Introduction
Surfactant molecules are by nature bipolar. They usually contain a polar head group which, in water, is hydrophyllic, and a hydrophobic hydrocarbon chain tail. When placed in water these amphiphiles will form clusters of molecules known as micelles, which allow the molecules to exist an equilibrium state(3). Self-assembling surfactant systems have many applications in industry, and are also essential parts of biological systems. Lipid bilayers which make up cell membranes are self-assembling systems that are the primary control for transport of materials into the cell. In industry, self-assembling surfactants can be used in the design and fabrication of material composites, because of the intrinsic order that the self-assembling properties have (4).The purpose of this paper is to show the behavior of a minimal model self- assembly simulation.
The behavior being studied is the possibility of a phase transition in the solvent-
surfactant system. There are many theories which recommend a way in which one can detect a
phase transition in surfactant self-assembly, but the theory we are looking at, by
Blankschtein, Thurston and Benedek, relates to the moments of the cluster distribution
function. Blankschtein et al. have theorized that the critical point and coexistence curve
of the surfactant solution can be determined by low order moments of the cluster
distribution function. The moments of order a have been defined as
Ma(X,T,p)=En^aXn (1)
where Xn is the concentration of the clusters. In our system M0 is the number of clusters
in the system, and M1 is the number of chains. Blankschtein et al. also developed a
recursive relationship
M(a+1)=M2 (dMa/dX) (2)
which shows that the second moment (M2) provides an essential link to the physics of the
system. In C.M. Care's paper, a quantity called the average cluster size fluctuation was
presented which is defined as:
sigma^2=Ei^2ni/Eni -
where n=Eini/Eni =M1/M0 is the average cluster size and ni is the number of clusters of
size i. This s is a measure of the polydispersity of the solvent-surfactant mixture.
Because of the relation between sigma and M2 we believe that the peaks shown in Care's
data represent a phase transition. Our simulations were done using the same techniques as
Care's simulations to determine if we could replicate his results and to determine whether
or not there is a phase transition present in micelle formation. Because we believe that
there is a phase transition, we are also interested in locating where it occurs in our
simulation and whether or not that location is comparable to where we believe the phase
transition is in Care's results.
The model we used to simulate the surfactant self-assembly is discussed in the Method section along with the Monte Carlo simulation technique. The results of the simulation and their comparison to Care's results are presented in the Results and Discussion section, and the conclusions thus far about phase transitions seen in our model are presented in the Conclusion section.
The model used to simulate self-assembly, previously used by Care, consists of a two dimensional non compressible lattice of an amphiphile and solvent mixture. This model was done on a lattice because off lattice calculations are much more complex (5). The lattice contained periodic boundary conditions to help minimize finite size effects. The surfactants were represented by chain of s lattice sites with the first i representing the head and the rest (s-i) representing the solvophobic tail. The self-assembly of the surfactants was studied on two different lattice sizes, 64x64 and 128x128, each containing 128 and 512 amphiphiles of length s=3, respectively, for an equal concentration of 9.4% on each lattice. The results from the 64x64 lattice simulation were less smooth than those from 128x128 lattice simulation, making them more difficult to interpret. Yet, the form of the data remained the same. Because the results from the larger lattice are smoother, only those results are shown in the results section. Some simulations were also done for amphiphiles of length s=6, and the results were also similar in form to the results for length s=3, but because the molecules were much longer, it was much easier for the molecules to get stuck in metastable states because of the limited movement they had. Molecules larger than s=3 will be investigated in more detail at a later time.
In this model only nearest neighbor interactions were used. The energy at any given
lattice site can then be written as
U=nHS*eHS +nTS*eTS +nHH*eHH +eBENDING (4)
In this equation, HS represents the head-solvent interactions, TS represents the
tail-solvent interactions and HH represents the head-head interactions, n represents the
number of nearest neighbor interactions of binding energy e. The last term is the bending
energy, which, in this model, is set to zero allowing the molecules to be completely
flexible. All other interactions between parts of the molecules are not introduced in this
model because the above interactions are the minimum needed for micelle formation at low
temperatures. As in ref.(2) the following notation will be used throughout
b=eTS/kT g=eHS/eTS (5)
In this model the head-head interactions are set to zero (eHH=0) and the tail-solvent
interactions are assumed to be solvophobic (eTS>0). Simulations were done for a variety
of different head- solvent interaction strengths. If both head-solvent and tail-solvent
interactions were solvophobic there would be complete phase separation at low
temperatures. As g is decreased (eHS<0) the head-solvent interactions become stronger, allowing micelles formation. If the head-solvent interactions are very strong and solvophyllic the strength of the interaction will be strong enough by itself that the repulsion from the tail-solvent interaction will not be enough to warrant micelle formation. In this case most of the molecules will remain monomers. At high temperatures the molecules are mostly monomers because the energy is high enough that the hydrophobic tail- solvent interactions do not warrant micelle formation.
To iterate the model to an equilibrium state a Monte Carlo method was used. For most of the data presented 1x10^6 Monte Carlo steps were made at each temperature each with the number of attempted moves equal to the number of amphiphiles on the lattice. For thermalization, 5x10^4 additional Monte Carlo steps were used. The surfactant molecules move by a snake like reptation across the lattice. Although, most data was taken by beginning with an initial high temperature random configuration, some data was also taken by using a low temperature configuration and heating it. There was no significant hysteresis.
At each temperature the total energy of the system and the average energy of the system
is calculated and from that the energy fluctuation is calculated along with the average
cluster size and its fluctuation. From the fluctuations of the energy
sigma = sqrt((E^2)-(E)^2) (6)
the specific heat can be calculated by sigma^2/T^2. The results from calculating the
specific heat by the fluctuations on the energy were also compared to the specific heat
calculated by taking the derivative of the energy; they were found to agree with one
another.
The average cluster size and its fluctuations were also calculated because of their relation to the second moment of the cluster distribution function. The equations were shown previously in the introduction section. Average cluster size distributions were also calculated to investigate the properties of the simulated system. In this simulation, the size of a cluster is defined by the number of molecules that have tails which are nearest neighbors. When the average cluster size distributions were calculated, the model was iterated another 1.5x105 Monte Carlo steps and averaged over both the last 25k and 50k Monte Carlo steps. The difference between averages was small as expected, because the average distributions should vary less for the higher Monte Carlo iterations, evening out the fact that there are less values to be iterated over (25k versus 50k).
The resulting curves for the average cluster size fluctuations for the
128x128 lattice are shown in Figure 1.The results shown are for four different values of
g, g=-0.2, g=-0.4, g=-0.6, g=- 1.0. The multiple runs for g=-0.6 and g=-1.0 are shown to
demonstrate the repeatability of our results. For each negative value of g there is a peak
in the results. For values of gł0 the average cluster size fluctuations diverged. Both
the simulated peaks and Care's peaks tend to shift towards lower temperatures as g is
decreased. Our results agreed well with Care's for g=-1.0 and partially agreed with Care's
results for g=-0.6. The location of these peaks is important because we believe that the
changing of the polydispersity represents a phase transition. The difference between our
results and Care's results for high values of negative g is mostly a result of the
instability of the fluctuations at the low temperature. Because only slight movements in
the molecules can can use drastic changes in the polydispersity, it is very easy for the
surfactants to get stuck in a metastable phase. The peak in the polydispersity of the
micelles implies that something else is also happening with the energy of the system,
since the energy of the system is directly related to the formation of micelles.

Figure 1(above). Graph of Average Cluster Size Distribution. The peaks are clearly visible
for the higher values of g but the peak for g=-1.0 is very small.
Figure 2 (below). Graph of Specific Heat. Notice how the peak gets sharper as g increases.
To further investigate this phase transition, the specific heat of the system was calculated. A graph of the specific heat as a function of temperature is shown in Figure 3. The specific heat results show a peak as well for each value of g>=0, which narrows as g decreased. The peak in the specific heat definitely suggests a phase transition. Even more interesting is that, althought the peaks have varying g, they are all in approximately the same location. However, the specific heat peaks are consistently located at a temperature higher than the peaks of the cluster size fluctuations. Because of this difference, and the difference between our results and Care's results we are continuing to investigate the surfactant system, in order to more precisely locate the phase transition temperature.
Another important observation to make is the double peak in the average cluster size fluctuations for g=-0.2. Upon closer inspection it may be possible that there are double peaks in all of the curves of the fluctuations. This would imply that there may be two transitions, one stronger than the other, since when looking at the specific heat curves, there only appears to be one peak.
Finally, the average cluster distributions about the peaks in the
fluctuation were investigated. Graphs of the cluster distributions at g=-0.6 and g=-1.0
about the peaks are shown in figures 3 and 4 respectively. The graphs simply tell us
visually about the degree of the polydispersity. For g=-0.6, the polydispersity is much
higher than for g=-1.0. About the peak for ^M g=-1.0, it is easy to watch the
polydispersity level off and then begin to decr ease again, because there is a plateau
that forms in the average cluster distribution.
Figure3. Average Cluster Distribution for g=-0.6

Figure 4. Average Cluster Distribution for g=-1.0
Theses cluster distributions were averaged over 2.5x104 Monte Carlo steps after an
additional 1.25x105 Monte Carlo steps were taken after the initial run.
To investigate the phase transition of a self-assembling surfactant system a simple Monte Carlo model was used. Both the specific heat as a function of temperature and the change in the fluctuations of the average cluster size as a function of temperature were calculated. The specific heat was investigated, because it should exhibit a sharp peak at a phase transition. The average cluster size fluctuations were investigated because they are related to the second moment of the cluster distribution function, which Blankschtein et al. believe relate the self-associative nature of the system to a phase transition.
A maximum in both the specific heat and the average cluster size fluctuations was
observed, suggesting that a phase transition is occurring. The location of the peaks
should help us to determine where the phase transition is. Currently the results of the
simulation show very similar specific heat peak location for different strengths of
head-surfactant interactions. The peaks also narrow as g becomes less negative. But, the
location of the average cluster size fluctuation is consistently lower than the specific
heat peak. Our results for the average cluster size fluctuation agree with Care's results
for g=-1.0 and agrees partially with Care's results for g=-0.6. In Care's results the peak
of the average cluster size fluctuations shifted slightly to lower temperatures with
increasing g, but in our simulation our result shifted much more towards lower
temperatures and away from the specific heat peak. This difference in results could be
from an intrinsic instability of the average cluster size fluctuations, as well as an
increased tendency towards metastable states at low temperatures for high values of g. The
difference between Care's results and our results is still being investigated. ue to the
difference in locations of the peaks in the specific heat and the average cluster size
fluctuation, the location of the phase transition is still unknown.
(1) Blankschtein, Daniel, George M. Thurston, and George B Benedek. Theory of Phase
Separation in Micellar Solutions. Phys. Rev. Letts. 54 955 (1985).; Blankschtein, Daniel,
George
M. Thurston, and George B Benedek. Phenomenological theory of equilibrium thermodynamic
properties and phase separation of micellar solutions. J. Chem Phys. 85 4558 (1986).
(2) Care, Christopher M. Cluster Size Distribution in a Monte Carlo Simulation of the
Micellar
Phase of an Amphiphile and Solvent Mixture. J. Chem. Soc. Faraday Trans. I, 83 2905
(1987).
(3) Israelachvili, Jacob N. Physical Principles of Surfactant Self-Association into
Micelles,
Bilayers, Vesicles, and Microemulsion Droplets.
(4) See articles in special issue of Science 277, Frontiers in Material Science, Aug 29
(1997).
(5) Bhattacharya, Aniket, Mahanti, S.D., and Amitabha Chakrabarti.Self-Assembly of Neutral
and Ionic Surfactants: An Off-Lattice Monte Carlo Approach. J. Chem. Soc. June (1998).