Topology, length scales, and energetics of surfactant micelles

In aqueous solutions, surfactant molecules spontaneously self-assemble into diverse geometrically complex and dynamically fluctuating morphologies ranging from spherical and cylindrical to very long, flexible wormlike micelles with or without branches^1–7,9–11 and topologically rich knotted structures.^9,10 The diversity in microstructure and rheological properties makes micellar solutions beneficial to numerous applications¹² as hydrofracking fluids in oil industry, turbulent friction drag reducing agents,¹⁰ thickening agents in consumer products, drug carriers in targeted delivery,¹³ and templates to create functional nanofluids.^14–16 

Since the early work of Debye,¹ microstructural transitions in micellar solutions have been investigated both theoretically^4,17–20 and experimentally.^1,2,21,22 It is now well-recog-nized that the molecular structure of co-surfactant or salt has a spectacular effect on the micellar morphology. In particular, aromatic organic salts have stronger binding affinity to the micelles and induce enormous growth of micelles with increasing salt concentration c_s. Consequently, micelles become long and flexible, and entangle even at relatively low surfactant concentrations c_D. Dilute solutions with spherical or short cylindrical structures exhibit Newtonian fluid rheology.^21–23 In contrast, solutions above the overlap concentration ϕ^∗ consist of very long thread-like structures with contour lengths that span from a few tens of nm to several μm and show viscoelastic behavior reminiscent of flexible polymer solutions.^4,21,22 However, unlike polymers, wormlike micelles (WLMs) can merge and undergo reversible breaking at time scales that are detectable both in scattering experiments and simulations. Under non-equilibrium conditions, such as under shear flow, the structure, dynamics, and the resultant rheological properties could change dramatically: two notable examples are shear induced structure (SIS) formation and shear banding that have been often observed in surfactant solutions.^21–29 Due to the presence of such dynamic processes, a quantitative description of the microstructure of micellar fluids is incomplete. While rheology of micellar solutions has long been studied in experiments, molecular-level simulation and modeling are necessary to address how interactions at multiple length and time scales lead to different morphologies and corresponding rheological behaviors. In this paper, we present a systematic simulation study of self-assembly (SA), emerging structures, length scales, energetics, and rheology of a model cationic surfactant solution in the presence of explicit solvent, electrostatic, and hydrodynamic interactions.

Despite decades of research on micellar fluids, a few fundamental questions still remain unanswered. First, although the scaling of average micelle length 〈L〉 for neutral or highly electrostatically screened micelles is known,^4,17,20 how 〈L〉 and micelle length distribution change with the salt to surfactant ratio R = c_s/c_D is not well understood. Contour length measurements of WLMs in the entangled regime by Cryo-TEM as well as light and neutron scattering experiments have not yielded a conclusive topological picture.³ Second, experimentally observed non-monotonic dependence, with multiple maxima, of zero-shear viscosity η₀ on c_s^31–34 cannot be rationalized within the framework of existing theories. Third, the end-cap energy E_c, defined as the excess energy of the surfactants in the hemispherical region compared to the energy of those on the cylindrical body of a micelle, is generally inferred indirectly from rheological measurements.^15,16,36 Various surfactant and salt systems have been studied experimentally to understand the non-monotonic variation of zero-shear viscosity as a function of salt concentration. Examples include cetyltrimethylammonium chloride (CTAC) and sodium chloride (NaCl), CTAC and sodium bromide (NaBr), CTAC and sodium salicylate (NaSal), cetylpyridinium chloride (CPyCl) and NaCl, and CPyCl and NaSal. As shown in our earlier studies,^6–8 organic counterions such as Sal⁻ strongly bind to the positively charged micelle corona. Specifically, the hydrophobic aromatic ring (represented by the SC₄ beads in Fig. 1) penetrates into the hydrophobic micelle interior while the negatively charged $Q a −$ moiety interdigitates between the positively charged surfactant head groups. In contrast, inorganic counterions such as Cl⁻ or Br⁻ bind weakly to the micelle surface. As a consequence, salt-induced morphological transitions (e.g., sphere to rod transition reported in Ref. 6) are significantly more pronounced in systems with organic salts. Further, double maxima in zero-shear viscosity with increasing salt concentration are observed in micelle solutions containing Sal⁻ and not Cl⁻ or Br⁻ ions. The first decrease in viscosity has been attributed to micelles branching. Yet, a clear picture of the underlying microstructures responsible for this remarkable effect is still lacking in the literature. In order to understand the underlying physics behind the non-monotonic variation of zero-shear viscosity with c_s, we simulate aqueous complexes of CTAC surfactants and strongly binding Sal⁻ counterions using a coarse-grained (CG) molecular representations shown in Fig. 1. Figure 2 (Multimedia view) shows the microstructure evolution of surfactant self-assembly in a typical simulation. From these simulations, we enumerate the emerging micelle morphologies, outline a phase diagram, quantify the end-cap energy and recombination time of micelles, and the relevant length scales that help to shed light on the observed anomalous viscosity variations with respect to R.

