Dioctyl sodium sulfosuccinate surfactant self-assembly dependency of solvent hydrophilicity: a modelling study †

The self-assembly of dioctyl sodium sulfosuccinate (AOT) model surfactant in solvent environments of diﬀering polarity is examined by means of dissipative particle dynamics (DPD) bead model parametrized against Hildebrand solubility parameters from atomistic molecular dynamics (MD) simulations. The model predicts that in hydrophobic solvents ( e.g. dodecane) the surfactant forms small ( N agg B 8) reverse micellar aggregates, while in a solvent corresponding to water lamellar assembly takes place, in good agreement with literature structural parameters. Interestingly, solvents of intermediate polarity lead to formation of large, internally structured aggregates. In these, the surfactant headgroups cluster within the aggregate, surrounded by a continuous phase formed by the hydrocarbon tails. We show that the partitioning of the headgroups between the aggregate surface layer and the inner clustered phase depends primarily on solvent polarity, and can be controlled by the solvent, but also system composition. Finally, we compare the DPD assembly response to simplified eﬀective interaction potentials derived at dilute concentration limit for the interactions. The comparison reveals that the simplified eﬀective potential descriptions provide good level of insight on the assembly morphologies, despite drastic, isotropic interactions simplification involved.


Introduction
][9][10] Intermolecular interactions drive surfactant self-assembly into diverse morphologies, including spherical, applications in, e.g., nanoparticle synthesis, 10,29,30 as microreactors, and as extraction systems. 31,32On the other hand, in dry hydrocarbon solvents, AOT forms monodisperse, slightly aspherical aggregates with aggregation numbers corresponding to 20-30 AOT molecules, depending on the molar volume of the solvent. 14,24,33,34Additionally, aggregation propensity of AOT has been shown to depend strongly on the polarity of the solvent, as well as the presence of trace water in the system. 35,36Particularly, the shape and size of the reverse micelles depend on the waterto-surfactant ratio. 37,38][41][42][43] AOT-water self-assembling systems have received comparably less attention.5][46][47][48] The lamellar mesophase has been most well documented, 45,[48][49][50][51] however, bicontinuous cubic phases and inverse hexagonal phases have also been observed at high AOT concentrations. 33,45Additionally, small AOT micelles have been reported in experiments corresponding to isotropic solutions at low concentration, [52][53][54] although there remains a lack of consensus on micellar shape in these conditions.][57] Complementing the experimental characterization works reviewed above that provide information on the average AOT aggregate assembly structure, computer simulations, particularly molecular modelling techniques, are a powerful tool for determining the assembly morphology and dynamics of surfactant aggregates at a molecular level.Extensive computational studies of AOT self-assembly, with particular focus on AOT reverse micelles, 38,[58][59][60][61] exist in literature.The shape and fluctuations of small and larger preformed AOT reverse micelle aggregates, at dry and with various water loadings, have been extensively studied at both atomistic 13,38,58,[61][62][63] and coarse grained 59,60,[64][65][66] detail.Capturing larger scale self-assembly behaviour, such as the formation of AOT lamellar phases in water, is often beyond the time and length scale limitations of classical MD.On the other hand, atomistic MD simulations have been used to capture water and surfactant dynamics, and interface fluctuations in individual bilayers. 50,67Additionally, some work on overall bulk self assembly structures at atomistic detail level exist. 680][71] The trade-off is the abstraction of chemical features and molecular scale interactions, such as hydrogen bonding and solvent shell formation.3][74][75][76][77] Furthermore, the approach has also enabled establishing scaling laws. 78ere, we construct a coarse-grained (CG) DPD model for AOT surfactant and explicit hydrocarbon solvent.The model is based on a bottom-up parametrization approach which relies on atomistic MD simulation data.In the model construction, we target such degree of coarse-graining N CG that achievable assembly length and time scales become sufficient to examine large scale assembly changes.The simplified, CG surfactant can be considered to model AOT, because via systematic tuning of the model interaction parameters, the equilibrium assembly response of the AOT surfactant-solvent binary systems in aqueous and hydrocarbon solvents is achieved.We map the assembly response in varying solvent polarities with examined range covering solvents from an apolar hydrocarbon solvent, to polar water solvent, as well as intermediate solvent polarities.The resulting equilibrium self-assembly phases are characterised with comparison to existing experimental data to verify the model response.The self-assembly response is analyzed for assembly characteristics and system dependencies with the purpose of gaining solvent dependent control over the assembly response.In addition to the known apolar hydrocarbon solvent and water solvent responses, the characterization covers self-assembly phases in solvents on intermediate polarity, for which the AOT self-assembly response in largely unknown.Finally, to map the limits of model simplification, we extract the corresponding dilute limit effective interactions potentials for selected AOT systems and compare them and their finite concentration assembly response in the DPD simulations.

