Micellar Crystals in Solution from Molecular Dynamics Simulations
Abstract
Polymers with both soluble and insoluble blocks typically selfassemble into micelles, aggregates of a finite number of polymers where the soluble blocks shield the insoluble ones from contact with the solvent. Upon increasing concentration, these micelles often form gels that exhibit crystalline order in many systems. In this paper, we present a study of both the dynamics and the equilibrium properties of micellar crystals of triblock polymers using molecular dynamics simulations. Our results show that equilibration of single micelle degrees of freedom and crystal formation occurs by polymer transfer between micelles, a process that is described by transition state theory. Near the disorder (or melting) transition, bcc lattices are favored for all triblocks studied. Lattices with fcc ordering are also found, but only at lower kinetic temperatures and for triblocks with short hydrophilic blocks. Our results lead to a number of theoretical considerations and suggest a range of implications to experimental systems with a particular emphasis on Pluronic polymers.
I Introduction
Micelles of multiblock polymers are finite aggregates, typically around fifty polymers or less, where the insoluble blocks shield the soluble ones from contact with the surrounding solvent. Depending on control variables (temperature, polymer concentration, pH, etc..) micelles may selfassemble into gels that exhibit longrange order such as bcc, fcc, hcp or other, more unusual, crystals. Micellar crystals exhibit a number of unique properties that have made them extremely attractive for fundamental studies as well as for applications (1); (2). Extensively studied experimental systems include aqueous solutions of Pluronics (also known as Poloxamers), ABA triblocks where A is Polyethylene oxide (PEO) and B is Polypropylene oxide (PPO) (3); (4); (5) and inverted Pluronics, where the A blocks are PPO and the central block B is PEO (6); (5) as well as nonaqueous systems such as PolystyrenePolyisoprene (PSPI) diblocks in decane (7) and other solvents (9); (10); (11); (12).
Micelles in solution are highly dynamical entities with polymers continually being absorbed and released through time. Therefore, a micellar crystal has a considerably intricate structure, where the long range order remains stable as the individual polymers are constantly hopping from one micelle to the next. Theoretical approaches such as density functional or mean field theory (8); (13); (16); (14); (15); (2); (17) directly study the ordered micelles and ignore the dynamical degrees of freedom of the polymers. Studies using molecular dynamics (MD) have the advantage of providing a reasonably realistic description of the dynamics, thus allowing the investigation of the role of single polymer degrees of freedom. In contrast with other approaches, MD also offers the important advantage that no assumptions need to be made about what is the possible thermodynamic state of the system.
The goal of this study is to predict phase diagrams of triblock polymers using MD simulations and to gain an understanding of the dynamics of micellar crystal formation. Because of our ongoing interest in Pluronic systems in aqueous solutions (18), we examine systems of triblocks, where the beads are hydrophilic and beads are hydrophobic. Although there have been a number of previous studies of multiblock copolymers in solution using MD, see Ref. (19) for a recent review, the prediction of crystalline structures presents substantial difficulties. Experimentally, it is well known that the approach to thermodynamic equilibrium is slow in these systems, with time scales of the order of minutes or hours. Therefore, even with suitably coarsegrained models, the long time scales involved provide a considerable challenge for MD simulation studies.
In this paper, we provide a detailed investigation of self assembled micellar crystals using MD. We aim to understand the mechanism by which they form, predict their range of stability and elucidate their static and dynamic properties. We also present an indepth study on the challenges associated in reaching the thermodynamically stable state using MD and a successful strategy to overcome them.
Ii Model and Simulation Details
ii.1 Simulation Details
Systems of polymers are modeled by coarsegrained beads in an implicit solvent. Ref. (18) provides a more detailed justification than the outline provided here. Individual polymers are symmetric triblocks, where beads are hydrophilic and beads are hydrophobic. All systems in this work are monodisperse with fixed values for and . Nonbonded pair potentials consist of a standard attractive LennardJones potential for hydrophobic interactions
(1) 
and a purely repulsive potential for hydrophilic interactions
(2) 
Pair potentials are cutoff to zero at . All beads have the same mass . , and are arbitrary and uniquely define the units of all numbers in the simulations. Neighboring atoms in the polymer chain are connected together with a simple harmonic potential. The packing fraction of the polymers is given by
(3) 
where is the number of polymers in a box of linear dimensions and is the number of beads in each polymer. Simulation boxes are cubic with periodic boundary conditions. For this work, we use a time step of . All simulations were performed using the LAMMPS software package (20) in the NVT ensemble via NoséHoover dynamics (21).
The temperature sensitivity of Pluronic systems is a reflection of the underlying strong temperature dependence of the hydrophobic effect. To model temperature dependent phases in such systems, the solvent quality for the beads must change (18). In this paper, we keep the solvent quality fixed and vary only the kinetic temperature when needed to equilibrate the systems properly.
Recently, a different implicit solvent model has been developed for describing Pluronic systems (22), where the pair potentials for the coarsegrained simulation are fitted to results obtained from allatom MD and quantum chemistry simulations. Although the coarsegraining is different (each monomer , represents one PEO or PPO monomer) than in our model, the resulting potentials are quite similar to the ones in this work, Equations 1 and 2, with the only significant difference that the values of are different between the two and monomers, and shows a minor maximum.
ii.2 Observables
A number of observable quantities are monitored for every recorded time step during a simulation run. The micelles themselves are identified by the same algorithm used previously in Ref. (18). Any hydrophobic beads within a distance of of one another are identified as belonging to the same micelle. Identified micelles containing less than 3 polymers are typically free polymers in the process of being transferred from one micelle to another and are removed from further consideration. Observables such as statistics of micelle aggregation number, gyration tensor, center of mass, and micelle lifetime are calculated and examined for every simulation performed in this work. Methods used to calculate these are described in Ref. (18).
The structure factor is calculated over the center of mass coordinates of all micelles in the system using the formula
(4) 
with the components of as multiples of due to the use of periodic boundary conditions. The peaks in are then used to reconstruct the full 3D real space lattice basis, if it exists. In this manner, is not being used to simulate real scattering intensities as may be obtained by Xray experiments, but as a mathematical order parameter to discriminate between the different ordered structures that may be present. For convenience, the normalization is chosen so that . A more sophisticated treatment that allows continuum values of , suitable for quantitative comparisons with Xray experiments, has been recently introduced (23), but we do not use it here.
Individual polymers are constantly hopping from one micelle to another. This is quantified over the entire simulation box as an overall rate of polymer transfer, , by examining contiguous simulation snapshots. Sets of indexed polymers belonging to each micelle are compared between the snapshots to find the number of polymers transferred. The rate is then expressed as a fraction of the number of polymers in the box transferred per one million time steps. A polymer that is transferred out and back to the same micelle between snapshots will not be counted by this analysis, so snapshots are recorded every 100,000 time steps to minimize undercounting. At , a typical micelle only loses/gains one polymer per ten snapshots recorded.
Radial distributions of the beads surrounding micelles are also of interest. These are calculated by creating a histogram with bin width and then counting the number of beads belonging to a micelle that fall between and , where is the distance of the bead from the center of mass of the micelle. Good statistics require averaging this histogram over all micelles in the simulation and over all time steps after the micellar crystal has formed. The average histogram is transformed into a radial density distribution of beads around a micelle by calculating the packing fraction of beads in each bin
where refers to either or beads. We use to balance smooth graphs with the need for long simulation runs to obtain detailed statistics.
Iii Micellar crystals studied using MD
There are two major challenges faced in using MD to determine equilibrium phases. First, the simulation must last longer than any of the relaxation times in the system. Second, the simulation box size must be chosen properly to avoid finite size effects. The first problem is related to the kinetic temperature at which the system is run. The second problem can become particularly severe for simulating crystals with three dimensional order, where the incorrect choice of an even large box size can force the system into very distorted ordered phases. These two issues are addressed systematically using simulations of the polymer. It provides a coarsegrained description of one of the most extensively studied Pluronics, F127 (3). The conclusions of this study lead to a general methodology valid for any other polymer.
iii.1 Micelle crystallization and kinetic temperature
Initial simulations are performed with a system size of at a kinetic temperature and with concentrations of , , , and run over million time steps each. In all cases, the initially randomly placed polymers aggregate into micelles quickly in a few thousand time steps.
Visual examination indicates that concentrations at and above are strong candidates for micellar crystals. Micelles are locked in place and only move a few from their average positions. The analysis of the order parameter , however, indicates no longrange ordered structures exist in any of these initial simulations. An inspection of the polymer transfer shows that it remains negligible throughout all simulations. This is confirmed by calculating the mean squared displacement averaged over all beads in the system. Figure 1 clearly shows nondiffusive behavior at . These results suggest that the individual micelles are quickly frozen in a configuration that is not representative of equilibrium, thus preventing the entire system from reaching thermal equilibrium.
The lack of equilibration of micellar degrees of freedom suggests subsequent runs at a larger kinetic temperature . A single simulation run at formed a textbook fcc lattice after about 10 million steps, and is shown in Figure 2. Throughout the duration of the run, the rate of polymer transfer was substantial, about of the polymers in the box every million time steps, even after equilibrium is reached. Playing a movie of the simulation shows that after the lattice formed, micelles do not appear to move except by vibrating about their average positions. However, while the micelles appear static, polymers are constantly being exchanged among the micelles, so that any individual polymer will eventually explore the entire simulation box. This is independently confirmed by the analysis of the mean squared displacement of beads in Figure 1, which shows a classic diffusion result at .
The behavior of the number of micelles in the box as a function of simulation time is of particular interest. In simulation runs where micellar crystals are found, the system reaches a plateau where remains constant, see Figure 3 for an example. Despite their dynamic character even at equilibrium, then remains constant for the duration of the run. This correlation is a general feature in all simulation runs performed. Every single one that leads to a stable plateau in as a function of time formed a micellar crystal confirmed by peaks in the order parameter .
No. configs  Ordering  
1  500  32  textbook fcc  15.625 
3  500  unstable  none  
1  600  36  none  16.67 
3  600  36  distorted bcc  16.67 
5  800  48  distorted bcc  16.67 
5  900  54  textbook bcc  16.67 
1  900  55  distorted bcc  16.36 
2  1000  unstable  none  
4  1000  60  distorted bcc  16.67 
No. configs  Targeted  Ordering  