The paper is organized as follows. In Sec. II, we describe the coarse-grained molecular model and simulation methods. In Sec. III, we outline a phase diagram of micelle morphologies as a function of surfactant and salt concentrations at ambient temperature and pressure. In Sec. IV, we present an analysis of the micelle shape parameter using differential geometry tools. In Sec. V, we discuss the distribution of micelles length at different salt concentrations. We also discuss coarsening and recombination dynamics of the micelles. We numerically quantify the persistence length and the equivalent 1D bending modulus κ and the end-cap energy of micelles. In Sec. VI, we present the effect of salt concentration on the zero-shear viscosity of a concentrated solution and discuss the mechanisms underlying the non-monotonic variation of zero-shear viscosity. In Sec. VII, we summarize the conclusions of this study. Finally, in the  Appendix, we present an algorithm developed to measure the contour length of the micelles and provide details of supplemental animations.

In the simulations, we represent the surfactant (CTAC), salt (NaSal), and water by CG beads. We adopt the MARTINI force field³⁷ which has been successfully employed previously to model shape transitions⁶ and shear-induced dynamics⁸ in the same system. In this scheme, CTA⁺ ions are represented by five beads, Sal⁻ by four beads, and water by a single bead as illustrated in Fig. 1. We simulate the dissociated state of CTAC and NaSal in water, and charge neutrality in the system is imposed by adding one Cl⁻ ion per CTA⁺ and one Na⁺ ion per Sal⁻ ion. The non-bonded CG beads i and j interact through a Lennard-Jones potential

truncated and shifted at r_ij = 1.2 nm, where $r i j =| r → i − r → j |$⁠, σ_ij and ϵ_ij are the effective diameter of the bead and depth of the potential well, respectively. We use the effective Lennard-Jones parameters σ_ij and ϵ_ij as defined in Ref. 37.

To ensure the connectivity of the component beads that are bonded together [see Fig. 1], we impose a harmonic potential between two neighboring beads as

where K_bond = 1250 kJ mol⁻¹ nm⁻² is the spring constant. An equilibrium bond distance r₀ = 0.47 nm is used for the CTA⁺ beads. The bond length between the CG beads of Sal⁻ within the ringed part is constrained to r₀ = 0.27 nm, while that between the charged bead and the ring particle adjacent to it is fixed at 0.47 nm. The LJ interaction between the bonded atoms is turned off in all simulations.

The rigidity of a chain is controlled by an angular interaction between triplets of atoms. The corresponding potential is

where the force constant K_angle = 25 kJ mol⁻¹ and the equilibrium angle θ₀ = 180°. We model the electrostatic interaction between charged species q_i and q_j as

where ϵ₀ is the permittivity of free space and ϵ_r is the relative dielectric constant of the medium set to 15.0 following the standard practice. Detailed description of the force field parameters can be found in Ref. 37. Simulations were performed using the LAMMPS³⁸ molecular dynamics package. We use the Particle-Particle Particle-Mesh (PPPM) solver to compute long-ranged electrostatic interactions available as a built in package in LAMMPS. Our simulations start from a random distribution of surfactant, salt, and water molecules in the simulation box. The initial configurations are prepared using the Packmol software.⁴³ The density is adjusted by simulating the system in an NPT ensemble (P = 1 atm) for 10 ns at 300 K. After the initial equilibration, the equations of motion are integrated for a constant NVT ensemble with the temperature controlled via a Nose-Hoover thermostat with a time step Δt = 15 fs. To avoid the freezing of the solvent beads, in each simulation, 10% of the water beads were modeled as antifreeze particles as described in the MARTINI force field.³⁷ Since the hydration level in the simulations is high, we observe various disordered and topologically rich self-assembled structures at different R as shown in Fig. 3.

