J. Phys. Chem. B 2008, 112, 5423-5427
5423
Low-Radii Transitions in Co-assembled Cationic-Anionic Cylindrical Aggregates Marc Michael Del Rosario Lim, Yuri S. Velichko, Monica Olvera de la Cruz, and Graziano Vernizzi* Department of Materials Science and Engineering, Northwestern UniVersity, EVanston, Illinois 60208 ReceiVed: October 31, 2007; In Final Form: January 17, 2008
We investigate the formation of charged patterns on the surface of cylindrical micelles from co-assembled cationic and anionic amphiphiles. The competition between the net incompatibility χ (which arises from the different chemical nature of oppositely charged molecules) and electrostatic interactions (which prevent macroscopic segregation) results in the formation of surface domains. We employ Monte Carlo simulations to study the domains at thermal equilibrium. Our results extend previous work by studying the effect of the Bjerrum length lB at different values of the cylinder’s radius R and χ and analyze how it affects the transition between helical, ring, and isotropic patterns. A critical surface in the space (lB, R, χ) separating these three phases is found, and we show how it corresponds to a first-order phase transition. This confirms that the Bjerrum length lB is a significant parameter in the control of the helical-ring transition; the ring pattern is strongly associated with short-range forces, whereas the helical pattern develops from dominant long-range electrostatic interactions.
1. Introduction Many heterogeneous molecules, such as surfactants or blockcopolymers, self-organize in solution into finite size aggregates such as micelles, vesicles, bilayers, or more complex structures.1,2 The unique properties, functionality, and striking potential applications of such nanoaggregates have attracted much attention.3-8 The main driving forces toward self-assembly in these systems are long-range electrostatic interactions, shortrange van der Waals interactions, and the molecular structure of the heterogeneous molecules. The competition among those forces results in two different levels of organizational complexity. First, it leads to the spontaneous formation of stable finitesize aggregates of different morphologies.9 Second, it drives the formation of surface domains of the two molecular species.10 The latter effect has been observed and is well explained for flat lipid bilayer membranes.11-14 A similar description of domain formation induced by electrostatics on curved surfaces is less-understood, because of the complex interplay between geometric curvature of the surface and its charge density distribution. In this paper, we study immiscible species of cationic and anionic amphiphilic molecules that co-assemble into a cylindrical micelle. As has been shown in related works on cylindrical micelles,15-17 such surface domains tend to organize in regular patterns including isotropic, ring, and helical/ lamellar symmetries.16 It has also been shown that the different phases are controlled by few parameters, namely the cylinder radius R, the average dielectric constant , the net degree of compatibility χ (the standard Flory-Huggins parameter), and salt-induced inverse screening length κ.18 In particular, the competition of electrostatic forces versus immiscibility forces yields the following general limit behavior: if no electrostatics is present, the system macroscopically segregates in two components. If no immiscibility is present, the system organizes in a regular ionic crystalline structure. If both forces are * To whom correspondence should be addressed. E-mail: g-vernizzi@ northwestern.edu.
vanishingly small, an entropic random mixing of the two components dominates. However, in all former studies of this model, the precise role of the average dielectric constant has only been partially included by using scaling arguments. In this work, we want to shed light on how the phase diagram changes by varying the average dielectric constant . In particular, we are interested in understanding what happens to the transitions between the isotropic, helical, and ring phases as a function of . However, because an overall dielectric constant can be “hidden” in continuum models by adopting dimensionless units, we will explicitly consider a discrete coarse-grained model with a lattice constant that fixes all dimensionless units. We note that the problem on which we are focusing is not the one of changing the relative dielectric constant of the medium inside the cylinder with respect to the one outside. This would require a more complex treatment, involving both the solution of the Poisson equation inside and outside the cylinder, and the matching of the solutions with suitable boundary conditions at the cylinder’s surface.19 Our aim is instead focused on understanding the extent to which long-range electrostatic interactions are responsible for the particular symmetry patterns and phase transitions among them. From this point of view the average dielectric constant plays the effective role of damping the long-range Coulombic forces with respect to the short-range interactions controlled by χ. Because throughout this article we keep both the temperature and the charge valency constant, we consider variations of by means of the Bjerrum length lB ) q2/(4πKBT), where q is the valency of the charges, KB is the Boltzmann constant, and T is the absolute temperature. Our results are therefore limited to the effects of varying the Bjerrum length and do not include general variations of the dielectric constant independently from the temperature and charge valency. The outline of this article is as follows: the next section is devoted to the description of the model, the details of the simulation we use to solve it numerically, and the order parameters we use to analyze different phases. In Results and Discussion, we present our results, particularly the details of
10.1021/jp7105132 CCC: $40.75 © 2008 American Chemical Society Published on Web 04/08/2008
5424 J. Phys. Chem. B, Vol. 112, No. 17, 2008
Del Rosario Lim et al. N
lB
∑ i>j
qiqj rij
)
{
2
γlBN Lz
+
lB Lz
the transitions among isotropic, ring, and helical/lamellar symmetries. All of the results are summarized in a phase diagram, which allows a comparison of our findings with closely related studies of this system. We conclude in the last section by highlighting the physical consequences of our findings and by discussing possible applications. 2. Model and Simulation Technique A model for studying the charged domains over cylindrical micelles has been introduced in ref 15. Its basic assumptions are that each aggregate is a stable cylindrical structure, while all possible shape fluctuations such as bendings, bulgings, twistings are neglected. It also assumes that the surface number density is constant, and that the system is composed of oppositely charged units (q, with the constraint of global electroneutrality. Such a model aims to describe mainly the local order of the surface domain patterns and their thermodynamical properties. The charges are assumed to be confined on a cylindrical monolayer of radius R and infinite length. For simulational purposes the cylindrical micelle is placed in the center of a box of size L × L × Lz, with the z-axis being the cylinder axis, and Lz the cylinder length. Periodic boundary conditions are imposed along the z-direction. The charges on the surface of the cylinder are positioned at the vertices of a triangular lattice (see Figure 1), with one lattice vector orthogonal to the cylinder axis. Thus, each charge has six equidistant neighbors, i.e., the charges are in a fully packed configuration. The degree of incompatibility between the two molecular species is controlled by the Flory-Huggins parameter, χ ≡ (ω++ + ω- - - 2ω+-)/(2KBT), where ω(( are all possible net pair-interaction energies between positive and negative charges. In particular, we choose ω- - ) ω+- ) ω-+ ) 0 and ω++ g 0. The Hamiltonian of the system is the sum of a Coulomb term and an immiscibility term, i.e.
H KBT
N
) lB
∑ i>j
qiqj rij
+
1
N
∑ Wij, 2 i,j
(
ω 0 W ) ++ 0 0
)
(1)
where qi is the sign of the charge of the i-th unit (i.e., qi ) (1), rij ) |r bi - b rj| is the distance between the i-th and j-th units at position b ri and b rj, respectively, and the sum over {i, j} is over all first-neighboring pairs. Because periodical boundary conditions are applied only in the z-direction, the Lekner summation technique20 can be used to calculate the electrostatic energy more efficiently, i.e.
qiqj × ∑ i*j
∞
∑ cos(2πkzjij)K0(2πkFjij) - ln k)1
-ψ(zjij|) -
Figure 1. Snapshots of charge distributions for the (a) isotropic phase (R ) 2.5, ) 25, χ ) 3), (b) helical phase (R ) 2.5, ) 25, χ ) 12), (c) intermediate phase (R ) 2.0, ) 30, χ ) 8.5), and (d) ring phase (R ) 2.0, ) 30, χ ) 12).
N
π 2
() Fjij 2
, Fjij * 0
cot(π|zjij|),
(2) Fjij ) 0
where γ ≈ 0.577 is the Euler constant, K0 is the modified Bessel function, ψ is the digamma function, jrij ) {xij, yij, zij}/Lz, Fjij )
xx2ij+y2ij/Lz, and for m ) 0 the terms with i ) j are omitted. In
our numerical simulations, the infinite sum in eq 2 is truncated by following the indications given in.20 Specifically, the electrostatic energy can be computed with an accuracy that grows with increasing Lz/R. For instance, an accuracy of 1 × 10-5 for Lz/R ≈ 1 × 103 requires about 2000 terms. In our case Lz/R ≈ 1 × 102, which requires about 300 terms to achieve a precision of 1 × 10-5. The equilibrium thermodynamics of this system can be computed by standard Monte Carlo simulations. The positions of the charges are fixed, only the type of charge (+q or -q) at each lattice site constitutes the degrees of freedom of the system. Monte Carlo moves in the phase space are performed by simply exchanging two randomly chosen particles according to the Metropolis scheme.21 The system reaches the equilibrium typically over 5 × 103 MC steps per particle and afterward another 1 × 105 MC steps are used to perform measurements. The equilibration process is accompanied by a gradual decrease of temperature (temperature quenching) from Tmax ) 5 to Tmin ) 1, in temperature units with KB ) 1. Our simulations promptly find several charge patterns with different morphologies, i.e., isotropic, helical, ring, or transitional (see Figure 1). At thermal equilibrium, the system is characterized by its internal energy 〈E〉 and the heat capacity
CV )
1 〈(E - 〈E〉)2〉 KBT
(3)
We run a series of simulations for different values of ω++ ∈ [0, 15] (which corresponds to χ ∈ [0, 30]), ∈ [10, 30], R ∈ [1.2, 2], and Lz/σ ) 100, with σ being the triangular lattice constant). In order to determine the exact value χ ) χ* at which the transition between the helical and ring phase occurs (at a given cylinder radius R and dielectric constant ), we adopt the asphericity parameter A introduced in ref 22. It is defined as follows. First, one makes a list of all clusters C forming connected domains on the cylinder surface. Second, one computes the inertia matrix TRβ for each independent cluster C
TRβ(C) )
1
NC
∑ (rRi - rβi )(rRj - rβj )
N2C i>j
(4)
where NC is the number of charges forming the domain, ri ) (rxi , ryi , rzi ) is the position vector of the i-th charge in threedimensional Cartesian coordinates. The asphericity parameter A(C) for each individual cluster is defined as follows
A(C) )
2 1 3 〈Tr[T (C)]〉 2 〈(Tr[T(C)])2〉 2
(5)
Low-Radii Transitions in Cationic-Anionic Aggregates
Figure 2. The specific heat CV as a function of χ, at different values of the Bjerrum length, and for (a) R ) 1.5 and (b) R ) 2.
where Tr[T] is the trace of the matrix T. Finally, the asphericity parameter A is averaged among all clusters of a configuration. For spherically symmetric objects (isotropic phase) the inertia matrix is proportional to the identity matrix, and the asphericity parameter A is equal to zero. It is also easy to verify that for thin annular objects (ring phase), one has A = 1/4. Long thin stripes (lamellar phase) correspond to A = 1. The asphericity parameter A is therefore a good order parameter for analyzing the transitions between different phases. 3. Results and Discussion To characterize the phase behavior, we calculated the internal energy, specific heat, and surface domain asphericity parameter at several values of the Flory-Huggins parameter χ and the Bjerrum length lB. All the results are summarized in a phase diagram, where phases with domains having different symmetries are plotted. The specific heat CV as a function of χ (Figure 2) clearly shows a broad peak around a region between the isotropic and helical phases. This is a signal of a possible transition, corresponding to the topological rearrangement of the charge patterns on the cylindrical surface (order-disorder phase transition).23 The peak moves slightly toward larger values of χ with increasing Bjerrum length lB (i.e., decreasing dielectric constant). However, there is no significant dependence on the value of the radius R. Moreover a second, sharp peak appears where the helical phase meets the ring phase. Such a transition seems to be strongly dependent on both the radius of the cylinder R, and the Bjerrum length lB (see Figure 2). It is well-known that the phase transition from isotropic to periodic structures cannot be continuous in three dimensions.24 In two-dimensional binary systems with 1/r Coulomb interactions, a classical phase transition from isotropic to periodic structures is uncertain because the long-range order in such systems is questionable.13,25-27 It was suggested that the melting of striped structures may occur via the Kosterlitz-Thouless mechanism.28 We expect to recover the two-dimensional limit when R f ∞. The second peak seems to soften at larger values of lB. This might indicate a possible phase where the helical and the ring phases merge, and the equilibrium configuration is an ordered
J. Phys. Chem. B, Vol. 112, No. 17, 2008 5425
Figure 3. Average asphericity of the domains of positive charges as a function of the Flory-Huggins parameter χ, and for (a) R ) 1.5 and (b) R ) 2. A dashed line is drawn at A* ) 0.60 to better illustrate the sharp transition between the helical and ring phase.
ionic-like lattice. On the other hand, at small values of lB (i.e., large value of the dielectric constant), electromagnetic interactions are mostly screened and a macroscopic segregation among charges occurs. The Monte Carlo equilibrium configurations we obtained in this region confirm this behavior. Also the asphericity parameter shows a sharp transition between helical and ring symmetries. In Figure 3, starting from small values of χ, the helical symmetry approaches a maximal asphericity value of 1, though predominant helical character is noted when asphericity is greater than 0.8. The asphericity for the ring symmetry varies between 0.2 and 0.4, with larger rings corresponding to lower asphericity values. We thus define the transition from the helical to ring symmetry at the point where A* ) 0.6. The isotropic symmetry at small values of χ inevitably centers on an asphericity value of 0.65 (Figure 3). This value correspond to a systems of small (planar) ellipsoids with eccentricity ∼0.919. This number can be numerically obtained by averaging the asphericity parameter over all small random clusters on a triangular lattice. Therefore, it is a factor that depends on the lattice local geometry, and it does not depend on the cylindrical geometry (R) or thermodynamics (lB, T). By a direct inspection of the equilibrium configurations from the Monte Carlo simulations, we observed that stronger electrostatic interactions favor helical patterns, while strong short-range interactions favor ring patterns. Such a qualitative observation has been confirmed by recent analytical investigations,17 on the electrostatic origin of spontaneous chiral configurations. In this respect, we also analyzed the distribution of the angles formed by the helical stripes with respect to the cylinder axis. However, the data we collected were not sufficient to identify a possible preferred chiral angle and its dependence on the dielectric constant. On the other hand, our results are also fully consistent with ref 17 at large value of the dielectric constant, where electrostatics is mostly suppressed and the system is macroscopically segregated, with half of the cylinder being of positive charges and the negative ones in the remaining half. We define the critical value χ* at which A ) A*. In all our simulations, we find that the critical value χ* increases as R increases. This effect is also enhanced by the fact that with increasing R, there are more charged molecules in the system; the consequence is that larger charged domains require higher values of χ to be
5426 J. Phys. Chem. B, Vol. 112, No. 17, 2008
Del Rosario Lim et al.
Figure 4. Phase diagram of charged cylindrical patterns in the space (lB, χ). The transition line between helical phase (center) and ring phase (top), has been computed from the jump in the asphericity parameter for R ) 1.5 and it corresponds to a first-order phase transition. The line limiting the isotropic region (bottom) has been computed by the peak in the specific heat and it does not depend on R.
balanced. As we mentioned in the introduction, the whole scenario can be simply understood in terms of the competition between immiscibility and electrostatic interactions. In particular, the helix-ring transition occurs precisely because increasing electrostatic interactions (which favors a isotropic phase) require increasing degree of immiscibility (which favors a segregated phase). The periodicity and stripe length d of the ring phase has been already considered in ref 15. In that paper, it is shown how the free energy F of the ring phase can be written as
2 F ) + ds2 d
(6)
where d is the dimensionless distance among the rings, and s2 is the average over a cell of the electrostatic interaction with a potential of the form V(r) ) lB/r. That is
s2 )
1 2
∫ ∫cell dxdy σ(x)σ(y)V(y - x - Λ) ∑ Λ
(7)
where the integration is over a unit cell of the lattice {Λ} generated by the periodic rings, and σ(x) is the charge density. We computed analytically the expression for s2 in the ring phase, and minimized the free energy F in eq 6 with respect to d. In the small-R limit, we find d ≈ 1/R. By scaling arguments we also obtain d ≈ 1/xlB in the ring phase. We find full agreement of our numerical simulations with such scaling behaviors. Figure 4 collects all our above results in a single phase diagram. A dominant feature is that the dielectric constant (through the Bjerrum length hindering electrostatic interactions) narrows the range over which the helical phase prevails. The line separating the helical phase from the isotropic phase does not depend on the radius R sensibly. In fact, for values of χ < 5, the entropy of the configurations dominates the free energy. Therefore, the effect of the Bjerrum length lB is only marginal in that region. On the other hand, the line separating the helical phase from the ring phase strongly depends on the radius R; its behavior for different values of R is plotted in Figure 5. A general feature of this phase diagram is that when is small, even for very high values of χ > 10, the helical phase persists. On the other hand, for values of χ < 4, the isotropic phase dominates at small values of lB, whereas at large values of lB
Figure 5. Radii dependence of the critical line separating the ring phase (above the line) from the helical phase (below the line). The critical values χ* are computed from the data in Figure 2.
(larger than the ones explored in our simulations), a phase with full ionic order of the domains is expected. Finally, in all of our simulations, the helical phase disappears at small values of R. A word of caution is in order here. The finite-size of the simulation box introduces a spurious periodicity along the cylinder axis. This amplifies further the effect of the commensurability constraint of the patterns with the cylindrical geometry. Such a phenomenon is relevant mainly at large values of the domain size with respect to Lz, where commensurability constraints are dominant. We directly verified the effects of a finite-size simulation cell, by changing the value of Lz. While we did not find any qualitative difference for the behavior of ring and isotropic patterns, we noticed sensible differences for helical patterns for small Lz values. 4. Conclusions In this paper, we considered a simple model for charged patterns formation over a cylindrical geometry. We analyzed its equilibrium thermodynamics by Monte Carlo simulations. We focused on a specific region of the parameter space, where ring, helical and transitional patterns meet one another. Different phases are characterized by different values of the asphericity of the charged domains, which is assumed as the order parameter. Our numerical investigations find a quantitative dependence of the asphericity of the equilibrium domains on the average dielectric constant of the surrounding medium. The helical phase is favored when and χ are small, and R is large. The ring phase is favored where short-range, nonelectrostatic forces are dominant, that is at large χ, large , or low R. In general, the model we considered in this article allows us to analyze the changes of the surface patterns of charges by changing the dielectric properties of the medium, for instance by simply using mixed solvents. These findings may be of interest for bio- and nanotechnological applications. We believe our results might be used for designing cationic-anionic micelles with specific helical, ring, or isotropic charged domains. For instance, these structures could be used for changing the conductivity of nanofibers if metal ions are adsorbed on the surfaces of the fiber; in this case, the ring structure will be an insulator and the helical structure a conductor.
Low-Radii Transitions in Cationic-Anionic Aggregates A potential limitation of our analysis is the assumption of fixed micelle length along the z-axis. In the case of long fibers, we cannot exclude the possibility of positional defects and coexisting not homogeneous phases.28 We postpone the exploration of this aspect and its implications to future studies. Another aspect which might deserve further attention is the inclusion of screening length effects, due to the presence of salt and counterions in the solution.29,30 Acknowledgment. M.M.D.R.L. and G.V. acknowledge the support by the Nanoscale Science and Engineering Initiative of the National Science Foundation under NSF Award Number EEC- 0647560. Y.S.V. and M.O.C. acknowledge the support of the Materials Research Science and Engineering Center program of the National Science Foundation DMR0520513. References and Notes (1) Rubinstein, M.; Colby, R. H. Polymer Physics; Oxford University Press: New York, 2003. (2) Hamley, I. W. Introduction to Soft Matter: Polymers, Colloids, Amphiphiles and Liquid Crystals; John Wiley & Sons: New York, 2000. (3) Veatch, S. L.; Keller, S. L. Phys. ReV. Lett. 2002, 89, 268101. (4) Laradji, M.; Kumar, P. B. S. Phys. ReV. Lett. 2004, 93, 198105. (5) Kaler, E. W.; Herrington, K. L.; Murthy, A. K.; Zasadzinski, J. A. N. J. Chem. Phys. 1992, 96, 6698. (6) Niece, K. L.; Hartgerink, J. D.; Donners, J. J. J. M.; Stupp, S. I. J. Am. Chem. Soc. 2003, 125, 7146. (7) Velichko, Y. S.; Stupp, S. I.; Olvera de la Cruz, M. J. Phys. Chem. B 2008, 112, 2326-2334. (8) Palmer, L. C.; Velichko, Y. S.; Olvera de la Cruz, M.; Stupp, S. I. Philos. Trans. R. Soc. London, Ser. A 2007, 365, 1417.
J. Phys. Chem. B, Vol. 112, No. 17, 2008 5427 (9) Israelachvili, J. N. Intermolecular and Surface Forces, 2nd ed.; Academic Press: London, 1992. (10) Meyer, E. E.; Lin, Q.; Hassenkam, T.; Oroudjev, E.; Israelachvili, J. N. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 6839. (11) Andelman, D.; Brochard, F.; Joanny, J. J. Chem. Phys. 1987, 86, 3673. (12) Zhang, D.; Carignano, M. A.; Szleifer, I. Phys. ReV. Lett. 2006, 96, 028701. (13) Loverde, S. M.; Velichko, Y. S.; Olvera de la Cruz, M. J. Chem. Phys. 2006, 124, 144702. (14) Loverde, S. M.; Solis, F. J.; Olvera de la Cruz, M. Phys. ReV. Lett. 2007, 98, 237802. (15) Solis, F. J.; Olvera de la Cruz, M.; Stupp, S. I. J. Chem. Phys. 2005, 122, 054905. (16) Velichko, Y. S.; Olvera de la Cruz, M. Phys. ReV. E 2005, 72, 041920. (17) Kohlstedt, K. L.; Solis, F. J.; Vernizzi, G.; Olvera de la Cruz, M. Phys. ReV. Lett. 2007, 99, 030602. (18) Debye, P. W.; Hu¨ckel, E. Phys. Z. 1923, 24, 185. (19) Jackson, J. D. Classical Electrodynamics. 3rd ed.; John Wiley & Sons: New York, 1998. (20) Grzybowski, A.; Bro´dka, A. Mol. Phys. 2002, 100, 635. (21) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. J. Chem. Phys. 1953, 21, 1087. (22) Rudnick, J.; Gaspari, G. J. Phys. A 1986, 19, L191. (23) Chaikin, P. M.; Lubensky, T. C. Principles of Condensed Matter Physics; Cambridge University Press: Cambridge, U.K., 1995). (24) Brazovskii, S. A. SoV. Phys. JETP 1975, 41, 85. (25) Mermin, N. D. Phys. ReV. 1968, 176, 250. (26) Baus, M. J. Stat. Phys. 1980, 22, 111. (27) Loverde, S. M.; Olvera de la Cruz, M. J. Chem. Phys. 2007, 127, 164707. (28) Nelson, D. R. Defects and Geometry in Condensed Matter; Cambridge University Press: Cambridge, U.K., 2002). (29) Henle, M. L.; Santangelo, C. D.; Patel, D. M.; Pincus, P. A. Europhys. Lett. 2004, 66, 284. (30) Menes, R. R.; Pincus, P. P.; Stein, B. B. Phys. ReV. E 2000, 62, 2981.