Dissipative particle dynamics
DPD, first formulated by Groot and Warren, 79 is a common CG particle-based simulation technique.The chemical system is described as interacting particles each corresponding to groups of atoms within a molecule or a region of solvent (i.e.cluster of solvent molecules).In addition to soft interaction potentials, DPD relies on dissipative forces calculated pairwise additively between the particles.The time evolution of the system of particles, commonly referred to as DPD beads, is obtained by integration of Newton's equations of motion.The DPD approach leads to conservation of momentum and preservation of hydrodynamics.The total force acting on a single DPD bead i 69,75 is formulated as In this, F C ij is the conservative force, F D ij the dissipative force, and F R ij the random force.Each of these forces is pair-wise additive, which means that the summation j a i runs over all neighbouring beads j within a certain cutoff radius r c .Hence, , where r ij measures the distance between beads i and j.
The repulsive soft-core potential, corresponding to the conservative force, allows using a significantly longer time step than in classical MD simulations.This translates into increased computational speed, as well as temporal and spatial resolution.The use of a soft-core potential is justified by the softness of effective interactions between groups of several hard-core (Lennard-Jones type) atoms. 80The DPD conservative force is given by where a ij is the conservative force repulsion coefficient between beads i and j.The unit vector corresponding to r ˆij = r ij /r ij is referred to as r ˆij .In this, r ij = r i À r j and r ij = |r ij |.The random F R ij and dissipative F D ij forces are formulated as respectively, where v ij is the relative velocity between beads i and j, and g and s are the friction coefficient and noise amplitude.The parameters o D (r ij ) and o R (r ij ) are the distancedependent weight functions for the dissipative force and random force, and x ij is a random number with zero average and unit variance (normalized Gaussian white noise).The dissipative and random forces of the DPD approach incorporate the effect of Brownian motion into the larger length scale and are responsible for momentum conservation in the system.Their coupling also forms the DPD thermostat via a pair-wise Brownian dashpot In this, T is the absolute temperature and k B the Boltzmann constant.Notably, the friction coefficient g influences viscosity and hydrodynamic modes dissipation in DPD; if modelling dynamical properties of the different solutions, g becomes an important parameter to consider in the model, see e.g.ref. 81-85.Here, the focus was on equilibrium structure, and a value of g = 4.5 was used. 79Assembly structure is mainly dependent on the conservative force parameters. 86In standard DPD, all interacting beads, regardless of bead type or model chemical species, have the same mass m i and diameter r c .Typically, and without loss of generality, DPD simulations employ reduced units, so that r c = m i = k B T = 1.

Surfactant model parametrization
In the parametrization, we focus on tuning the repulsion parameters a ij , which relate to the chemical nature of the modelled system and govern the equilibrium assembly response.Notably, if dynamical viscoelastic properties of the solutions are targeted, also the friction coefficient g and the fluctuations decay should be considered.The CG DPD representation of the AOT molecule consists of two hydrophobic tail beads (T) linked by a hydrophilic head bead (H), see Fig. 1.The solvent is described by a monomeric solvent bead (S).
The self-repulsion parameter a ii is typically obtained by matching to isothermal compressibility k T .For water at T = 298 K, with coarse-graining degree N CG = 1, and DPD density r = 3, a ii = 25 k B T. The same value is commonly used for also polymer blends, see e.g., ref. 79.To avoid losing generality of the model, we use the same compressibility value for the entirety of the examined range.a ii scales linearly with N CG , as does the bead diameter r c = (rN CG v m ) 1/3 , where v m is the molecular volume of -in most cases -water.Therefore, the repulsive potential becomes steeper for larger beads, and the potential approaches hard-sphere-like behaviour.Following this approach, at an upper limit -in practice at N CG,limit o 10 87 -the solvent freezes.It should be noted that DPD mappings with high N CG disregard the compressibility matching while matching other parameters, such as density and pressure.This allows DPD approaches to model also significantly higher N CG values than strictly allowed by the compressibility consideration.Since here we have N CG = 8, the choice a ii = 25 k B T is well motivated.The cross-repulsion parameters a ij were obtained via mapping the interactions to the Flory-Huggins mixing parameter w ij as Here z = 0.101 AE 0.001 for DPD bead density r Z 3. In the current work, the Flory-Huggins mixing parameter is calculated from Hildebrand solubility parameters 88 as In this, V bead is the DPD bead volume, while d i and d j are the Hildebrand solubility parameters for components i and j, respectively.The Hildebrand solubility parameter d is defined based on cohesive energy E coh , that is, the energy difference between the condensed phase and the gas phase of a component.
Here, E coh was determined based on all-atom MD simulations.In these, the CHARMM-C36 force field 89,90 and the compatible Tip3p all-atom water model 91,92 were used.4][95] simulation package.For the atomistic detail simulations values, separate condensed phase simulations consisting of 180 octane, 172 2-methylheptane, 122 sulfosuccinic acid, and 1265 water molecules, were first energy minimized for 1000 steps using the steepest-descent method.
The molecular species were chosen such that parameters for DPD parametrization corresponding to apolar hydrocarbon solvent (octane), polar solvent (water), and the AOT molecule composed of a head bead and two tail beads can be obtained.Notably, as all bead types in standard DPD are identical in mass and volume, the resulting beads correspond to the massequivalent of 8 water molecules of the parametrized species.In the model, the headgroup of the AOT molecule is considered as bound with the Na + counterion.
Equilibration was achieved via 5 ns NPT ensemble simulations with timestep Dt = 2 fs.In these, temperature was controlled using the stochastic velocity rescale thermostat developed by Bussi et al., 96 with reference temperature T = 300 K and t T = 0.5 ps.A target pressure of p = 1.0 bar was maintained by isotropic Berendsen barostat 97 with t P = 2.0 ps and compressibility of 4.5 Â 10 À5 bar À1 .
Next, the non-bonded energy terms (Coulombic and Lennard-Jones contributions) were sampled during a 2 ns NVT ensemble simulation at T = 300 K. Additionally, single molecule gas phase, i.e. vacuum (in vacuo) simulations of individual sulfosuccinic acid, 2-methylheptane, and octane molecules, as well as, the 8 water molecules considered for single DPD bead, were set up using the exact simulation box dimensions of the condensed phase NVT run.This choice was made so that the system volume dependent part of the particle mesh Ewald method, 98 used for electrostatic interactions, would be for same sized box.Non-bonded energy terms were sampled during a 2 ns NVT ensemble run at T = 300 K. E coh was then calculated as the difference in non-bonded energies of the condensed phase and the in vacuo simulations.Notably, reliable values of E coh require a well equilibrated condensed phase simulation.The solubility parameter can then be calculated as , where V m is the molar volume.For tabulated nonbonded energies and solubility parameters, see Table S1 in the ESI.† The resulting DPD conservative repulsion parameters a ij are presented in Table 1.
Intramolecular T-H bonds in the three bead CG representation of the AOT molecule, Fig. 1, are described by a simple harmonic potential À1 and the inter-bead bond equilibrium length r 0 = r c .

