Amino acid-based surfactant SLSar (30 wt%) containing 0.2 wt% NaCl (Crodasinic™ LS30) was used along with CAPB (30 wt%) and 4.2 wt% NaCl (Crodateric™ CAB 30). The agents used to adjust pH were NaOH and citric acid, from Sigma Aldrich, both above 99% purity. Distilled water was used in all experiments.

Solutions with eight different mass ratios of CAPB/SLSar were prepared as detailed in Table 1. Samples were mixed using a digital magnetic stirrer until total homogenisation. The solutions varied from pH 3.5 to 6.5. After preparation, the solutions were allowed to equilibrate for at least 24 h before any experiments were performed to ensure equilibration.

To explore the rheological properties of solutions, measurements were taken using an Anton Paar MCR 302 rotational rheometer. Measurements were performed with a cup and bob stainless steel device of 18 mm diameter and 17 mm diameter, 25 mm length, respectively, with a fixed gap. A solvent trap was used to prevent solvent evaporation during measurements. Samples were allowed to relax and acclimatise for 5 min after loading, before measurement. All measurements were carried out at 25°C in duplicate and the average was reported in the results.

The flow properties of the samples were measured by recording shear stress (σ) and viscosity (η) values when shearing the samples at ascending and descending shear rates ranging from 0.01 to 1,000 s⁻¹ then 1,000 to 0.01 s⁻¹ (logarithmic scale). Data was gathered at 10 points per decade. The parameter collected during this experiment was the zero-shear viscosity, denoted η ₀ and oftentimes obtained in the quasi-Newtonian (plateau) region.

Table 2 shows the results of the flow shear rate sweep tests. In regions that were both low pH and low mass ratio of CAPB/SLSar turbid systems were observed. These liquids were white and cloudy and, under cross polarised microscopy, showed an absence of crystalline material. Shortly after blending these liquids had extremely high zero-shear viscosities of tens of thousands of mPa s. After a few hours these liquids transitioned to low viscosity. It is likely that this behavior is related to the protonation and precipitation of the acid form of SLSar at a pH close to or below its pK _a of 4.5 (Wallach et al., 1992). Similar effects were encountered by Vu et al. (2020) and we agree with their assessment that these systems are best discarded due to their unsuitability for formulation. The remaining systems listed in Table 2 were clear and colorless liquids with differing viscosity as reported.

We observe a steadily increasing viscosity peak as the mass ratio of CAPB/SLSar is increased up to 1.5. Further increasing to the ratio to 4.0 results in a dramatic drop in recorded viscosity. We observe the largest peaks in viscosity at close to the experimentally recorded value pK _a for SLSar of 4.5 where the molar ratio of anionic to acid form of SLSar is approximately 1:1. Similar behavior was observed by Vu et al. (2020) who hypothesized that the hydrogen bonding between the acid and anionic forms of SLSar and the zwitterionic component (CAPB here) at the effective pK _a is a contributing interactions leading to viscosity enhancement.

In our DPD model (which follows the approach of Anderson et al. (2017)) both surfactant molecules are defined by a set of linearly bonded beads comprised of multiple chemical groups. These beads contain 1–5 “heavy atoms” (i. e. carbon, oxygen, nitrogen, sodium). For this study we use eight surfactant bead types; three corresponding to alkyl groups [CH₃CH₂] [CH₂CH₂] and [CH₂]; a secondary amide [CH₂ − C (=O)NH − CH₂]; a tertiary amide [CH₂ − C (=O)N(CH₃) − CH₂]; a quaternary ammonium cation [ CH 2 N + ( CH 3 ) 2 CH 2 ]; and the deprotonated carboxylate [COO⁻] or protonated carboxylic acid [COOH] alternatives. These beads should of course be bonded as appropriate to the chemical stoichiometry of the molecules. Three supramolecular solvent beads are also used: [2H₂O] [Na⁺ ⋅2H₂O], and [OH⁻⋅H₂O], representing water, sodium, and hydroxyl ions respectively. These latter beads are all the same in the model apart from the charge, in other words free ions in this approach are represented by charged water beads. In developing this model we follow the approach outlined in prior works (Anderson et al., 2017; Anderson et al., 2018; Anderson et al., 2023; Panoukidou et al., 2019; Del Regno et al., 2021; Bray et al., 2020; Bray et al., 2022; Wand et al., 2020), where the bead volumes are based on partial molar volumes (Durchschlag and Zipper, 1994) and the interactions between DPD beads upon log  P and density measurements. As this method is now well established we do not give more extensive details here, rather we point to the Supplementary Material. Following this model, SLSar in the deprotonated state takes the structure.