As discussed before, molecular structure of the counterions and their interactions with the surfactant molecules play an important role in determining the micelle shape. Therefore, we first study the effect of charge distribution on the morphology of micelles in solution, which can be used to construct an approximate phase diagram for such systems.

To understand the distribution of the various ions in the solution, we calculate the radial distribution functions (rdfs), g(r), between charged beads for c_D = 0.2M and R = 0.8. As shown in Fig. 4(a), the probability of finding Sal⁻ ions in close proximity to the CTA⁺ ions (which constitute a thin film separating the water-hydrocarbon interface) is much greater than that for other charged species. This clearly shows that the binding affinity of Sal⁻ is very large. The probability of finding Cl⁻ near the micelle surface is the second largest as shown in Fig. 4(a). This implies that most of the Cl⁻ ions remain close to the micelle surface although with lower density as compared to that of Sal⁻. This is due to the electrostatic repulsion arising from the Sal⁻ ions.

Figures 4(b) and 4(c) plot the radial distribution functions, g(r), of CTA⁺ and CTA⁺, CTA⁺ and Na⁺, Sal⁻ and Sal⁻, Sal⁻ and Na⁺, as well as Sal⁻ and Cl⁻. The first peak in Figs. 4(b) and 4(c) is approximately at ∼0.47 nm and signifies the closest distance between the CTA⁺ ions and Sal⁻ ions on the micelles surface. In Fig. 4(b), we have shown the radial distribution function of the CTA⁺ with itself. The peak at 1 nm corresponds to the CTA⁺ ions residing on the micelle surface. Similarly, the pronounced peak in Sal⁻ and Sal⁻ rdf in Fig. 4(c) corresponds to the Sal⁻ adsorbed onto the micelles surface. We can infer the distribution of Na⁺ in the solution from the dashed line plot in Fig. 4(c). The relatively smaller peak in the rdf plot of Sal⁻ and Na⁺ in Fig. 4(c) suggests that Na⁺ ions are mostly distributed within the solution although there is a finite probability of finding them near the micelle surface next to Sal⁻ due to their electrostatic attraction with the adsorbed Sal⁻ ions. Figure 4(e) shows a typical configuration obtained from the simulation which is consistent with the above inference.

Figure 5(a) outlines the phase diagram as functions of c_D and c_s. We determine the boundary separating the dilute and semidilute regime from the overlap volume of the micelles. Assuming each micelle as a rigid cylinder, we compute the overlap volume fraction $ϕ= ∑ i N m L i 3 /V$⁠, where L_i, N_m, and V are the length of the micelle, number of micelles, and the volume of the box. The dilute regime (ϕ < 1) consists of mostly spherical or short cylindrical aggregates. In the opposite limit, above the overlap concentration ϕ ≥ ϕ^∗ ≈ 1, the solution becomes semidilute and consists of wormlike micelles and shows remarkable changes in the contour length and dynamical as well as rheological properties. This is illustrated in Fig. 6, where we have shown 2D views of the microstructure at three different concentrations: (i) ϕ ≪ 1, (ii) ϕ ≈ 1, and (iii) ϕ ≫ 1. At higher concentrations of surfactant and salt, branched or multiconnected structures form. The boundary between the semidilute and branched regimes is determined from the number of nodes $N n = ∑ i N m N n i$ in the solution. We consider solutions that consist on average of 30% branched micelles as the branched regime; see Fig. 7 (Multimedia view) and the discussion in the  Appendix on contour tracking that is utilized to count the number of branches in a micelle. In a limited range of intermediate concentrations, we find micelles with handles [Figs. 3(e) and 3(f)]. These topologically rich toroidal micelle structures form via end-cap attachment of a flexible cylindrical micelle as illustrated in Fig. 8 (Multimedia view).

The mechanisms of micelle shape transformation can be quantified by analyzing the free energy of an interfacial surfactant film when the film thickness is small compared to the neighboring layers. Mathematically, a surface embedded in ℝ³ can be fully characterized by the mean H, and the Gaussian K curvatures. In the small curvature limit, the free energy of the film is given by

where κ and $κ ̄$ are the bending moduli and H₀ is the spontaneous curvature of the film. For an amphiphile with surface area a and liquid hydrocarbon volume v, the packing parameter v/al that minimizes the free energy in Eq. (5) is given by¹⁸ 