System set-up and analysis
Simulations were performed in a cubic box of (60 r c ) 3 with r c E 0.9 nm.Periodic boundary conditions were applied in all dimensions.With r = 3, this corresponds to 648000 beads per simulation box.For the initial configuration of each system, 14547 AOT (equivalent to a molar concentration of 150 mM or 6.7 w/w%) were randomly distributed within the simulation box.This particular concentration is selected for this study because it is an extremely common concentration used in various different types of experimental setups involving AOT, both in fundamental research and technological processes, see e.g.ref. 99-102.Additionally, the selected concentration exceeds critical micellization concentration (CMC) throughout the examined solvent polarity range. 47,103,104The remainder of the box was subsequently filled with monomeric solvent.Simulations were performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) MD program. 105Molecular visualization were produced using VMD. 106Each simulation was run for 5 Â 10 6 steps with a time step Dt = 0.01.For data analysis, the last 5 Â 10 5 simulation steps were used.
During the beginning part of the trajectory, the structural and assembly characteristics equilibrated, as measured by cluster size distributions and order parameter time evolution analysis.Solvent was varied in the simulations by varying a TS and a HS , such that their values were scaled in steps, to correspond to a solvent environment from an apolar hydrocarbon solvent (system a 1 ) to polar water (a 12 ).The corresponding AOT miscibility parameter pairs a ij are summarized in Table 2. Sensitivity of the assembly response to balance between a TS and a HS in the scaling was examined for systems a 1 , a 3 , a 6 , and a 9 , by variation of the solvent-surfactant interactions.Table 3 summarizes the range of sensitivity analysis and provides the sampled parameter values.All other interactions were kept constant.The detailed results of the sensitivity simulations can be found in the ESI.† The number of aggregates and aggregate sizes were determined using the DBSCAN clustering algorithm, 108 which works on the assumption that clusters are dense regions in space separated by regions of lower density.The sklearn DBSCAN clustering function 109 was used with parameters eps = 2 r c and min_samples = 25.eps defines the maximum distance between two beads to be considered a part of a cluster, accounting for periodic boundary conditions.min_samples defines the minimum number of points to form a core point.The internal structure of the assemblies was characterized using the order parameter where x is the number of beads of type j around beads of type i (i a j) within a defined cutoff radius (r cut = r c ). x 0 is calculated for the initial state (first simulation frame at t = 0), where beads i and j are considered to be in the most mixed state, while x is calculated based on the equilibrated system.c i varies between 0 and 1, with c i = 0 corresponding to a random mixture and c i = 1 to a fully phase separated system within the length scale defined by r cut .