[CH₃CH₂]–[CH₂CH₂]₄–[CH₂ − C (=O)N(CH₃) − CH₂]–[COO⁻], with [Na⁺ ⋅ 2 H₂O] as counterion. The deprotonated form of the zwitterionic CAPB follows similarly,

except that no counterion is needed because in this state the molecule is charge neutral. In the case that either surfactant is protonated the terminal carboxylate bead is replaced by a carboxylic acid bead, and in the case of CAPB a compensating number of [OH⁻⋅ H₂O] counterions are added to maintain overall charge neutrality. The surfactant structures are illustrated in Figure 1.

We calculate the CMC of the surfactants (following Anderson et al., 2018) to validate our model. For SLSar we obtain a calculated CMC of ≈ 10 –14 mM, versus an experimental value of ≈ 12 –14 mM (Patra et al., 2018). Likewise, for CAPB we find a CMC of ≈ 2 –4 mM, versus an experimental value ≈ 3 m M (Dai et al., 2014; El-Dossoki et al., 2020).

For a given pH both CAPB and SLSar surfactants can co-exist in their given acid and conjugate base form, i. e., in carboxylic acid (protonated)/carboxylate (deprotonated) form, respectively. If the pH is decreased, the fraction of protonated surfactants increases as determined by the dissociation constant pK _a of the carboxylic acid group. Changing the ratio of protonated to deprotonated surfactant can have dramatic changes on the self-assembled structures and therefore the properties of formulations (Vu et al., 2021a; Vu et al., 2021b; Xu et al., 2019). We represent pH in our simulations by adjusting the ratio of protonated to deprotonated surfactants following the classic Henderson-Hasselbalch model in Eq. 1 (Henderson, 1908). From the reported pK _a of the surfactants, i. e. pK _a = 4.5 for SLSar (Wallach et al., 1992) and pK _a = 1.83 for CAPB (Wood, 1987), the ratio of protonated and deprotonated surfactants at each pH is calculated from log 10 [ AH ] / [ A − ] = p K a − p H . (1) Here [AH] and [A⁻] are the concentrations of the protonated and deprotonated surfactant, respectively. We assume that both surfactants are monoprotic acids (i. e. carboxyl group) and ignore the highly basic ammonium ion of CAPB which has a pK _a ≳ 9, since we do not work in such high pH conditions.

In reality the degree of protonation depends upon the surroundings of the surfactant, i. e. the effect of self-aggregation on the acid-base equilibrium of one of these surfactants is a known effect resulting in a p K m (apparent pK _a for micellar surfactant) shifted from the aqueous counterpart by factors such as hydrogen bond availability, dielectric constant changes, surfactant distribution in the micelles, and overall surfactant concentration (Maeda and Kakehashi, 2000; Goldsipe and Blankschtein, 2006). Vu et al. (2020) determined the effect of the co-surfactant cocoamidopropyl hydroxysultaine (CAHS) on SLSar by potentiometric titration, observing a shift from pK _a ≈ 5 to 4.5 when going from pure SLSar to a CAHS/SLSar mass ratio of 60:40. Whilst we recognise that a shift can occur we adopt pK _m = pK _a in this work as the offset may not be known for de novo surfactants in future studies, and we wish to trial the simplest approach possible.

Figure 2 shows the shift in populations of acid and conjugate base as the pH changes. We can see over the pH range studied (i. e. between 3.5 and 6.5 in our study) that CAPB is in near-complete acid form and changes in populations of the SLSar forms is the main driving force behind any pH effect observed (when using fixed CAPB/SLSar surfactant mass ratio).

When building our models, we do not consider water dissociation. The concentrations of OH⁻ and H₃O⁺ ions from such are negligible compared to the counterion concentrations unless the solution becomes very acidic with a pH below 3. The presence of counterions in solvent (corresponding to [Na⁺ ⋅2H₂O] or [OH⁻⋅H₂O]) is used to maintain charge neutrality.