where l ≈ 2.0 nm is the length of the hydrocarbon tail. Focusing on the micelles shapes, we numerically calculate v/al by mapping a small section of the interface to a Monge patch as illustrated in Fig. 3(g). A Monge patch is a local surface of the form x(u, v) = (u, v, h(u, v)), where (u, v) is the parametric position of a point and h(u, v) is the height function that is differentiable with respect to u and v. With this parametrization, the mean curvature is defined as

where R₁ and R₂ are the principal radii of curvatures of the patch and the subscripts refer to the partial derivative with the respective coordinates. Similarly, the Gaussian curvature is defined as

For each micelle, we select a vertex (head group) and make a patch of surfactant head groups and the associating counterions within a radius of 2r_hh, where r_hh is the average distance between the head groups [Fig. 3(g)]. In our analysis, we define the normal at each vertex as the average orientation of the surfactant molecules belonging to that patch. It is interesting to note that the normals at the vertices are unambiguously defined from the orientations n of the molecules even at the umbilical points. We rotate the patch to align the normal vector along the z direction. Then, we fit⁴⁴ the positions of surfactant head groups and the counterions with a surface to the following functional form:

where A, B, C, D, and E are constants. From the extracted coefficients, it is straightforward to calculate the principal curvatures at each point, and hence the packing parameter using Eq. (6). Our discussion is limited to the dilute cases thus excluding the patches with saddles that might be present in branched micelles at higher concentrations. An animation of the micelle surface with the Monge patch is shown in Fig. 9 (Multimedia view).

The average value of v/al from these calculations for c_D = 0.10M over a range of salt concentration 0 ≤ c_s ≤ 0.2M (0 ≤ R ≤ 2) is shown in Fig. 5(b). In this regime, the surfactant concentration is not large enough for micelles to form networks even at higher salt concentrations. Without added salt or at lower salt concentrations, v/al ≈ 0.33 and the thermodynamically most favorable shape is a sphere. The surfactant packing parameter v/al can be derived for simple micelle shapes from geometric arguments.¹⁸ N amphiphiles can be packed into a spherical aggregate of radius R_sphere if the following two conditions are satisfied: $4π R sphere 2 =Na$ and $4 π R sphere 3 3 =Nv$⁠. The maximum radius of a spherical micelle without voids is l. From the above equations, we get v/al = 1/3 as the upper limit of the packing parameter for a spherical micelle. Herein, we quantify this parameter by extracting the mean and the Gaussian curvatures of the micelle surface using Eq. (6). Our results agree well with those obtained from geometric arguments. Similar analysis can be done to the cylindrical aggregates leading to the inequality 1/3 ≤ v/al ≤ 1/2 for the range of their packing parameter. Upon increasing c_s, Sal⁻ ions interdigitate into the micelles surface as shown in Figs. 3(a)-3(f), thereby effectively screening the electrostatic interaction between the charged head groups [Fig. 3(d)]. This effect gradually increases with increasing c_s as shown by the decreasing effective surface area per surfactant in Fig. 5(b). In other words, the energy cost for the splay configuration $∼ ▽ ⋅ n 2$ is large as compared to the uniform orientation of the molecules. This results in a more tightly packed structure leading to growth along the major axis, and cylindrical micelles are energetically favorable. The value of the salt to surfactant ratio, R, at which salt-induced sphere to rod transition is observed in our present work is quantitatively comparable to that reported in Ref. 6, which used a more detailed polarizable water model.

To provide insights into micelle length distribution, we simulate systems consisting of N = 16 000 surfactants (c_D = 0.16M) for R = 0.2, 0.8 and 1.6. These simulations consist of   ≈ 1 × 10⁶ atoms in a cubic box of length 54 nm and were conducted long enough (t > 700 ns) to reach equilibrium states. Figures 10(a)-10(c) shows the distribution of the monomer aggregation number N for different R. When R > 1, the distribution p(N) is an exponentially decaying function as shown in Fig. 10(c). However, p(N) is log-normal rather than exponential for R < 1 as evidenced from Figs. 10(a) and 10(b). The exponential distribution p(N) is the signature of the breakage and recombination kinetics of micelles as shown by Cates and Candau⁴ using a mean field theory. For R < 1, the Sal⁻ ions condense non-uniformly over the micelles. It is possible, therefore, that electrostatic interactions in some micelles is only partially screened as compared to certain others. These results can be compared with the existing theories of Ref. 4, which showed that the interplay between entropy and the end-cap energy gives a broad, exponential distribution of lengths $p ( N ) ∝exp − N / 〈 N 〉$⁠. Our results clearly show that this is not the case for charged micelles when R < 1.