Results and discussion
Fig. 2 provides a visual summary of the AOT surfactant selfassembly response observed in different solvents.In the modelling, the solvent is varied by the differing miscibility of the AOT surfactant headgroup and tails with the solvent, as captured by the a HS and a TS parameters.Data shows significant assembly response changes due to solvent polarity ranging from the a 1 (hydrocarbon solvent) system, that forms oligomeric aggregates and has free AOT dissolved in the solvent phase, to the aqueous solvent a 12 , which results in formation of well-defined lamellae.
For the a 1 (hydrocarbon solvent) system, very small AOT aggregates with average size N agg = 7.8 form, see Fig. 3 for the aggregate size distribution.The observed aggregate size is smaller than the N agg E 10-30 previously reported for AOT reverse micelles in near water-less hydrocarbon solutions, 15,24,33,[110][111][112] including benzene (N agg = 9-12), 111,112 iso-octane (N agg = 16-21), 112,113 cylohexane (N agg = 14-22) 15,112 and dodecane (N agg = 29). 110However, the notable variance in reported aggregation number for reportedly ''water-less'' systems, reflects the challenge of accurately quantifying the sizes of small non-hydrated reverse aggregates experimentally.Also, the presence of even trace amounts of water in the system greatly affect the size, morphology, and dynamics of the formed reverse micelles. 112,114The exponential distribution of aggregate sizes, see Fig. 3, resulting from open association aggregation process, makes the average aggregate size a poor descriptor of the system. 115,116Additionally, a significant portion of the AOT remains monomerically dissolved in the solvent phase, see Fig. 4.
Further analysis of the formed aggregates by the radial distribution function for AOT head-beads g HH (r) in Fig. 7 indicates that the formed aggregates are spheroidal with an average diameter of the AOT headgroup core of approximately   2.4 r c E 2.2 nm.This is slightly larger than the 1.7-1.9nm previously reported based on scattering experiments 14,15 and scaling laws 15,117 for low hydration AOT reverse micelles.The discrepancy is likely explained by the difference in the actual molar volumes corresponding to the molecular regions mapped to the head and tail beads, respectively.As mentioned, beads in DPD are assumed to be equal in volume and in mass.Another source of disagreement is the soft nature of the DPD potential, which allows for bead overlap.Further comparison of the AOT apolar solvent self-assembly response is complicated by most existing recent characterizations of AOT reverse micellization in solvents comparable to the apolar hydrocarbon solvent, here mainly concerning systems where water is present, at 3 to 20 water molecules per AOT concentration. 13,58,63These studies indicate that AOT reverse micelles undergo morphological changes in shape and size from elongated, highly dynamic filament-like aggregates at low water content (B3 water molecules per AOT), to stabilized spherical reverse micelles at high water content (47 water molecules per AOT).In the current work, we considered only dry hydrocarbon solvents as DPD level description of water leads to significant coarse-graining which can be expected to decrease modelling accuracy.Typically, an atomistic detail or hybrid modelling approach is needed, see e.g.ref. 118.
Although not covered by this work, as trace water is commonly present in AOT-apolar solvent systems, it is worthwile to consider literature further: In spherical reverse micelles with a sizeable water core, most of the water behaves as bulk water, and remains distanced from the AOT headgroup-tail interface.At lower water content, most of the water is interfacial, and aggregate growth occurs primarily via water bridging along one major axis, resulting in elongated aggregates. 13,27,38,119Previous modelling studies also reveal that the degree of micelle asphericity and elongation at low water content is larger for aggregates allowed to self-assemble from isotropic mixture than those relaxing from pre-built initial aggregates. 13,58,63he findings can be expected to generalize to AOT aggregation in the absence of water, especially if other polar impurities or additives are present.For example, several non-aqueous polar solvents, such as short-chain alcohols, form reverse micelles [120][121][122][123] and act as cosurfactants. 123,124For example, AOT/non-polar solvent/alcohol tertiary systems have revealed that alcohol molecules, such as methanol, not only bind to AOT headgroups, but can also solvate AOT hydrocarbon tails and disperse in the apolar solvent phase. 122,123Similar findings for alcohol acting both as cosurfactant and cosolvent exist for other surfactant species. 125ncreasing polarity of the solvent results first into reduced miscibility of the AOT surfactant.This shows as the sharp reduction of free monomers in Fig. 4, but also the formation of large, well-defined individual aggregates, as can be seen from the simulations configurations visualizations in Fig. 2. The phase separation resulting from reduced miscibility of the AOT surfactant is also apparent in the sudden increase in the tail bead order parameter value c T , see Fig. 5. Specifically, the higher head and tail order parameter values compared to the reverse micelle phase suggest compartmentalization and ordering within the larger aggregate.This internal ordering is evident from the aggregate cross-sections presented in Fig. 6, where the formation of AOT headgroup rich regions within a continuous phase of AOT tails is evident.The change in head bead order parameter c H is less pronounced, due to similar clustering of the AOT head beads being present already in the reverse micellar phase corresponding to a 1 .Both tail and head order parameters plateau at a maximum value for solvents corresponding to a 4 and higher polarity.Furthermore, the assembly response taking place in intermediate polarity solvents allows a sparser AOT headgroup packing than in the tightly packed reverse micelle aggregates of the apolar solvent.This shows as the initial peak of head-head radial distribution function g HH (r) shifting to higher values of r in Fig. 7. Formation of the large aggregates also introduces a second, broad peak in the g HH (r) data at r B 2.8 r c .This is related to the inter-spacing of the AOT headgroup rich clusters within the aggregate.
As solvent polarity is increased further, the hydrophilic AOT headgroups migrate to the surface of the aggregate.This is Fig. 4 Fraction of free AOT surfactant monomers in solvents of different polarity.The solvent polarity variation is modelled by a i , see Table 2 for corresponding a HS and a TS parameter combinations.
Fig. 5 Order parameter c calculated for head (c H ) and tail (c T ) beads of the AOT assemblies formed in solvents of differing polarity indicated by a i , see Table 2 for corresponding a HS and a TS parameter combinations.
evident from the shift of the first peak in the head-solvent radial distribution function g HS (r) to smaller values of r, and a notable increase in the peak height, see Fig. 8.This behaviour is reversed in the tail-solvent radial distribution function g TS , with the AOT head beads shielding interaction of the hydrophobic tails with the hydrophilic solvent, see Fig. 9. Furthermore, the order parameter dependencies on the solvent polarity in Fig. 5 reveal that partitioning of the AOT headgroups between intra-aggregate clusters and the aggregate surface causes a sharp decrease of the head order parameter c H in the solvent polarity range corresponding to that between a 10 and a 11 systems.This corresponds to the assembly response switching to a layered, onion-like surface structure for the elevated polarity solvents, see the cross section visualizations for parameter set a 11 in Fig. 6.Finally, in highly polar solvents, such as system a 12 , AOT forms a lamellar phase, resulting in further decrease of head order parameter c H value, Fig. 5.As the AOT tailgroups form the inner, well-organized region of the lamella, the tail order parameter c T retains its high value.Additionally, at high solvent polarities, a slight increase in the fraction of free AOT surfactant can be observed.However, the amount of free AOT could in the simulations be affected by the compressibility characteristics of the simulation and the periodicity of the simulation box.
Phase diagrams for AOT/water systems are well established in literature: 55,56,126,127 AOT/water binary systems at ambient temperatures undergo phase changes from isotropic mixture, to lamellar phase, to bi-continuous cubic phase, and finally, to inverse hexagonal phase with increasing AOT concentration.In our simulations, at 150 mM/6.7 w/w% AOT, the system should be in the meta-stable coexistence region of the isotropic Fig. 6 Cross-section cuts of representative self-assembled AOT aggregates at different solvent conditions a i , see Table 2 for corresponding a HS and a TS parameter combinations.The AOT head beads are shown in cyan and tail beads in purple and all cuts are made along the (x,y)-plane through the center of mass of the major visualized aggregate.Solvent beads are omitted in the visualizations for clarity.Fig. 7 Radial distribution function g HH (r) calculated between the AOT head beads.The variable a i refers to solvent polarity variation, see Table 2 for corresponding a HS and a TS parameter combinations.
Fig. 8 Radial distribution function g HS (r) calculated between the AOT head beads and the solvent beads.The examined systems a i differ in solvent polarity, see Table 2 for corresponding a HS and a TS parameter combinations.
Fig. 9 Radial distribution function g TS (r) calculated between the AOT tail beads and the solvent beads.The examined systems a i differ in solvent polarity, see Table 2 for corresponding a HS and a TS parameter combinations.
and lamellar phases. 55,126Indeed, we observe a rippled, yet stable, lamellar structure with some free AOT monomers in the hydrophilic, water-like solvent phase.Similar rippled AOT bilayers in the metastable region of the phase diagram have been previously reported, both experimentally 45,[128][129][130] and computationally. 50Experimentally, values of 1.8-1.96nm for AOT bilayer thicknesses have been measured in water. 45,49,128isually, the simulations assemblies correspond to approximately two DPD beads, which is consistent with the experimental values.Deviations can, however, be expected due to the lack of angle potentials within the DPD AOT model employed here combined with the ability of the DPD beads to overlap.Especially at high packing density, such as in the lamellar phase, this may result in nonphysical molecular conformations when mapped to an atomistic representation.We note that the employed DPD approach omits, e.g., electric double layer effects and solvation layer changes in the presence of charges.Spacing between the AOT headgroups in each lamella leaflet is B0.75 r c E 0.68 nm, which is in good agreement with values of 0.64-0.67nm from previously published experimental measurements. 45,128This means that even in the absence of proper electrostatic interactions description in the DPD bead model, particularly the model lacking accuracy in describing the repulsive interactions between the negatively charged sulfonate headgroups, the model presented here correctly reproduces the packing characteristics of the AOT bilayer.
To summarise the assembly characterization in terms of varying solvent, the here presented DPD parametrization for AOT surfactant in solvents of differing polarity reproduces in good agreement with literature structural observables for the isotropic reverse micelle phase (a 1 ) and the lamellar phase in water (a 12 ), corresponding to the two extrema of the examined solvent polarities.The agreement covers aggregation numbers and AOT bilayer topology, bilayer thickness, and headgroup spacing.This means that even a simplified DPD parametrization, with large enough degree of coarse-graining to achieve assembly phase level modelling, can capture response typical to model surfactant such as AOT.
For the solvent phases of intermediate polarity (a 2 -a 11 ) literature sources on matching AOT-solvent binary systems are sparse and therefore accuracy of the model is challenging to judge.In these systems, we observe AOT agglomeration into singular large aggregates with internal compartmentalization of AOT headgroups within a continuous phase formed by the AOT hydrophobic tails.The response is a direct outcome of the competition between the miscibility and immiscibility characteristics in a dual-natured molecule, such as the surfactant.The response is well known for e.g., block-copolymers.With increasing polarity of the surrounding solvent medium, AOT headgroups form a hydrophilic outer layer for the aggregates.Notably, solubility of AOT monomers in bulk solvent decreases sharply with increased solvent polarity.
To further quantify the differences rising from the interactions with different solvents, we also performed simulations consisting of only two CG AOTs and the solvent beads in the DPD simulation box.In this low-density (LD) limit, it is possible to calculate the effective interaction potentials between the centres of mass (CM) of the AOTs, f eff CM , via the LD limit CM radial distribution function g LD CM as bf eff CM (r) = Àln(g LD CM (r)), (10)   where b = (k B T) À1 .The essential simplification here consists in assuming the interactions between AOTs to be isotropic when extracting the effective interaction potential.The results are obtained by averaging g LD CM (r) over 20 independent runs, each one lasting 8 Â 10 6 steps with an integration time step Dt = 0.05 in a cubic box of (60 r c ) 3 .Fig. 10 shows the resulting effective interaction for the cases a 1s2 , a rep , and a 11s1 , see Table 3, as well as the corresponding equilibrium configurations from finite density DPD simulations.Interestingly, all extracted effective interactions measured for the centres of mass of the AOTs exhibit a Lennard-Jones-like form, i.e. composed by a shortranged repulsion and medium-ranged attraction that decays relatively fast with increasing separation.The phase response of such systems is known to exhibit, e.g., gas-liquid coexistence. 131Additionally, the a 11s1 system also exhibits a second, now long-ranged, repulsive barrier.Such form, often referred to as short-ranged attraction and long-ranged repulsion (SARL), is known to give rise to microphase-separation, 132 mixed particles in binary mixtures, 133 and clustering in both passive 134 and active systems. 135Note that in the following discussion we assume that the common phase response of Fig. 10 Low-density limit effective potentials between pairs of AOTs, bf eff CM , for the cases a 1s2 , a rep , and a 11s1 , and the corresponding finite density equilibrium configurations.
Lennard-Jones-like particles is valid also for the surfactant systems.
Comparing the effective potentials and the DPD simulated 150 mM AOT concentrations in Fig. 10 reveals that the strongest attraction, a rep , leads, as expected, strong aggregation and well-defined assembly structure.The intermediate attraction well, combined with the second long ranged repulsion barrier, a 11s1 , rises from the system assembling to deformed, flexible aggregates with significant portion of the AOT remaining soluble as monomers and oligomers.Most interestingly, the weakest effective attraction, with no barrier, corresponds to the small clustered aggregates and high solubility of the AOT.The comparison shows that such effective potentials can by their form and attraction strength indeed be connected with the equilibrium assembly structures.For example, the a 1s2 system results in a significant attraction well region with depth EÀ1.8k B T in the effective potential.Such form clearly indicates that at low enough temperature, the effective system will phase separate according to basic thermodynamic considerations.However, the corresponding equilibrium configuration exhibits small clusters instead.Furthermore, the probability distribution of the cluster size decreases quite monotonically as the function of aggregate size, data presented in the ESI † as Fig. S3.This signifies that the binary AOT-solvent system is in a thermodynamic state outside the binodal line (assuming that such line exists in the surfactant system) in the miscibility region, but forming correlated clusters due to the attraction.On the other hand, the case a rep shows the highest attraction strength with a minimum of the potential of EÀ4.2k B T. At the concentration used to run the bulk DPD simulations, the AOT system exhibits a strong phase separation, in which the AOTs aggregate to a single large cluster in the bulk DPD simulations.This means that the thermodynamic state point chosen here lies below the miscibility line and inside the phase coexistence region of the phase diagram.Within the cluster, due to the symmetric interaction between head and solvent, and tail and solvent, heads and tails are well mixed (see Fig. 10).Finally, in a 11s1 , the relatively strong attraction combined with the second repulsion range results in a rather interesting assembly response: here the attraction strength is enough to enforce phase separation, as shown in Fig. 10.Additionally, compared to a 1s2 , a much higher number of free AOTs exist in solution.This is expected, since for decreased attraction strength, at fixed temperature, the concentration corresponding to the coexisting dilute phase moves to higher values (respectively, the dense phase concentration is decreased), see e.g.ref. 131.Interestingly, the cluster also exhibits an internal layered structure.These features correspond to a second length-scale emerging in the system.For the assembly, the second length scale characterises the spacing between the lamellae.The response can be connected with the repulsive tail in the effective interaction potential, on the order of E2 r c , but also depends on other system conditions.
The three example cases, for which we considered the assembly response via effective isotropic interactions based model, show that the functional form of f eff CM clearly gives relevant insight on the expected assembly phase for surfactant systems.Hypothetically, knowing precisely the effective pair potential allows the calculation of an effective phase diagram.This can be done in different ways, including simulations of spherical particles interacting via f eff CM (r) or via energy minimization of the corresponding free energy.The level of phase response predictability of this CG approach can be expected to correlate with the degree of symmetry of the real interactions between the surfactants.
Further insight to the properties of f eff CM can be obtained by considering the corresponding second virial coefficient B in the quadratic virial expansion, reading where p is pressure in the particle system and r is the number density of the CG surfactants.For our surfactant systems, rAOT = 0.067 r c À3 .B can be calculated directly via Evaluating this gives, for all three examined sample cases, a negative value for B. This indicates predominantly attractive interaction between pairs of CG particles.Here, the largest negative value is B a rep E À134 r c 3 , followed by an intermediate value B a 11s1 E À33 r c 3 , and finally B a 1s2 E À14 r c 3 .Note that the numerical integration has been truncated at r = 4 r c .The truncation is done since the small deviation of g LD CM (r c r c ) from unity becomes a source of numerical error when weighted by the factor r 2 in eqn (12).The trend of the values of B matches with the propensity of the surfactant systems to self-assemble (as discussed above).To summarise, the significance of this is that even modelling considerations more simplistic than the mesoscale DPD approach here can extract useful information on the assembly.Additionally, the effective potentials demonstrate the origins of dual length scale response observed here.
Next we present a critical overview of the DPD model and possible improvements to it.In this work, the DPD conservative repulsion parameters were derived based on determining Hildebrand solubility parameters from atomistic detail MD simulations data and employing the Flory-Huggins solubility theory to back-map the solubility parameters to DPD interactions.The Hildebrand solution theory is viable for non-polar molecules, but fails to accurately describe molecules with high polarity, hydrogen bonding capability, or charged moieties, such as the AOT sulfonate headgroups.For such cases, use of alternative solubility parameter descriptions, such as the Hansen solubility parameter, have been proposed. 72,136owever, these do not have a similar physical basis as the Hildebrand solubility theory, and are therefore limited by the accuracy of group contribution methods mapping or availability of suitable experimental data for empirical estimation of the cohesion energy contributions.Anderson et al. 137 have proposed tuning of DPD interaction parameters to reproduce experimental water-octanol partition coefficients.While this approach accurately reproduces, e.g., CMC for polyoxyethylene derived surfactants in water, average aggregate size is not accurately reproduced.This is mainly due to the dependence of the DPD parametrization on molecular bond length.Transferability of the approach by Anderson et al. to surfactants in apolar solvent without distinct CMC behaviour remains currently unexplored.Hence, for this work, we chose the Hildebrand approach and obtained with it assembly response that matches to a good degree to available literature data and expected response.Further advances could be made by expanding the work to other solubility estimates.
We focused the examination of the model performance on the literature reported structural assembly characteristics.The current parametrization could clearly be fine tuned to match experimental observables such as CMC, partition coefficients, 137 or structural data, such as radial distribution function of binary systems. 84However, additional tuning may deteriorate the response in some other sense, and we decided against this in a general model.For example, increasing DPD bond stiffness often leads to better match to MD derived structural data at the cost of significantly reducing time step length and therefore limiting available simulation time. 84,1380][141][142] While incorporation of headgroup charge through, e.g., charge smearing approaches would be a viable option, 143,144 this implies a uniform dielectric screening environment, which is not the case in a highly structured system, especially in internally structured aggregates such as those described in this work.Smearing of the headgroup charge also affects equilibrium selfassembly structures. 78Finally, most surfactants with strongly polar or charged headgroups are hygroscopic, which implies trace amounts of water present even in ''dry'' binary systems. 145xplicit description of the charged headgroup and its counterion would require addressing the degree of hydration of the counterion and surfactant headgroup.
Furthermore, the DPD model described here utilizes bonded potentials only for the AOT head-tail bonds, similar to previously published DPD surfactant models. 75,146Increased conformational accuracy could clearly be achieved by matching the parameters for the bonded potentials, especially the spring constant and equilibrium bond length, to, e.g., atomistic simulations data.Most commonly, radial distribution function of the centers of mass of a group of atoms is targeted. 84Furthermore, introduction of additional angle and dihedral terms may be important for accurate modelling of, e.g., surfactant systems, for which properties such as CMC and partitioning have been shown to depend on chain rigidity and equilibrium bond length. 147,148Here, the lack of angle terms in the AOT structure could contribute to the difference in thickness of the AOT bilayer in comparison to priorly reported values.The CG description also results in failure to capture hydrogen bonding related contributions to solvation and aggregate formation, demonstrated in our previous works to influence reverse micellar aggregation in apolar solvents with hydrogen bonding ability. 115,116 the model parametrization, the Na + counterion is considered effectively bound to the headgroup of AOT surfactant, resulting in a lack of charged interactions and an explicit description of ionic species in the DPD bead model, but also allowing larger degree of coarse-graining which gives access to extended length and time scales, Fig. 2.This means that electrostatics interaction-mediated aggregate growth, as well as AOT sulfonate headgroup repulsion, are captured only marginally by the effective parameters in this level of CG model.The anionic sulfonate headgroup dominates the water/AOT interface interactions in both the AOT reverse micelles and lamellar structures, with a full solvation shell of B6 water molecules per sulfonate headgroup. 48,149Especially for dry reverse micelles, such as those corresponding to system a 1 of this work, AOT sulfonate headgroup repulsion plays an important role in stabilising aggregates of finite size.Additionally, reverse micelles may gain net charge from micelle collision events. 150Such phenomena are not captured by this model.While attempts to include electrostatic interactions between charged beads via distributed charges exist, 78,143,151,152 the parametrization of these interactions for specific DPD systems remains challenging.