We sample the same eight CAPB to SLSar mass ratios as covered experimentally but chose to explore a reduced number of pH levels covering the same range–this was done to provide good sampling whilst reducing computation cost. Simulations were performed at seven pH values equivalent to 3.9, 4.3, 4.7, 5.1, 5.5, 6.0 and 6.5. As described above, the pH was used to fix the ratio of protonated to deprotonated surfactant, as a “quenched” variable. For each mass ratio and pH, we ran simulations in box sizes of 40³ (22.6 nm³) to ensure enough material of each form of surfactant was present in the box and that large (relatively) aggregates have the space to form. We discuss the consequence of altering the box size in Section 4. For each system, we ran an initial simulation to determine a resultant N _agg value. Here, we adopt a cut-off of N _mic = 6, i. e., aggregates with more than 6 surfactant molecules present are considered as micelles. Anything below this point is considered a submicellar aggregate. We base this on previously reported works (Anderson et al., 2018). Where the final time-step N _agg from the simulation was greater than seventy-five, a total of ten repeats simulation of the CAPB/SLSar mass ratio/pH combination were conducted, each with a different random seed. This was done to obtain average, reliable, details of the self-assembly process Where N _agg < 75 was observed, we did not perform further repeats as random sampling of systems fulfilling this criterion showed little variation (less than ±5). This is discussed in more detail in Section 4.

For the simulation time adopted, we choose to go against convention where the aim is usually to determine equilibrium structures, and here explore self-assembly and micelle behavior over shorter time-scales. As such we sample self-assembly over one million time-steps We do this for two reasons.

• The first is computational efficiency. We have demonstrated in our previous work (Anderson et al., 2023) that equilibration of micelle structures can take a very large simulation time (at least 10–40 million time-steps in the case presented in the cited literature). Whilst it is technically feasible to simulate such long time-scales given enough resources, to be of practical use it is of our view that the computational screening protocols ought to be more cost efficient and less time consuming than the alternative of running the laboratory trials.

Assessing the utility of short time-scale simulations (exploring early stage self-assembly) is core to our research question in performing this work. For selected systems we explore longer, four million time-step, simulations in order to explore the consequences of extended run times.

DPD simulations were performed using the dl_meso simulation package (Seaton et al., 2013) and we adopt a DPD time-step of 0.01 DPD time-units. Hereafter we shall, for the most part, express elapsed times in DPD time-units. One DPD time-unit corresponds to ≈ 50 p s of real elapsed time, according to previously reported mappings based on matching diffusion coefficients (Fraaije et al., 2018; Anderson et al., 2023). Hence the adopted DPD time-step of 0.01 DPD time-units corresponds to ≈ 0.5 p s . This is a couple of orders larger than the typical femtosecond time-step in MD, which goes some way to explain the preference for using CG methods such as DPD for simulating surfactant self-assembly.

A constant pressure ensemble (NPT), following the Langevin piston implementation by Jakobsen (2005) was used, and we set the system pressure to p = 23.7 (DPD units) to match the pressure of pure water in the model defined by Groot and Warren, (1997). The electrostatics of the charged DPD beads were defined by the model of González-Melchor et al. (2006) which assumes a uniform dielectric constant (chosen to be equivalent to pure water by setting the electrostatic coupling parameter Γ = 15.4) in the simulation box. The electrostatics were solved using a smoothed-particle mesh Ewald (SPME) algorithm (Essmann et al., 1995). All simulations were started from randomly dispersed initial configurations.