In contrast to the dilute systems, where surfactant exchange plays a significant role in micelle growth, coalescence is the primary mode of growth in highly concentrated systems as illustrated by the molecular dynamics (MD) trajectory in Fig. 11 (Multimedia view). We analyze the coarsening rate by tracking the number of micelles N_m shown in the inset of Fig. 10(b) for two different solutions with R = 0.2 and R = 0.8. We observe a power law scaling: $N m ∼ N m 0 t − α ( R )$ with α(0.2) = 0.3, α(0.8) = 0.6, and α(1.6) = 0.8 (not shown). The observed power law in the growth kinetics suggests that there are at least two stages in the growth process. The faster initial stage is characterized by nucleation and coalescence of smaller micelles and micelles fragments. Longer (and fewer) micelles thus produced tend to recombine less frequently, hence the growth rate is reduced. The increase in the growth rate for higher salt concentration may be attributed to increased electrostatic screening by the salt. In order to gain insights into the time scale of micelle growth, we track the MD trajectory of each micelle for ≈300 ns after an initial run time of 500 ns for R = 1.6: see Fig. 11 (Multimedia view) where we have explicitly visualized two successive combination events of a tagged micelle. In Fig. 10(d), we have shown successive combinations of the micelle by plotting N vs. time. The recombination time τ_c, which is defined as the lifetime of a micelle between two such consecutive events, is ≈65 ns. The distribution of τ_c [inset of Fig. 10(d)] in the solution is nearly exponential as it is the case for the micelle length (Fig. 10(c)).

The global flexibility of a micelle $x= 〈 L 〉 l p$ determines the viscoelastic properties of WLM solutions,³⁹ where 〈L〉 is the contour length and l_p is the distance beyond which orientational correlations are lost. Flexible behavior is observed for x ≫ 1, whereas for x ≈ 1, the micelles behave as rigid rods exhibiting I-N transition under shear flow.^29,30 To the first approximation, 1D bending modulus κ ∼ l_pk_BT. Here, we explored the flexibility of micelles for R = 0.8. From the constructed contour path of the WLM chain (see Fig. 17), we define the orientation vector along the contour. Figure 5(c) shows the tangent to tangent correlations $〈n ( s ) ⋅n ( s + s ′ ) 〉=exp − s ′ / l p$⁠, and l_p at different T. It is remarkable that the predicted values l_p ≈ 25 nm and κ ≈ 250k_BTÅ are both consistent with experiments. Furthermore, the inset in Fig. 5(c) shows that l_p is a constant over the temperature range of 280–315 K.

To quantify the micelle end-cap energy, we conducted simulations at different temperatures and measured 〈L〉 for R = 0.8. Figure 12(a) shows 〈L〉 vs. $1/ R ̄ T$⁠, where $R ̄$ is the molar gas constant and confirms the Arrhenius dependence as shown in Refs. 16 and 17. Using these simulation results, we next calculate the end-cap energy of micelles using an existing theory of micelle growth by MacKintosh et al.¹⁷ which predicts that

where ϕ is a dimensionless surfactant concentration, $E c ′ = 1 2 E c − l B v ∗ 2 r ϕ$ is the renormalized end-cap energy, E_c is the energy required to form two hemispherical end-caps of a neutral micelle which is also known as the scission energy, $E elec = l B v ∗ 2 r ϕ$ is the electrostatic contribution to the scission energy, and r, l_B, and v^∗ are the radius of micelles ≈2.5 nm, Bjerrum length, and the effective charge per unit length of micelles, respectively. From the slope of 〈L〉 vs. $1/ R ̄ T$ plot, we estimate the effective end-cap energy $E c ′ ≈9.8 k B T$ in good agreement with experiments.^15,16,22 In addition to this, we have used an alternative method to calculate $E c ′$ as illustrated in Fig. 12. In the inset, we show the distribution of surfactant interaction energy minus the electrostatic energy (E_int − E_elec), for the micelle shown in Fig. 12(b). The distribution of ΔE = E_int − E_elec is approximately normal as shown by the solid red line, i.e., $∼ A σ exp − ( Δ E − Δ E m ) 2 2 σ 2$ with A ≈ 60.21, ΔE_m = 39.22, and σ ≈ 3.51. Note that the surfactants in the micelle are colored based upon their interaction energies. As expected, the energy cost to pack the surfactants in the spherical cap is larger compared to the energy required to pack them in the cylindrical part. A measure of $E c ′$ is therefore given by the width at half maximum of the energy distribution. This results in the estimate of $E c ′ ≈8.04 k B T$ when averaged over all the micelles in the system, in good agreement with the Arrhenius estimate.