Conclusions
We have presented, and verified against literature, a bottom-up DPD parametrization for sodium sulfosuccinate (AOT) surfactant in solvents of varying polarity, from apolar hydrocarbon solvent (octane) to polar solvent (water).The parametrization of the model is based on mapping Hildebrand solubility parameters, obtained based on atomistic MD simulations, and mapped to DPD conservative force repulsion parameters via Flory-Huggins solution theory.Additionally, we assessed the sensitivity of the derived DPD parametrization and compared the ability of simple, effective isotropic modelling approaches to capturing the AOT assembly response, reporting that, although the finite density phase response bears clear dependency on the solid concentration, such effective interaction potential derived from the low density limit response can be used as a simplified tool to predict the possible assembly structures.
Despite considering a relatively large degree of coarsegraining, N CG = 8, the DPD model presented here reproduces well structural phases of AOT in apolar and polar solvents, where comparable literature data is available.Additionally, the model allows large-scale assembly modelling and characterization of the changes, not accessible at lower N CG .While additional tuning of the model is likely to be required for investigation of, e.g., viscoelastic response and dynamical properties, more complex ternary systems, or systems perturbed via an external field (e.g.electric field), the current model describes binary AOT-solvent systems on time and length scales currently unattainable via atomistic models.The predicted AOT selfassembly response can be used to design surfactant-solvent systems of desired assembly response.Furthermore, the diverse