1  16 micelle bcc  267  16  textbook bcc 
4  108 micelle fcc  1800  unstable  none 
4  128 micelle bcc  2134  128  textbook bcc 
3  250 micelle bcc  4168  unstable  none 
Box size effects
The influence of the simulation box size choice is assessed by running additional simulations with , , , , and with fixed concentration and . Several different random initial configurations are used at each system size to ensure repeatability. Table 1 summarizes the results of all these simulation runs.
Interestingly, the additional three simulation runs at do not exhibit fcc structures. Instead, each of them remains unstable with never achieving a constant plateau and is devoid of peaks. Larger system sizes do reach stable ordered structures with constant and numerous peaks in . Most of the structures that occur appear bcc when examined visually, but a detailed analysis of the lattices indicates that many of them are distorted. One obviously distorted structure is depicted in Figure 4 where the lattice is body centered tetragonal with .
Other distortions include body centered tetragonal with for and a lattice that appears to be almost exactly bcc when , except that the central micelle in the unit cell is shifted slightly from the true center. Lastly, a textbook bcc lattice is formed in the simulation runs with .
Examining the average aggregation number leads to a very illuminating result. It is found that for almost all simulations these numbers are identical. This indicates that the average micelle aggregation number is independent of the box length (at a fixed ) and is only a function of the polymer structure, kinetic temperature and concentration.
Assuming without proof that either the fcc or the bcc lattice represents the real thermodynamic equilibrium of the system, the previous observation then suggests a way to generate magic numbers of polymers to simulate equilibrium states free of finite size effects. In a bcc lattice, each cubic unit cell contains micelles, while the fcc lattice contains . Therefore, in order to obtain a bcc or fcc lattice with by by unit cells in a cubic simulation box, the number of polymers needed to achieve this is given by
(5) 
iii.2 MD Simulations without finite size effects
An algorithm to simulate micellar crystals without finite size effects follows very naturally from the previous results, and is summarized in the following steps.