To understand the morphological changes underlying the non-monotonic variation of η₀ with c_s, we simulate systems with N = 4000 (c_D = 0.34M) for 0 ≤ R ≤ 2. Figure 13(a) plots 〈L〉 and η₀ vs. R. Broadly, we classify these variations into three regimes. Region I is the micelle growth regime, which is further divided into two sub-regimes: those below and above ϕ^∗. For ϕ < ϕ^∗, there is a gradual rise of 〈L〉 with increasing R. However, 〈L〉 increases sharply for ϕ ≥ ϕ^∗ in the semidilute regime. In region II, as branched and interconnected structures form, 〈L〉 is widely distributed. At the onset of branching (at the interface between regimes I and II), the system consists of fairly elongated linear micelles with relatively short branches. As we move deeper into regime II, branches with comparable lengths are observed. For branched structures, the length is interpreted as the distance between two adjacent nodes or a node and an adjacent end-cap. Hence, there is an overall reduction in 〈L〉 as the morphology changes from mostly linear structures with short branches to highly branched structures (note that the total number of surfactants is kept constant). Finally, in region III, the increased entropy of free Sal⁻ ions in the solution favors more end-caps and 〈L〉 decreases again.

We probe the mechanisms that cause the anomalous variations of η₀ with respect to R shown in Fig. 13(a). We extract η₀ from reverse non-equilibrium molecular dynamics (RNEMD) simulations.⁴⁰ Briefly, by imposing a momentum flux in the system, we measure the shear strain as the response. This is different from the conventional method in which one imposes shear strain, and the resulting stress is calculated. For a moderate perturbation, which is controlled by varying the frequency of swaps and the number of swapped atoms, we identify the region where the velocity profile across the simulation box is linear. This region is found to be 0.04 ≤ Wi ≤ 0.40, where $Wi= γ ̇ τ$ is the Weissenberg number defined as the product of shear rate and a characteristic stress relaxation time τ which is calculated from the shear stress autocorrelation function for R = 0.4 at the lowest shear rate (Wi = 0.04). This gives τ ≈ 3 ns. Figure 14 shows the velocity profile and the shear viscosity at different Wi for c_s = c_D = 0.34M and illustrates that the viscosity remains approximately constant for small deformations. We emphasize that this particular choice of τ does not affect the final results for the zero-shear viscosity η₀ since it is calculated by extrapolating the shear viscosity at low shear rates $η ( γ ̇ → 0 )$⁠.