Simulation trajectory files were analyzed using the ummap package (Bray et al., 2019), which was modified for this work to enable determination of branching in micelle structures following the method presented in Conchuir et al. (2020) (see below). The trajectory data comprise simulation snapshots taken every 2000 DPD time-steps (i. e. 20 DPD time-units ≈1 ns). The N _agg (number average) values are extracted as previously described in Anderson et al. (2018) after partitioning the surfactants into sub-micellar and micellar aggregates. N _agg was plotted as a function of time to understand self-assembly behavior. Shape metrics determined to understand the average morphology of micelles were calculated as reported in the original ummap paper using cut-off values ϵ = 0.333, ϵ _rod = 0.5, ϵ _disc = 1.0, N _rod = N _disc = 100. These values were determined by trial and error to ensure that observed structures (i. e. from visualizing simulation outputs) corresponded to the shapes described through automatic identification using the ummap package. In this work we combine spherical, oblate and prolate spheroids resulting from the shape analysis into one “spheroid” class as the distinction between these structures does not not provide detail relating to viscosity (all structures in the spheroid class would result in low viscosity). Micelles with gyration tensors that reported as rod or ellipsoid by ummap are classed as “rod-like”. This was decided upon after visualising the resulting ellipsoids resulting from the gyration tensor which were found in reality to be exclusively “bent” rods. All other classifications remain the same as the original citation except lamellar which refers to a structure with two principle extents (defined using the gyration tensor) longer than the bounding box size. In this work, these structures are generally branched micelles (determined via visual inspection) and we refer to these as 2D extended. Those with three extents larger than the box dimensions are referred to as 3D extended. It should be noted that for aggregates that self-interact across the boundaries only the unreplicated structure (suitably unwrapped) is used for shape analysis (and aggregation size) which limits their total extents. Hence care must be taken in interpreting the 1D to 3D extended structures where replicas can be staggered and stacked within the box. For these a clearer interpretation is given using the below “skeleton” method which also indicates connections between the replicated structures.

The degree of branching in micelles was determined by identifying the “skeleton” (central path) of each micelle, using a modified version of the Conchuir et al. (2020) method, and locating the points at which a junction occurred in the resulting micelle skeleton graph. From this method we can identify both micelle branches and end caps. Examples of the resulting skeleton for branched micelles are presented in Figure 3 where the skeleton is shown running a central path through the micelle. Each “red” ball is a segmented fragment of the micelle with the “blue” ball providing additional links obtained by studying the connections between segments. The line then provides the skeleton connections. The branch points are located at red points where three or more lines connect and the end caps where only one line connects. A finite branched WLM (2 branch points, 4 end-caps, no closed paths) is seen in mass ratio 0.67 at pH3.9 from a 40³ simulation box (Figure 3A). Branched micelles that wrap around the box and self-link to form an infinite length (e. g. 2D and 3D extended structures) have more branch points than micelle end-caps as did cases where closed paths formed and each junction is treated as a separate branch point. Examples are Figure 3B where we see a heavily interconnected structure (7 branch points, 3 end-caps, several closed paths) from a single simulation snapshot from the mass ratio 0.0 at pH3.9 in a 40³ simulation cell. In Figure 3C we show a heavily interconnected (13 branch points, 7 end-caps, many closed paths) structure resulting from a simulation snapshot from the system with mass ratio 0.0 at pH3.9 in a larger, 50³, simulation cell.

In the text we report the average number of branch points per micelle ⟨N _branch⟩ as indicative of the overall system. However, in practice some resulting micelles were highly branched whilst several remaining shorter and rod-like micelles which lower the reported average value.

In Figure 4 we show the aggregation behavior as a function of time for the CAPB/SLSar mass ratio 0.67 mixture at four different pH values (Figures showing N _agg evolution for all sampled ratios and pH levels are supplied in the Supplementary Material). Figures 4A–D shows the behavior at pH 3.9, 4.3, 4.7 and 5.1, respectively. The plots show the average aggregation number (⟨N _agg⟩ hereafter) (black line), resulting from the ten repeat simulations, in addition to the best fit (red dashed line) obtained by fitting to ⟨N _agg⟩. In each case, the fit was determined with a simple nonlinear least squares procedure (via Xmgrace). Calculated fits show very high correlation with the determined value of ⟨N _agg⟩. By using replicate simulations we have been able to determine three distinct growth profiles for self-assembling micelle aggregates. These are logarithmic, linear and cubic (referred to as “cubic” hereafter) micellar growth. These behaviors are discussed more throughout the article. Results from the ten individual simulations are shown as colored dotted lines.