Concentration selection: The concentration must be chosen high enough so that micelles pack into a potential crystalline state.

Temperature selection: The temperature should be chosen large enough to ensure a significant rate of polymer transfer, but low enough that the micellar crystals are not in a disordered phase. As a rule of thumb, we have been using a polymer transfer around of the polymers in the box every million steps.

System size selection: Calculate the average micelle number via test simulations and use (5) to determine the final system sizes to perform simulations on.

Ensure reproducibility: The formation of micellar crystals is a stochastic process so several simulations with different initial configurations must be run.
The advantage of this algorithm is that steps 14 can be accomplished with relatively modest computer resources on small system sizes, leaving the production runs with large polymer numbers as the only computationally intensive calculations.
iii.3 Micellar crystals of general triblocks
Further simulation runs are carried out on the polymer system to test the algorithm. These additional simulations target 16, 128 and 250 micelle bcc configurations, along with 32 and 108 micelle fcc ones. The simulations are summarized in Table 2. None of the simulation runs targeting fcc phases ever reach a stable number of micelles and correspondingly, indicates there is no longrange order. On the other hand, both the 16 and 128 micelle bcc configurations are perfect and completely reproducible with the expected ordering confirmed by the structure factor. Figure 5 shows a snapshot of the 128 micelle bcc phase after it forms. Further discussion on the calculated structure factors are presented below. Larger simulations attempting the formation of a 250 micelle bcc crystal never resulted in a stable plateau.
The polymer is also studied as it provides a closer representation to the real Pluronic F127, since there are twice as many monomers in a polymer and the level of coarsegraining is less, as discussed in Ref. (14). Applying the algorithm developed above to search for micellar crystals, the concentration was selected at and the kinetic temperature at . Initial simulation runs establish that . The results, summarized in Table 3, are similar to those for the smaller polymer, with the bcc structure being the most commonly found. The fcc phase again only appears in a single simulation run and is not reproducible.
Given the prevalence of bcc lattices so far, a polymer with shorter blocks (), expected to form more crewcut micelles and hence more prone to assemble into a fcc lattice, is also investigated. Applying the algorithm developed above, the concentration is selected at , the kinetic temperature at , and the average aggregation number is found to be . Again, the bcc micellar crystal is most commonly found, as shown in the results in Table 4. No fcc phases are found at this kinetic temperature.
However, the rate of polymer transfer at is larger than per million steps, so some additional simulation runs are performed at a slightly reduced temperature , where polymer transfer is slightly reduced. In this case a completely reproducible fcc structure is found for relatively small system sizes, but no fcc lattices could be stabilized for larger ones as shown in the results in Table 5.
No. configs  Targeted  Ordering  
2  16 micelle bcc  360  16  textbook bcc 
1  32 micelle fcc  720  32  textbook fcc 
1  32 micelle fcc  720  35  distorted bcc 
2  32 micelle fcc  720  unstable  none 
1  54 micelle bcc  1215  54  textbook bcc 
3  54 micelle bcc  1215  55  distorted bcc 
2  108 micelle fcc  2430  unstable  none 
No. configs  Targeted  Ordering  
5  16 micelle bcc  222  16  textbook bcc 
5  32 micelle fcc  444  unstable  none 
5  54 micelle bcc  750  54  textbook bcc 
5  108 micelle fcc  1500  unstable  none 
3  128 micelle bcc  1778  128  textbook bcc 
1  128 micelle bcc  1778  136  distorted bcc 
3  250 micelle bcc  3472  unstable  none 
4  432 micelle bcc  6035  unstable  none 
No. configs  Targeted  Ordering  
4  32 micelle fcc  444  32  textbook fcc 
1  32 micelle fcc  444  unstable  none 
2  108 micelle fcc  1500  unstable  none 
2  108 micelle fcc  1500  103  distorted bcc 
1  108 micelle fcc  1500  104  distorted bcc 
Iv Dynamics of crystal formation
Molecular dynamics simulations not only allow the prediction of equilibrium phase diagram, but also describe the dynamics of micelles as they evolve towards thermal equilibrium. In the simulations presented previously, micelles quickly form from a completely random configuration of polymers, and then after a long simulation time (10 to 20 million time steps) order into a micellar crystal. We now present a quantitative picture of the dynamics of the polymers and micelles as the system approaches equilibrium. All results in this section are obtained from the analysis of the simulations of the polymer.
iv.1 Polymer transfer is an activated process
Polymer transfer plays an important role in achieving the formation of micellar crystals in the simulations discussed above. If there is too little, the single micelle degrees of freedom do not reach equilibrium. Too much pushes the system into a disordered state. Moreover, the rate of polymer transfer is extremely sensitive to the kinetic temperature. An equilibrated bcc micellar crystal was taken as an initial configuration for additional simulations that continued with ranging from to . Figure 6 shows the results. The rate of polymer transfer starts near 0 at and increases exponentially, following an Arrhenius form
(6) 
That is to say that polymer transfer is an activated process. Curve fitting, we find
(7) 
Figure 6 also includes results from simulations performed using a Langevin thermostat (25). It controls the temperature by adding an additional force to every particle where the magnitude of the random force and set the temperature through the fluctuation dissipation theorem (25). The results for two different values of , which are plotted in Figure 6, show the dependence of the polymer transfer for different drag coefficients.
In Figure 6, the slope of the lines in the inset plot are universal for both thermostats and both values of . The universality of this value implies that the calculated is the free energy barrier between a polymer in a micelle to a transition state between micelles. While their slopes are universal, the yintercepts (related to ) should depend on the diffusion coefficient of the hydrophilic beads, which, from the Einstein relation (26) is inversely proportional to the drag coefficient. This is reflected in the offsets of the various plots in Figure 6, which are clearly different.
The difference in free energy between a polymer within a micelle and in the transition state entirely surrounded by solvent and hydrophilic beads can be roughly estimated as , where is the number of hydrophobic beads exposed to solvent. If the length of hydrophobic block is low, . Thus for the polymer analyzed here, it is expected that , which is consistent with the measured value. Hydrophobic beads in polymers with a longer hydrophobic block are expected to form a globule in the transition state, leading to a slower increase of going as . Assuming that remains constant as increases, then maintaining the same rate of polymer transfer will require increasing the kinetic temperature by the same factor. These considerations agree remarkably with the kinetic temperature of selected for the polymer simulations,
.
However, as the size of the hydrophobic block grows, the transition state should also have an additional contribution due to the free energy cost of passing a globule through the brush formed by the hydrophilic coronas. So this simple estimate should eventually break down.
Implicit in our arguments is that polymer transfer accounts for the diffusive behavior in Figure 1. Therefore, the diffusion coefficient can be roughly estimated from the time it takes a polymer to travel the nearest neighbor micelle distance . Comparing this estimate with the fit to Figure 1, a good agreement is obtained:
(8)  
iv.2 Dynamics of micellar crystal formation
The structure factor is sufficient as an order parameter to determine if the entire system is in an ordered state, but it reveals no information about the dynamics before that final state is formed. For this, we turn to the bond order analysis (27) and apply it to all the micelles at every time step in simulation runs performed on the polymer with . In short, at any given time step where a micellar crystal may have partially formed, the bond order analysis identifies those micelles belonging to the ordered and the disordered phases.
One simple way to examine the results is to count the number of ordered micelles in the simulation box at every time step. An example from one simulation run is shown in Figure 7, plots for all other simulation runs performed are qualitatively very similar, though the ordered phase appears at different times.
After starting the simulation shown in Figure 7 from a random configuration, micelles quickly form in only a few thousand time steps and single micelle degrees of freedom are equilibrated before one million time steps have passed. The system then explores configuration space as polymers are transferred and micelles travel anywhere from a few up to without any micelles appearing ordered until time step 10 million. At this point, the number of ordered micelles grows very quickly over the next one million steps until all 128 micelles are in the bcc lattice and remain there for the duration of the simulation.
During this short time span of growth, only of the polymers in the box are transferred between micelles, such a small amount that it cannot fully account for the ordering of all the micelles. During the same time interval, a detailed analysis shows that some micelles move significantly (up to ) before the ordered phase finishes forming. This may suggest that micellar crystal formation is a two step process, where first individual micelles are equilibrated by polymer transfer followed by a second step where polymer transfer becomes irrelevant and the actual crystal grows via the movement of micelles.
However, some additional simulations performed to test this hypothesis do not support it. Single micelle degrees of freedom in these tests are first equilibrated at for 5 million time steps and then the kinetic temperature is quenched to to significantly reduce polymer transfer and allow micelle movement. None of these simulation runs resulted in the formation of an ordered phase. We therefore conclude that a significant amount of polymer transfer remains a critical component in the actual growth of ordered micellar crystals.
V Properties of micellar crystals
v.1 Structure factor
We calculate the structure factor ((4)) for the equilibrium state of every simulation run performed and use it as an order parameter to determine if a micellar crystal is present. In those systems where the crystal does form it is a perfect single crystal and, correspondingly, a large number of peaks can be identified in as shown in Figure 8. Peaks occur in reciprocal space at discrete points , where are the reciprocal vectors of the corresponding lattice.
Peaks evident in Figure 8 decrease in magnitude for larger values of . This damping is expected to follow a DebyeWaller factor (28) approximately described as
(9) 
where is the mean square displacement of the micelle center of mass from its ideal lattice position. A curve fit, shown in Figure 8, is in excellent agreement with (9). The Lindemann ratio is defined as
(10) 
where is the nearest neighbor distance between micelles. Using the parameters of the curve fit yields .
The Lindemann parameter can alternatively be computed directly by measuring the mean square displacement of micelles about their average lattice positions. The results, shown in Figure 9 agree remarkably well with the estimate from the DebyeWaller factor. Furthermore, it has been established empirically that approximately at solids melt into disordered states (29), a result that is also supported from our simulations as we observed no micellar crystals form at kinetic temperatures greater than
v.2 Structure of the lattice of micelles
Micelles in the polymer system are arranged with the centers of mass sitting on a bcc lattice with a nearest neighbor spacing of . This number remarkably agrees with the length of the polymer if stretched completely into a straight line across the diameter of a circle from the center of one nearest neighbor micelle to the opposite one, , where is the equilibrium bond length. This suggests that polymers are maximally stretched, either along the diameter or in a bent configuration. Detailed visual examinations of simulation snapshots indicate that there are a significant number of polymers in the bent configuration, but there are also some in the linear extended configuration. The micellar cores are liquidlike, and over time, a given polymer constantly switches from linear to bent configurations.
The previous arguments also suggest a significant amount of overlap of the coronas between two neighboring micelles in the lattice. To examine this more quantitatively we calculate the micelle density distribution as a function of the radius averaged over all micelles and all time steps in the simulation after the micellar crystal has formed. Figure 10 shows the results, confirming the overlap. Hydrophilic beads from one micelle explore the solvent until occasionally bumping into the hydrophobic core of one of the nearest neighbor micelles.
A remarkable aspect of micellar crystals found in this work is the stability of the average aggregation number. If the polymers within the micelle are maximally stretched, the core of the micellar radius is , so the aggregation number can be estimated from
(11) 
For the this yields , in good agreement with the simulation results in Table 1. Aggregation numbers for the other polymers simulated do not agree, implying that the polymers in those systems are not maximally stretched.
Vi Conclusions
vi.1 Summary of results
Cubic micellar crystals of polymers form in MD simulations at sufficiently high concentrations. In order to form the crystals, high enough kinetic temperatures are needed to enable polymer transfer between micelles, which is critical for equilibrating the system. The polymer transfer process is activated and described by transition state theory. It results in an apparent diffusive behavior of the individual polymers while the lattice of micelles remains stable. Excessive polymer transfer at even higher kinetic temperatures triggers a disorder phase transition to a micelle liquid.
In the process of forming the ordered phase, the system spends a long time in a micellar liquid phase equilibration period where no ordered nucleates are present. Eventually, a large nucleation event takes place and the micellar crystal then grows very quickly, filling the entire simulation box in a relatively short time span. During this growth period, polymer transfer and movement of micelles are both crucial in the formation of the final micellar crystal. The preferred ordering near the disordered transition for all triblocks studied in this system is the bcc lattice. Only at lower kinetic temperatures and for polymers with short hydrophilic groups ( at ) is there some evidence for a stable fcc phase. These results are summarized in Figure 11.
All micellar crystals obtained in this work are perfect single crystals displaying a high degree of order. Even in a periodic simulation box with four unit cells along a side, a single additional micelle can disrupt the resulting lattice significantly. The number of micelles is controlled by changing the number of polymers in the box, as the system displays a remarkable stability in the average aggregation number of the micelles that form. In the equilibrium lattice, micelles are closely packed with a significant amount of overlap between the coronas of neighboring micelles. Polymers in the micelles are highly stretched across the liquid hydrophobic core through the solvent and qualitatively well described within the strongly stretched approximation (30); (17).
vi.2 Stability of the bcc lattice near the disorder transition
Our simulations results show a strong preference for bcc lattices near the disorder transition. A similar result has been experimentally observed for PSPI diblocks (10), where it has been attributed to the fact that near the disordered phase, micelle aggregation numbers are small. The phase diagram of fstar polymer systems shows that fcc lattices are only stable for large number of arms , while bcc lattices are favored when the number of arms is small (36). By considering polymeric micelles as fstar polymers, where the number of arms is given by then bcc lattices are favored when the aggregation numbers are small . This argument, however, hinges on two key assumptions: the dynamic nature of micelles does not play a significant role and that the hydrophilic blocks are sufficiently long. The first assumption is already somewhat problematic given the importance of polymer transfer found in this work. As for the second assumption, a criteria establishing how long hydrophilic blocks should be has been put forward in Ref. (8), where it is shown that if the size of the corona is and the core radius , bcc lattices are favored for . Our results for the and yield and . It is therefore not possible to attribute the stability of the bcc lattices observed in our simulations as being a consequence of the small aggregation numbers .
It is tantalizing to interpret the stability of the bcc lattice in terms of the AlexanderMcTague (AM) scenario (37), where it was argued that bcc should be generally expected to be the stable phase near a (weakly first order) disorder transition. Subsequent analysis however, showed that cubic lattices other than bcc cannot be ruled out near the disorder transition (38); (39). In Ref. (38); (39) it is shown that the characteristic property of bcc lattices is that their free energy is closer to the disordered state, thus suggesting that bcc lattices follow an Ostwald step rule (40), where the solid phase that nucleates first is the one whose free energy is closest to the disordered (or fluid) state. In this case, the complete crystallization process would require an additional step, where after the bcc crystallites are formed, they gradually evolve towards the stable thermodynamic phase. Certainly, it follows from our results that fcc lattices are difficult to obtain by MD for the systems we simulate, but at least in one system ( at ) the fcc lattice has been obtained reproducibly, and did not proceed through an intermediate bcc step. In addition, for this very same system, closer to the disordered state (), no fcc structure was found to be stable. Although not completely conclusive, our results are more consistent with the bcc as being a stable thermodynamic phase near the disordered phase.
There are serious limitations in identifying micelles as simple particles, because as the disordered phase is approached micelle aggregation numbers decrease and polymers become essentially free. So the disordered phase is not a simple liquid as it is assumed by AM and all subsequent work. Similar analysis in polymer melts (41), shows that in the vicinity of the spinodal, the only possible phase with cubic symmetry is bcc. Beyond the spinodal, other structures are favored (42). Based on the previous discussion, we attribute the stability of the bcc phases observed in our simulations as a reflection of the admittedly nonrigorous statement that bcc phases are usually favored near the disordered transition. We defer to future work to establish this result within a rigorous framework, where all the nuances involved in micelle formation are properly taken into account.
vi.3 Implications for Pluronic systems
The and polymers discussed in this paper provide coarsegrained descriptions of Pluronic F127. All simulations have been carried out in very good solvent conditions for the beads and at total volume fractions of 15–20%. Experimental results in this region of the phase diagram are surprisingly disparate. Simple cubic (31), bcc (32) and fcc (33) have all been proposed as the structure in this region. Very recently, on the basis of new experimental data, the situation has been thoroughly reviewed by Li et al. (34) (although for significantly higher temperatures), see also Ref. (35), but without clear conclusive results. Our theoretical analysis clearly favors the bcc lattice close to the disorder transition. The comparison of our results with the experimental F127 system is more accurate at low temperatures (at 2025 C), where the water can be considered as a good solvent for PEO, as discussed in Ref. (18).
vi.4 Outlook
We have shown that MD allows a detailed investigation of both the dynamics as well as the thermodynamic equilibrium of micellar crystals. Many studies have been performed by modeling micelles as point particles, where the complex structure of the micelles is accommodated through refined twobody potentials, either derived analytically or empirically. While successful in many situations, twobody potentials do not account for the dynamic nature of micelles, which play a critical role in determining the phases of the system, particularly near the disorder transition.
There are a few areas where further work needs to be performed. First, all simulations in this work have been performed in good solvent conditions. MD with implicit solvents of different quality are also of great interest, especially to determine the phases of Pluronic systems over a wide range of temperatures (18). Next, the largest micellar crystal formed in this work contains 2134 polymers (57,618 beads). We were unable to find any order in larger systems, even after running as many as 50 million time steps. It is possible that significantly longer simulations may be required for larger systems. Also, the range of applicability of MD is restricted to a relatively narrow range near the disorder transition, as schematically shown in Figure 11. Finally, the role of the solvent will need to be investigated in more detail, as it is found in the Rouse vs. Zimm dynamics for simple homopolymers (26). Future studies will be necessary to completely clarify and expand all these issues.
The relevance of our study goes beyond pure systems. In Ref. (43) for example, Pluronic polymers have been used to template the growth of an inorganic phase of calcium phosphate, aimed at creating new polymer nanocomposites with lightweight/high strength properties or that mimic the structure of real bone. An understanding of the pure systems is clearly a prerequisite for accurate models of polymer nanocomposites. We hope to report more on this topic in the near future.
Acknowledgements.
We acknowledge interest and discussions with S. Mallapragada, K. SchmidtRohr and G. Grason. We also thank J. Schmalian for clarifications regarding the AlexanderMcTague paper. This work is supported by DOEBES through the Ames lab under contract no. DEAC0207CH11358 and partially supported by the NSF through grant DMR0426597.References
 P. Alexandridis and B. Lindman Edtrs., Amphiphilic block copolymers : Selfassembly and applications, Elsevier, Amsterdam (2000).
 C.N. Likos, Soft Matter 2, 478 (2006).
 G. Wanka, H. Hoffmann and W. Ulbricht, Macromolecules 27, 4145 (1994).
 B. Chu, Langmuir 11, 414 (1995).
 K. Mortensen, Coll. and Surf. A 183, 277 (2001).
 K. Mortensen, W. Brown and E. Jorgensen, Macromolecules 27, 5654 (1994).
 G.A. McConnell, A.P. Gast, J.S. Huang and S.D. Smith, Phys. Rev. Lett. 71, 2102 (1993).
 G.A. McConnell and A. Gast, Phys. Rev. E 54, 5447 (1996).
 J. Bang, T. P. Lodge, X. Wang, K.L. Brinker and W.R. Burghardt, Phys. Rev. Lett. 89, 215505 (2002).
 T.P. Lodge, B. Pudil and K.J. Hanley, Macromolecules 35, 4707 (2002).
 T.P. Lodge et al., Phys. Rev. Lett. 92, 145501 (2004).
 T.P. Lodge and J. Bang, Phys. Rev. Lett. 93, 245701 (2004).
 B. van Vlimmeren et al., Macromolecules 32, 646, (1999).
 Y. M. Lam and G. GoldbeckWood, Polymer 44, 3593, (2003).
 X. Zhang, S. Yuan and J. Wu, Macromolecules 39, 6631 (2006).
 P. Ziherl and R.D. Kamien, J. Phys. Chem. B 105, 10147 (2001).
 G.M. Grason, J. Chem. Phys. 126, 114904 (2007).
 J. Anderson and A. Travesset, Macromolecules 39, 5143 (2006).
 V. Ortiz et al., J. of Poly. Sci. B 44, 1907 (2006).
 S.J. Plimpton, J. Comp. Phys. 117, 1 (1995).
 W. G. Hoover, Phys. Rev. A 31, 1695 (1985)
 D. Bedrov, C. Ayyagari and G.D. Smith, J. Chem. Theory Comput. 2, 598 (2006).
 K. SchmidtRohr, J. Appl. Cryst. 40, 16 (2007).
 DeLano, W.L. The PyMOL Molecular Graphics System, 2002, on World Wide Web http://www.pymol.org
 D. Frenkel and B. Smit, Understanding Molecular Simulations, Academic Press, San Diego CA, 2002.
 M. Doi and S. edwards, The theory of Polymer Dynamics, Clarendon Press, Oxford, UK, 2001.
 P. J. Steinhardt, D. R. Nelson and M. Ronchetti, Phys. Rev. B, 28, 784 (1983).
 C. Kittel, Introduction to Solid State Physics, John Wiley and Sons, Hoboken NJ, 2005.
 J.P. Hansen and I. R. McDonald, Theory of Simple Liquids, Academic Press, London UK, 1986.
 A. Semenov, Sov. Phys. JETP 61, 733 (1985).
 R.K. Prud’homme, G. Wu and D.K. Schneider, Langmuir 12 4651, (1996).
 K. Mortensen and Y. Talmon, Macromolecules 28 8829, (1995).
 T. Liu and B. Chu, Journal of Applied Crystallography 33 727, (2000).
 Y. Li et al., J. Phys. Chem. B. 110 26424 (2006).
 D. Pozzo and L. Walker, Coll. and Surf. A, 294, 117 (2007).
 M. Watzlawek, C.N. Likos and H. Lowen, Phys. Rev. Lett. 82, 5289 (1999).
 S. Alexander and J. McTague, Phys. Rev. Lett. 41, 702 (1978).
 B. Groh and B. Mulder, Phys. Rev. E 59, 5613 (1999).
 W. Klein, Phys. Rev. E 64, 056110 (2001).
 W. Ostwald, Z. Phys. Chem. (Leipzig) 22, 289 (1897).
 L. Leibler, Macromolecules 13, 1602 (1980).
 C. Marques and M. Cates, Europhys. Lett. 13, 267 (1990).
 D. Enlow et al., Mater. Chem. 17, 1570 (2007).