We compute the shear viscosity for c_D = 0.34M at various salt concentrations (0.4 ≤ R < 1.8) and extrapolate η₀. Interestingly, η₀ exhibits two maxima as observed in experiments for similar systems.^31–33 To further understand the underlying structural changes, we show the average micelle charge Q, fraction of unbound counterions $c s f / c s$⁠, and the node density N_n in Fig. 13(b). We can correlate these data and micelle configurations with the viscosity changes. Initial increase in viscosity and 〈L〉 in region I can be attributed to the transition from spherical or short cylindrical to wormlike micelles. The first maximum in viscosity occurs approximately at R = 0.7, in quantitative agreement with the experimental results for CTAC/NaSal and CPyCl/NaSal systems.^31–33 Beyond this maximum, the non-uniform distribution of counterions over a micelle surface (similar to the one shown in Fig. 12, where we have shown a non-uniform distribution of counterions on the cylindrical and spherical parts of a micelle pictorially) may create regions with unequal energy distribution. Regions with smaller counterion surface density penalize more energy like the end-caps of a rodlike micelle. At these points on the micelle surface, the end-cap of another micelle may attach to minimize the energy cost leading to the formation of branched structures as hypothesized in Ref. 5. Consequently, as shown in Fig. 13(a), viscosity decreases for 0.8 < R < 1.0 due to a transition from linear to branched micelles as observed in experiments.^31–35 Closer inspection of the microstructure reveals that some of the micelles are simply branched (Y junctions) or have multiple branches while some are even cross-linked (X junctions). This fact is illustrated by an increase of N_n shown in Fig. 13(b) and corresponds to the point where the effective micelle charge (Q) changes its sign. Such microstructures have previously been reported for several surfactant systems.^2,5,42 In fact, these structures are less viscous than entangled networks⁵ of WLMs and offers a faster mechanism for stress relaxation by sliding the cross-links along the contour. Beyond the viscosity minimum at R ≈ 1.0, micelle branches start to merge forming loops, resulting in a multiconnected network as shown in Fig. 15 (Multimedia view) of approximate mesh size ξ ≈ 21 nm. The viscosity increases upon network formation and passes through a second maximum at R ≈ 1.5 when the network is fully saturated as signified by a maximum in N_n in Fig. 13(b). The increase in viscosity from the branched state to the interconnected one is likely due to the increased rigidity of the network. The viscosity changes from 0.010 Pa s at the minimum to 0.016 Pa s at the second maximum, i.e., by a factor of 1.6. This is comparable to experimental observations for the CTAB/NaSal system.³⁴ It has been conjectured in the literature that micelle solutions with branched structures are less viscous as the branches provide an alternative (sliding) mechanism for stress relaxation.⁵ For multiconnected structures, there are two competing factors that determine the solution viscosity: network formation tends to increase the viscosity, while sliding of the junctions along the micelle contour tends to reduce it. The resultant change in viscosity is hence moderate. Further, our system sizes are orders of magnitude smaller as compared to experimental ones and finite size effects will be reflected in the MD predictions. In region III, the fraction of unbound counterions increases nonlinearly⁴¹ with increasing R as shown by open circles in Fig. 13(b), and the increased electrostatic attraction between the free counterions and the bound surfactants results in a gradual disintegration of the network as evidenced by a decrease in N_n. Consequently, micelles with branches or linear micelles are formed again. Hence, our simulations show that the following sequence of morphological transitions manifests as the double maxima observed for η₀ vs. c_s in Refs. 31–34: spherical → wormlike (η₀↑) → branched (η₀↓) → multiconnected (η₀↑) → branched/wormlike(η₀↓). Figure 16 shows typical micelle structures at different values of R. It is important to note that these are representative shapes obtained from our simulations. The solution itself consists of a collection of such microstructures or a distribution of different shapes.

In summary, we have performed large-scale MD simulations to study the topology, kinetics, and energetics of a model cationic surfactant micelle solution over time scales sufficiently long enough to probe micelle recombination dynamics. We have numerically computed the principal curvatures of the micelle-water interface and estimated the surfactant packing parameter as a function of salt concentration. In agreement with the geometric considerations, we have shown that the packing parameter changes from 1/3 to 1/2 upon increasing the concentration of highly binding aromatic counterions resulting in the sphere to rod transition of micelles. The simulations have unraveled a variety of self-assembled micelle structures such as spheres, cylinders, flexible worm-like, branched, networks, and handles. We have developed a novel contour tracking algorithm to quantify the topological features of linear, branched, and multiconnected structures. We have computed explicitly, for the first time from molecular simulations, the microstructural features such as contour length distribution, persistence length, mesh size, and branching density. Further, we have predicted reliably the end-cap energy from the dependence of the contour length on the solution temperature and estimated micelle recombination time from direct visualization of micelle coalescence events. Finally, we have correlated the growth and branching of micelles with the solution viscosity and rendered direct evidence to speculated mechanisms underlying the experimentally observed anomalous variations in viscosity with salt concentration. Specifically, electrostatic screening between the positively charged surfactant head groups by the negatively charged counterions adsorbed onto the micelle-water interface leads to an increase in the micelle length and hence the solution viscosity. At higher salt concentrations, micelle branching causes the solution viscosity to decrease. However, further addition of salt favors the formation of networks with increased viscosity until the micelle surface is saturated with adsorbed salt. Beyond this point, excess free counterions in the solution interact strongly with the oppositely charged surfactant head groups and the network disintegrates resulting in a less viscous solution. Given that the model captures experimentally observed trends, the simulation framework presented here can offer an ideal platform to probe non-equilibrium phenomena in ubiquitous micellar fluids.

A non-monotonic variation of viscosity with salt concentration Fig. 13 results from the modulation of the electrostatic interactions in the solution by strongly binding Sal⁻ ions. As such, the surface charge of the micelles is crucial in determining their morphology. Specifically, the onset of micelle branching is accompanied by a surface charge reversal due the condensation of Sal⁻ ions onto the micelle surface. Therefore, it is interesting to explore the effect of counterion valency on micelle morphology and solution rheology. Based on the insights gained from the current work, one might conjecture that for multivalent counterions, the onset of branching could occur at much lower values of the salt to surfactant ratio.