For the above system ⟨N _agg⟩ determined at the final time frame decreases with increasing pH from in excess of 400 at pH 3.9 to less than 80 at pH5.1. This trend is observed for each CAPB/SLSar mass ratio although the effect is more pronounced for the SLSar rich mixtures as this surfactant has a significantly stronger response to pH (as shown in Figure 2). Each of the presented four systems undergoes an early stage logarithmic growth region which persists to a N _agg of ≈ 70 –80. Following this, for systems that increase further in N _agg we observe a linear growth region, the gradient of which is correlated with pH. The lower the pH, the larger the gradient, and therefore the larger the resultant micelles at the end of the simulation. For the system sampled at pH 3.9, we observe a cubic growth region that commences at ≈ 7000 DPD time units that is not present for the other three pH values sampled for the CAPB/SLSar 0.67 mass ratio. We note that the appearance of the cubic growth regime seems to correspond with the formation of branches in (and between) micelle structures (see shape and branching metrics in Figure 5). Comparison of Figure 5B with Figure 4B shows that very small amounts of branching (leading to small amounts of extended structures as defined by shape metrics) does not prevent ⟨N _agg⟩ as being fitted to a linear growth relationship.

Inspection of the colored (dotted) lines in Figure 4, corresponding to individual simulations, shows the extent of variation in individual simulations. In the early, logarithmic growth regions, most simulations report aggregation behavior that is similar across repeats (Figure 4D). As systems present with strong linear or cubic growth, larger differences can occur. For example, the system corresponding to pH 3.9 (Figure 4A) has a discrepancy of ≈ 500 in the reported final frame N _agg value between the maximum and minimum determined values. For the same CAPB/SLSar ratio at pH 5.1 (Figure 4D) a discrepancy of ≈ 10 is observed.

In Figure 5 we present shape and branching metrics corresponding to the micelles formed in the CAPB/SLSar mass ratio 0.67 simulations (the same systems as Figure 4 for consistency). Additional figures showing shape metric evolution over time for all sampled ratios at pH levels from 3.9–4.7 are supplied in the Supplementary Material. Note that the presented images show the average shape and branching behaviors across the ten replicate simulations. For a CAPB/SLSar mass ratio 0.67 at pH3.9 we see an early majority of spheroid micelles decrease rapidly with time as both rod-like and then, slightly later, extended structures such as WLMs and 2D branched morphologies appear. This appearance of extended structures corresponds to the formation of branching in the simulations. The initial increase in rod-like micelles is reversed as the fraction hits ≈ 0.6 , at which point extended structures begin to develop. This occurs as rods extend into WLMs, with rods extending across periodic boundaries in predominantly one direction. More significant branching is then observed as WLMs form branches that shift the extended structures across periodic boundaries in multiple dimensions. This pattern of behavior is repeated for mass ratio 0.67 at pH4.3 (Figure 5B) where there is a slower decrease in spheroid micelles, a corresponding slower increase in the rod-like and WLM micelles, and hardly any 2D extended structures appearing (a very small degree of 2D extended shape is present in the final reported time frames corresponding to a small degree of observed micelle branching). Again, the onset of extended structures (WLMs in this case) appears as the fraction of rod-like micelles crosses ≈ 0.6 . This is again repeated in Figure 5C with the rates being slower once again with almost no WLMs (or branching) being observed. In the case of the mass ratio 0.67 at pH5.1 system the dominant shape remains the spheroid with a minority of rod-like micelles forming. No WLMs or extended structures appear in this case. Note that the rod-like fraction of all the systems presented display an early peak and then steep decay in the very early time-scales. Through visualisation of the morphology produced we attribute to the behavior of very small micelle aggregates just above the cut-off of N _mic = 6.

Whilst we choose to explore the early time self-assembly and aggregation behavior of surfactants in this study, it is natural to ask the question what happens if the simulations are extended beyond this time window. To answer this we have sampled several systems over a longer time period of 40,000 DPD time-units. Examples of the longer time aggregation behavior are presented in Figure 6. In these images ⟨N _agg⟩ is displayed for the full 40,000 DPD time-units along with best fit lines calculated over the initial 10,000 DPD time-units (red dashed lines). The systems studied in Figures 6B, C (corresponding to mass ratio 0.25 at pH4.3 and mass ratio 0.67 at pH5.1, respectively) remain broadly consistent with the behavior observed over the first 10,000 DPD time-units (i. e. dominated by logarithmic growth for mass ratio 0.67 at pH5.1 and linear growth for mass ratio 0.25 at pH4.3). As the simulation for the mass ratio 0.25 at pH3.9 system (Figure 6A) is extended beyond 10,000 DPD time-units (where ⟨N _branch⟩ = 0.5), initially the observed cubic relationship in ⟨N _agg⟩ continues up until approximately 20,000 DPD time-units. At this point we observe a sharp deviation from cubic behavior as ⟨N _agg⟩ begins to plateau. This plateau corresponds to ⟨N _agg⟩ ≈ 3,000 which is almost the entire surfactant content of the simulation box. The plateau is therefore a result of the finite material contained within the simulation box becoming depleted. At the final time-unit for this system ⟨N _branch⟩ ≈ 8 and the resulting morphology is heavily interconnected and continuous. As we extend the mass ratio 1.22 at pH3.9 system beyond the initial 10,000 DPD time-units we observe a change from the initial linear behavior (with ⟨N _branch⟩ = 0.1) to a much sharper cubic growth regime with ⟨N _branch⟩ ≈ 2.5 at final time-unit. Figure 6D shows this change which can be compared to the projection of the original linear growth (red dashes).