Fig. 1
Fig. 1 CG DPD bead representations of the AOT surfactant and the octane solvent used as the model hydrocarbon solvent in this study.

Fig. 2
Fig. 2 Simulation snapshots corresponding to AOT self-assembly (150 mM concentration) in the 12 different solvent compositions a i studied here.The corresponding a HS and a TS parameter combinations are summarized in Table2.The AOT head beads are shown in cyan and tail beads in purple.For all systems, the visualization corresponds to the final simulation frame.The solvent beads are omitted in the visualization for clarity, and the boundaries of the periodic simulation box are marked as blue lines.For all systems except a 1 , where only very small aggregates form, periodic images are used to visualize the aggregates.

Fig. 3
Fig. 3 Probability distribution for observed reverse micellar aggregates in a 1 solvent system.

Table 1
DPD conservative force repulsion parameters a ij for the different bead types.a TS and a HS represent the tail-solvent and head-solvent repulsive parameters that are varied in the work

Table 2
Summary of the repulsive parameter a ij values examined to capture the effect of solvent variation.The values for H-S and T-S DPD conservative force repulsion parameter have been scaled based on the apolar (hydrocarbon) solvent and polar (water) solvent parameters to model the interactions in solvent environments of intermediate polarity.In the parameter range, the system a 1 matches a hydrocarbon solvent (octane), and a 12 polar solvent (water).Back-mapping of the intermediate solvent range to some chemically specific solvents based on Hildebrand solubility parameters 107 results in the solvation conditions of system a 3 correspond to approximatively diethyl ether, a 5 to toluene, a 6 to benzene, and a 10 to methanolaij a 1 a 2 a 3 a 4 a 5 a 6 a 7 a 8 a 9 a 10 a 11 a 12 Head-solvent 71 67 66 61 56 51 46 41 36 31 25 21 Tail-solvent 25 30 31 38 44 51 58 63 70 76 83 89 This journal is © the Owner Societies 2023 Phys.Chem.Chem.Phys.

Table 3
Summary of the repulsive parameter a ij pairs for the headsolvent and tail-solvent interactions examined for sensitivity analysis and comparison with the purely scaling obtained values of Table2a ij a 1s1 a 1s2 a 3s1 a 3s2 a 6s1 a 6s2 a 9s1 a 9s2 a 11s1 a rep Head-solvent 60 80 56 76 41 61 29 33 25 80 Tail-solvent 25 25 31 31 51 51 76 76 80 80 This journal is © Owner 2023 Phys.Chem.Chem.Phys.assemblymorphologies in the solvents of intermediate polarities indicate that advanced functionalities from AOT assembly could be achieved by solvent variation.