This work used the computational resources provided by the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation Grant No. OCI-1053575. The authors acknowledged financial support by National Science Foundation under Grant Nos. 1049489 and 1049454.

In this  Appendix, we present the details of simulated systems, molecular dynamics trajectories, and an algorithm to track the contour length of floppy micelles.

We have enclosed the following animations:

Figure 2  shows the self-assembly (SA) of surfactants in aqueous solution starting from an initial random configuration.

Figure 7  illustrates the contour tracking algorithm.

Figure 8  shows the formation of a toroidal micelle from a floppy cylindrical micelle.

Figure 9  shows a micelle surface and the constructed Monge patch.

Figure 11  illustrates the recombination kinetics of a tagged micelle in MD simulation (only type C₁ and C₂ beads are shown for clarity).

Figure 15  shows three-dimensional views of a multiconnected network. Only type C₁ beads (see Fig. 1) of CTA⁺ are shown.

The results presented in this paper are collected from trajectories that span >0.5 μs in time. We simulate the systems described below.

System I

Number of CTAC molecules N = 800, number of water beads N_w = 85 000 at various salt concentration (0 < R < 2). Shape parameter in Fig. 5(b) in the main text is extracted from these simulations. The final simulation box without salt in this case is approximately 23.1 nm, which corresponds to c_D ≈ 0.1M.

System II

N = 16 000, N_w = 1 040 000 at various salt concentrations (R = 0.2, 0.8, 1.2 and 1.6). The final simulation box at the lowest c_s in this case is approximately 54 nm, which corresponds to c_D = 0.16M. The persistence length (l_p) in Fig. 5(c), length distributions in Fig. 10, and the end-cap energy (E_c) in Fig. 12 are calculated from these simulations.

System III

N = 4000, N_w = 120 000 at various salt concentrations (0 < R < 2). The non-monotonic variation of the viscosity and micelles length scaling presented in Fig. 13 in the main text is computed from these systems. The final simulation box at the lowest c_s in this case is approximately 27 nm, which corresponds to c_D = 0.34M solution.

The phase diagram in Fig. 5(a) is constructed from the above three sets and additional simulations conducted at intermediate concentrations.

To analyze various length scales and energetics, we developed algorithms for micelle tacking. We divide the microstructure into cluster of micelles from a nearest neighbor analysis. The algorithm is as follows:

• We first measure the nearest neighbor distance r_min between the hydrocarbon beads from the position of the first peak of the radial distribution function g(r).

• Two beads within a distance of d < r_min + δ, where δ = 1.0 in LJ units, belong to the same micelle. Search for the hydrocarbon beads belonging to the first micelle continues for N number of steps, and the following search excludes the previously assigned beads.

• After the first operation, we randomly pick a hydrocarbon bead that has not been previously assigned to a micelle as a starting point to identify a second micelle. The search continues until all possible beads are identified. Iterations (i) and (ii) continue until each hydrocarbon bead is assigned a unique identification number.

• We then collect the head groups, and the counterions belonging to each micelle by a simple distance analysis.

• Once we know the surfactants and the counterions belonging to a micelle, it is straightforward to calculate the effective charge (Q) of the micelle, and the fraction of the unbound counterions $c s f / c s$⁠.

Quantifying different length scales in micellar solution is a non-trivial problem. It is worth noting that Dijkstra’s algorithm⁴⁵ is not applicable to our case since the end-caps of micelles are not a priori known, because the flexible WLMs are in a state of dynamic equilibrium. Therefore, we develop an algorithm to track the direction of maximum density of the CG beads to construct the contour path of the wormlike chains. The algorithm is illustrated in Fig. 17 and consists of the following steps.

• First, we calculate the number of CG beads enclosed in a blob of radius l around each CG bead. The position of a bead corresponding to the largest number density is considered as the starting point as shown by the dark yellow bead in Fig. 17(a).

• In the following step, we place test blobs of radius equal to that of a spherical micelle around the original blob. We then compute the number density of beads in each test blob, and the blob with the maximum density is a new search position: red balls in Fig. 17(a). The iteration continues until 90% of the beads are enclosed in those blobs.

• Once we determine the contour path, we place more blobs between neighboring blobs and relax the system by employing the angle and bond potential between the neighboring beads to keep the original orientation. The accuracy of the distance measurement in these calculations is approximately 0.4 nm.