In addition to exploring the effect that longer simulation times can have on the behavior of the observed aggregates, it is natural to consider what would happen if the size of the simulation cell is modified beyond the default 40³ explored in this work. In Figures 7A–D we demonstrate the impact box size can have on ⟨N _agg⟩ for mass ratio 0.0 at pH3.9, mass ratio 0.67 at pH3.9, mass ratio 1.22 at pH4.3, and mass ratio 1.33 at pH3.9. For the first two, we explore the box size effect over 10,000 DPD time-units. For the latter we explore the behavior 40,000 DPD time-units.

For the mass ratio 0.0 at pH3.9 system (Figure 7A), simulation in both a 40³ and 50³ boxes result in cubic growth regimes at the end of the simulation. Here, the transition from a linear to cubic regime appears to occur at approximately the same region (marked by an arrow on the image). Whilst there is some difference between the reported ⟨N _agg⟩ over time, the results are broadly within error of each other. In Figures 8A, B we look at the shape and branching metrics for these two box sizes for the mass ratio 0.0 at pH3.9 system. For both, the general morphologies reported are similar but there is a key difference in the reported behavior. The behavior of spheroid micelles is broadly the same in both box sizes. Rods can be observed transitioning into WLMs after a slightly longer period of time in the larger box. This is intuitive as the definition of WLM is set by the box size. As the box gets larger, rods can grow longer before their classification is updated. What is perhaps less intuitive is the behavior of the 2D and 3D extended structures. Larger, 3D extended, structures seem to be preferred over the 2D alternatives as the simulation box gets larger. We speculate that the origins of this behavior stem from spatial requirements to form the 3D extended systems. The increased proportion of 3D extended structures reported from the shape metrics in Figure 8B is accompanied by a significant increase in the degree of branching in the system and the larger box presents almost double the average number of branch points per micelles when compared to the smaller box. Examples of morphologies resulting from the two box sizes for the mass ratio 0.0 at pH3.9 system are presented in Figures 3B, C. From the consideration of the mass ratio 0.0 at pH3.9 and mass ratio 0.67 at pH3.9 systems, we conclude that smaller boxes promote the 2D extended structures, whereas the larger boxes promote the 3D extended structures.

In the mass ratio 0.67 at pH3.9 system, we observe a micelle growth region following cubic growth behavior for the 40³ system (Figure 7B). Upon enlarging the box to 50³ the cubic behavior disappears, leaving a system that appears to grow linearly. Comparing Figures 8C, 5A, which explore the shape metrics for the mass ratio 0.67 at pH3.9 systems in 50³ and 40³ boxes, respectively, shows that the smaller of the two boxes presents with more WLMs and fewer rods, as can be expected based on previous discussion. For the larger box, the onset of 2D extended structures, as determined by the shape metrics, is later than that for the smaller box. Both systems present similar degrees of branching although there appears to be more branching earlier for the larger box. For both box sizes, the total fraction of 2D extended structures is comparable at the final time-unit. Deviation from the linear behavior of the 40³ box occurs at the point at which the 2D extended structures appear in the analysis of the shape metrics ( ≈ 6 000 DPD time-units). We postulate that the deviation is much smaller in Figure 8C enabling the behavior to be classed as linear for the entire duration. It is possible that larger boxes require more/larger degrees of branching to register significant deviation from linear behavior on the time-scales sampled here. This requires further testing.

Similar size effects to those already observed can be seen in the behavior from the mass ratio 1.22 at pH4.3 system (Figure 7C) in which we explore 20³, 30³ and 40³ simulation cells. Reducing the box size from 40³ (black line) to 30³ (violet line) results in similar behavior with both box sizes following the same best fit line (linear). Further reduction to a box size of 20³ results in different behavior, however. Here we see a micelle growth behavior that turns to a cubic just prior to 10,000 DPD time-units (blue line). After a short cubic region (demonstrated by the orange dashed line) we notice that the ⟨N _agg⟩ begins to plateau as was observed in Figure 6A. The resulting aggregate size for this box size corresponds to the total number of surfactants present in the simulation box. For the mass ratio 1.22 at pH3.9 system (Figure 7D) we appear to observe the same behavior on reducing the box size as surfactant is depleted in the smaller box.

The exploration of simulation times and box sizes presented here provides some degree of the behaviors that may be observed as changes are made to the simulation protocols. One could argue that longer simulations are required, or that larger boxes are necessary. The difficulty is where does this end? It is always possible (theoretically) to run longer/larger, and new behaviors may be observed. Practically, however, in designing a simulation campaign, one must balance simulation cost (resulting from box size and simulation length) with the desired accuracy and reliability in the derived result.

Finally, in Table 3 we present the ⟨N _agg⟩ value, the average number of micelles and the average number of micelle branch points, in addition to the overall maximum and minimum micelle sizes recorded over the ten repeats (where performed), for each of the sampled systems. We also present the functional form (logarithmic, linear or cubic) of a nonlinear least squares fit performed on the average micelle growth rate. We note that the different, dominant, growth behaviors are a continuum and therefore we assign the category based on somewhat arbitrary cut-off values. We are able to observe, in Figure 5, that cubic growth is a result of branching and the corresponding formation of extended structures. However, as we see from comparing Figure 5B with Figure 4B small amounts of branching can be tolerated in what is defined as linear growth.

From inspection of Table 3, we observe that low CAPB/SLSar mass ratio and pH (top left of Table 3) results in the largest micelle growth. This appears intuitive as here there is less charge repulsion in the resultant micelles as most surfactants appear in a protonated state, potentially leading to decreased curvature and therefore larger structures. In this region, for pH 3.9, strong cubic growth is observed, brought about by micelle branching and the formation of extended structures. As the pH increases more the table becomes dominated by linear growth (corresponding to rod-like species and much less branch formation), until at even higher pH values the table presents extensive logarithmic growth of smaller spheroid type micelles. Broadly, resulting ⟨N _agg⟩ decreases with increased CAPB content (left to right in the table) and decreases as pH increases (top to bottom in the table). However, at the bottom of the table ⟨N _agg⟩ is lower for SLSar rich systems. The decrease in size with increasing pH is more pronounced for the high SLSar content systems, which is a consequence of the larger changes to the protonation state of SLSar in response to pH versus that of CAPB. At the extreme, pure CAPB shows little change in behavior with increasing pH as the protonation state of this surfactant remains largely unchanged.

Table 3 highlights the difference in the maximum and minimum sizes obtained in the simulations and shows the danger of not performing several repeats in a study like this where we are looking to correlate self-assembly of micelle aggregates to the potential for viscosity build, especially for systems with resultant ⟨N _agg⟩ above 100. In the extreme case, if one was to perform only a single simulation, it would be very hard to make sense of the results and trends in behavior would not emerge cleanly. For systems that tend to form the largest aggregates, we see significant discrepancies between the maximum and minimum reported values at the final time-step of simulation.

We label certain systems ⟨N _agg⟩ in bold in the table. These are systems that we feel could be of interest experimentally based on two determinant cut-offs. The first is the resultant ⟨N _agg⟩ itself which we set as requiring to be ≳ 1.5× the typical spherical micelle size (enabling us to be well into the rod forming region at minimum), which also corresponds to the approximate ⟨N _agg⟩ where half of a system morphology is in spheroid configuration and half in rod-like configuration (See, for example Figure 5C). In this study, this value is ≈ ⟨ N agg ⟩ = 110 . The second criteria is that the system must be in a linear growth region as cubic growth regions could indicate large structures not conducive to appropriate viscosity build (See discussion in Section 4).