Quantitative prediction of oil–water interfacial tension in surfactant systems using dissipative particle dynamics

Abstract

We present a dissipative particle dynamics (DPD) model for surfactants at oil–water interfaces, parametrised directly against experimental data. The model is applied to pure and mixed systems of ionic and nonionic surfactants, including sodium dodecyl sulfate (SDS), dodecyldimethylamine oxide (DDAO), and polyoxyethylene alkyl ethers (C_iE_j). Our simulations reveal that exceeding the maximum interfacial packing density can lead to the spontaneous production of micelles from the interfacial surfactant layer into the bulk. This behaviour, together with our parametrisation scheme, enables quantitative prediction of interfacial tension as a function of surface concentration in excellent agreement with experiment. In addition, the model provides access to other interfacial properties such as maximum surface coverage and interfacial monolayer thickness. Moreover, we show that the model can be used to investigate synergistic effects in mixed surfactant systems, where combinations of surfactants yield lower interfacial tensions than either component alone. In particular, simulations of DDAO and SDS systems demonstrate a lower interfacial tension than the corresponding pure surfactants. Utilising these results, we propose that these synergistic effects result from a balance between the intrinsic ability of a surfactant molecule to lower the IFT and its achievable maximum packing at the interface.

---

1 Introduction

Surfactants are amphiphilic molecules that serve as primary components of cleaning detergents, as well as playing key roles in industries such as cosmetics,¹ food,^2,3 and medicine.^4,5 Their widespread industrial uses mean that surfactants have been the focus of numerous computational^6–17 and experimental^18–29 studies across the literature. This work includes phase diagram mapping;^21,27,29 characterising micellar properties^11–14 (e.g., size, shape and aggregation number); calculating the critical micelle concentration (CMC);^11,12 and rheological studies of surfactant solutions under different conditions.^22,23,28

The amphiphilic structure of surfactants drives their adsorption at interfaces, such as oil/water and air/water boundaries. This adsorption significantly reduces the surface (or interfacial) tension. This property is central for their use in cleaning products, where lower surface tension generally leads to increased wettability, and therefore more effective cleaning. A large body of experimental research has been dedicated to determining the interfacial tension (IFT)^18–20 of a wide variety of surfactants as a function of their concentration. Other key interfacial properties include the maximum packing density of surfactants at an interface,^18–20,24 where they achieve their maximum reduction in interfacial tension.

Computational methods, such as molecular dynamics (MD),^{9,15,17,30–33} can be used to study the behaviour of surfactants at interfaces. Such approaches can help provide insight into behaviours that are difficult to determine via experimental means alone. For example, simulations have been proven to be useful for determining conformational and orientational behaviour of molecules at an interface,^15,30,33 which is particularly difficult to infer from experiments. Computational approaches can also be used to make quantitative predictions about IFT as a function of surfactant surface concentration,^{7,9,15,17,30–32,34} or as a function of molecular properties such as tail length.^7,35

However, highly detailed MD studies are typically limited to small time and length scales and it is common to pre-assemble molecules at the interface at the start of the simulation, due to the prohibitively long time scales in letting molecules adsorb to the interface naturally. As a result, it can be difficult to determine properties such as the maximum packing density of surfactants at the interface with high accuracy.^9,36 To overcome these time- and length-scale limitations, coarse-grained approaches are frequently employed to study interfacial behaviour, taking advantage of the speed of approaches such as dissipative particle dynamics (DPD)^{6–8,37–44} or coarse-grained MD (e.g. models such as MARTINI^44–49).

In this work, we develop a detailed new DPD parametrisation scheme to study the interfacial behaviour of common surfactants at oil/water interfaces. We model nonionic, anionic, and zwitterionic surfactants, and address the known problem that surfactants containing ethylene oxide (EO) units are difficult to parametrise for DPD simulations^50,51 (due to polyethylene oxide being extremely soluble in water). We use an approach in which EO interactions are parametrised based on a combination of experimental data, including partitioning behaviour^52,53 and infinite dilution activity coefficients.^54,55 To model interfacial tension well, it is also essential to reproduce bulk densities and to account for the effects of bonding on DPD bead parametrisation. Our approach takes all these factors into account and yields excellent agreement with experimental IFT values for pure surfactants at water–oil interfaces. We are also able to study mixed surfactant systems, allowing us to predict synergistic effects, which are experimentally observed for such systems.^19,20,25

2 Dissipative particle dynamics (DPD)

2.1 Overview

Dissipative particle dynamics is a coarse-grained modelling method, where atoms are grouped into ‘beads’, and a single bead represents several atoms. By bonding beads together, one can represent large and complex molecules. DPD equilibrates far more rapidly than traditional MD because of the reduction in the number of sites and the design of pairwise interaction forces.

DPD beads primarily interact via a soft repulsion, where the strength of this repulsion is varied to reproduce the correct chemical behaviour of different bead types. This repulsion is referred to as the conservative force and acts up to a cut-off R_C. The conservative force takes the form

where a_ij is the magnitude of the repulsion, the vector between beads i and j is calculated as r_ij = r_i − r_j and the unit vector is defined as r̂_ij = r_ij/|r_ij|. Beyond the cut off R_C, the conservative force F^C_ij = 0. We note that this soft interaction potential allows for larger time steps than traditional MD.

Other forces are included to ensure that beads reproduce the behaviour of a realistic fluid. The total force acting on bead i is given by

where the sum is over all other beads j in the system. F^R_ij is the random force, F^D_ij is the dissipative force, and F^E_ij is the electrostatic force acting on charged beads. Two additional forces are included: F^B_ij for bond stretching between connected beads and F^A_ij for angle bending on bead triplets, which together enforce molecular connectivity and stiffness.

The random force introduces thermal fluctuations and is given by

F^R_ij = σω^R(r_ij)ζ_ijr̂_ijΔt^−1/2, (3)

while the drag force replicates viscous effects in the fluid, and is calculated as

F^D_ij = −γω^D(r_ij)(r̂_ij·v_ij)r̂_ij. (4)

The functions ω^D and ω^R are weight functions that vanish for distances above the cut-off distance |r_ij| > R_C. γ is a friction coefficient and σ is the noise amplitude, v_ij = v_i − v_j is the velocity between beads i and j, Δt is the time step, and ζ_ij(t) is a randomly fluctuating Gaussian variable with zero mean and unit variance. The weight function used in LAMMPS takes the form

To satisfy the fluctuation–dissipation theorem, the relationship between the weight functions must obey⁵⁶

ω^D = [ω^R]² (6)

and the relationship between the amplitudes is

σ² = 2γk_BT, (7)

where k_B is the Boltzmann constant and T is the temperature. In this work, we set values for the constants σ = 3 and γ = 4.5 and choose a time step of Δt = 0.01.

To model the electrostatic pair potential, we use Slater-type charge smearing, where the potential U_E between two charged beads i and j is given by

where q_i and q_j are the charges, Γ = e²/(k_BTε₀ε_rR_C) is a dimensionless electrostatic coupling parameter, and β* is the tuneable Slater parameter. We set β* = 0.929R_C⁻¹ in line with the work conducted by González-Melchor⁵⁷ and use ε_r = 78.2 for water at room temperature. The long-range interactions are computed by the particle–particle particle–mesh (PPPM) technique, where we set the real-space Ewald cutoff as 3.0R_C.

For the bonding forces F^B_ij and F^A_ij, we use harmonic spring potentials of the following form:

U^B_ij = K_B(r_ij − l₀)², (9)

U^A_ijk = K_A(θ_ijk − θ₀)², (10)

where l₀ is the equilibrium bond length, θ_ijk is the bond angle between beads i, j and k, and θ₀ is an equilibrium angle. The constants K_B and K_A define the strength of potentials, which will be defined in later sections of this article. We note that the usual ½ factor is included in K_B and K_A.

DPD convention means that simulations are conducted in reduced (dimensionless) units, where energy is measured in units of k_BT, length is in units of the cut-off radius R_C, and mass is measured in units of the mass of a single DPD bead. Therefore, in our simulations, we choose standard choices for the reduced units setting mass m = 1, the energy as k_BT = 1, and the cut-off distance is set as R_C = 1.

We perform all of our simulations using the LAMMPS software (Version 28 Mar 2023 – Update 1). In LAMMPS, bonded beads (by default) do not interact with each other through nonbonded interactions. Although it is possible to alter this through the use of the ‘special bonds’ setting, we choose to leave the interactions as their default behaviour. This is an important consideration when comparing with other DPD studies, as this choice has an impact on the IFT values calculated.

2.2 Defining DPD parameters

In the DPD framework, one can represent different bead types by altering the magnitude of repulsion a_ij in the conservative force given in eqn (1). Therefore, an a_ij value needs to be determined for each bead pairing i and j. One of the ways of defining a_ij is to tune its value to reproduce the correct infinite dilution activity coefficient (IDAC). This ‘target’ IDAC value can be obtained from experimental measurements or additional computational methods, such as using COSMO-RS software. Alternative approaches include using Flory–Huggins parameters, which are calculated via experimental solubility parameters.^8,37

In DPD, the IDAC (γ^∞) of a solute of type A at infinite dilution in a solvent of type B, can be calculated using the following relation (where the full derivation is provided in our previous work⁷),

where μ^A_ex is the excess chemical potential of the solution, μ^⊖_ex is the excess chemical potential of pure substance A, and ρ_A and ρ_B are the number densities of A and B. We note that there are multiple instances in this article in which we discuss the excess chemical potential associated with moving a single bead from one solvent to another. Therefore, henceforth we use the following notation: for a single bead i, the μ_ex of moving that bead from solvent j to k is denoted as Δμⁱ_ex(j → k).

Using eqn (11) to define a_ij requires us to know how the choice of a_ij value relates to a bead's excess chemical potential. In our previous work⁵⁸ we calculated this via a Widom insertion approach. We showed that for a single solute bead inserted into a solution of unbonded beads, the relationship between the excess chemical potential and a_ij takes the form

where we fit our calculated data to determine the constants as A_S = 2.8248 × 10⁻⁵, B_S = −6.391 × 10⁻³ and C_S = 0.6396, D_S = 0.12443 and E_S = 13.843.

In this study, we will also be modelling a solvent of dodecane, which consists of a series of four bonded beads with bond length l₀ = 0.6 (details of this solvent are provided in more detail in later sections). In our previous work,⁵⁸ we showed that the expression given in eqn (12) can alter as a result of changes to the solvent. Therefore, we also perform Widom insertion calculations with this alternative bonded solvent. Here, we find a slightly different relationship of the form

where A_B = 3.0234 × 10⁻⁵, B_B = −6.4998 × 10⁻³ and C_B = 0.6350, D_B = 0.13624 and E_B = 13.979. Fig. 1 shows the difference between the two relationships, highlighting that while the difference is negligible for very low repulsion values, the difference becomes significant at large a_ij.

Fig. 1 Excess chemical potential μ_ex calculated via Widom insertion for beads with different a_ij interaction values, where a_ij is the interaction between the solute j and solvent i. The relation is shown for two different solvents, where one is bonded (red), and the other is unbonded (black), where the a_ij values for the solvent beads are a_ii = 29.4 and a_ii = 25, respectively.

3 Parameterisation

In this section, we first discuss the coarse-grained mapping chosen for the molecules simulated, then discuss the determination of a_ij for the various bead types used, and finally, the choice of parameters defining bond and angle potentials.

All a_ij bead parameters are determined via matching to experimental data for studies performed at room temperature. The choice here to match experimental values, rather than alternative approaches such as using values calculated using COSMO-RS, is related to a known difficulty using COSMO-based approaches for parametrising EO groups.^51,59

3.1 Coarse-graining and unit conversion

In this paper, we model surfactants at a dodecane–water interface. We model multiple types of surfactant, including: nonionic polyoxyethylene alkyl ether surfactants H(CH₂)_i(OCH₂CH₂)_jOH (often abbreviated as C_iE_j), the anionic surfactant sodium dodecyl sulfate (SDS), and zwitterionic surfactant dodecyldimethylamine oxide (DDAO).

To represent the dodecane molecules, we model this as four bonded beads of the same type, such that each bead represents 3 carbon atoms. Water is modelled as a single bead representing a number of water molecules (see below)). We choose to model dodecane using all one bead type (as opposed to defining the tail beads separately) as this greatly simplifies setting up the solvent, and it is unlikely to have much impact on the behaviour.

For the nonionic surfactant, we coarse-grain slightly differently, such that we have four bead types, which are illustrated in Fig. 2, and are labelled C3T, C3, EO, EOT, i.e. in this case we treat the end beads as behaving differently to those contained within the tail or head chain. For SDS and DDAO, we treat the tail in the same way, also shown in Fig. 2. For SDS, the head group is represented as a single bead with a net negative charge (−e), and the sodium counter ion is treated as a single bead with a net positive charge (+e). Similarly, the head group of DDAO is represented as a single, neutral bead.

Fig. 2 Coarse-graining used for the surfactant molecules modelled in this work. Colours represent the different bead types used.

As previously highlighted, DPD simulations are conducted in reduced units (typically referred to in literature as ‘DPD units’) such that distances (e.g. simulation box dimensions) are presented in units of the cut-off R_C. It is possible to determine a value for R_C in ‘real’ units to relate our simulations to a real length scale. We do this by matching the mass density in our simulations to the experimentally known mass density at room temperature. Typically it is chosen to match to the density of the water, however, here we take a slightly different approach.

We choose to set up our simulations such that the bead density of the coexisting dodecane and water phases are both ρR_C³ = 3. Dodecane is defined as 4 bonded DPD beads. This means that we can define R_C in real units as that which reproduces the correct density of dodecane at room temperature. Taking the density of dodecane to be 750 kg m⁻³, we find R_C = 6.56 Å.

We now turn to the density of the water phase. There is no theoretical reason why a single water bead needs to represent an integer number of water molecules. Therefore, since the density of the water phase is related to the number of water molecules N_W a single bead represents (as well as R_C), we set the level of water coarse-graining to reproduce the density of water at room temperature with the previously determined R_C value. We find the coarse-graining of water to be N_W = 3.14.

We can now use our value for R_C to choose a bond length (l₀) for dodecane (eqn (9)). Given that C–C–C angles are tetrahedral (approximately 109.5°) and that the experimental C–C bond is 1.543 Å,⁶⁰ the separation between three carbon atoms ≈3.78 Å. Therefore we choose to set an equilibrium bond length of 0.6R_C for the beads representing dodecane in our simulation. Since our surfactant beads shown in Fig. 2 are defined with a very similar degree of coarse-graining, we choose the same equilibrium bond length for bonding between surfactant beads.

Having defined the coarse-graining we are now also able to define a conversion into real units for the interfacial tension calculated. Using that, at room temperature, k_BT = 4.11 × 10⁻²¹ J, the interfacial tension in DPD units γ_Sim can be converted into real units by using γ = γ_Sim × k_BT/R_C² ≈ γ_Sim × 9.55 mN m⁻¹.

3.2 Self-interactions

In our previous work,⁵⁸ we showed that bonding can have an effect on the density of a solvent, which must be corrected for when setting the self-interaction to produce the correct partitioning. Therefore, we define our self-interactions based on choosing the values of a_ii to reproduce the pressure of single-bonded beads, and hence producing the correct density. For our unbonded water beads, we set a_W,W = 25, which has a pressure of P = 23.7. Since we define dodecane using 4 bonded beads with a long length of l₀ = 0.6, we also showed⁵⁸ that this requires us to set a_Do,Do = 29.4. Since the surfactants we simulate in this work vary in length, we vary the self-interactions of different surfactants to reflect this. However, for a given bead interaction, in all surfactants we aim for the same change in chemical potential Δμ^ex. This results in slightly different a_ii values for different C_iE_j molecules.

The relationship between surfactant length and the self-interaction is provided in the SI. We note here that we choose to simplify, by defining the interaction a_C3,C3T to be the same as the self-interaction, i.e. a_C3,C3T ≡ a_C3,C3 ≡ a_C3T,C3T. We also make a similar simplification for the head group interactions, saying that a_EOT,EO ≡ a_EO,EO ≡ a_EOT,EOT.

3.3 Hydrocarbon/water cross-interactions

To set the interaction of the surfactant tail beads with water, and also dodecane beads with water, we use experimental IDACs. Theoretically, if one is using IDACs to define an interaction, for example between beads i and beads j, one could choose from the IDAC of solute i in solvent j, or could also choose solute j at infinite dilution in solvent i. It typically makes sense to choose the bead which has greater abundance in the simulation to be the solvent. Given that we expect that the surfactant molecules will partition into the water and dodecane phases, it makes sense to define the interaction between the tail beads and the water using the infinite dilution of alkanes in water, rather than the other way around. However, for the water/dodecane interaction, the optimal choice is less clear, which will be discussed in this section.

The experimental infinite dilution activity coefficient of dodecane in water is reported as⁵⁴ ln γ^∞ = 21.66. Taking into account the difference in the molar volume of dodecane and water, this corresponds to Δμ_ex = 19.12. Similarly, for hexane at infinite dilution in water⁵⁴ ln γ^∞ = 12.87, corresponding to Δμ_ex = 10.89. Therefore we find Δμ^C3T_ex(C3T → W) = 5.445 and Δμ^C3_ex(C3 → W) = 4.115. Using these values, along with the expression given in eqn (12) we can calculate the a_C3,W and a_C3T,W interactions. However, due to the fact that we plan to vary the self-interaction as a function of the length of the surfactant molecule, the exact values of a_C3T,W and a_C3,W will differ between different surfactants.

The cross interaction between the water and dodecane could be set with either the IDAC of water in oil or vice versa. We test the IFT calculated using both approaches, to determine the most appropriate choice. From the experimental IDAC of water in alkanes of varying length,⁵⁵ we find a value of water in dodecane to be ln γ^∞ = 6.9 for a single water molecule. This means that we require our water beads to produce Δμ^W_ex(W → Do) = 28.8, requiring us to set a_W,Do = 198.2. It is clear that this is significantly different from the values calculated earlier for the infinite dilution of hydrocarbon in water. We calculate the IFT using both interaction values to compare (see Section 5.2), where we show that a_W,Do = 198.2 produces an IFT value, which is much more comparable with experimental measurements for the IFT of the dodecane/water interface. Therefore, we choose to use a_W,Do = 198.2 for the water/dodecane interaction value. Even though the tail and dodecane beads interact differently with water, we still choose them to interact with each other through the same interaction: i.e. a_C3T,Do ≡ a_C3,Do ≡ a_C3T,C3T ≡ a_C3,C3 ≡ a_C3T,C3.

3.4 EO and EOT cross-interactions

We use a combination of different types of experimental data, in order to define the remaining interactions for the nonionic surfactants: EOT/W, EO/W, EOT/C3, EOT/C3T, EO/C3, EO/C3T, EOT/dodecane, EO/dodecane. The first simplification we choose to make is to assume that the head groups interact with any type of hydrocarbon bead equivalently, that is: a_EOT,C3 ≡ a_EOT,C3T ≡ a_EOT,Do, as well as a_EO,C3 ≡ a_EO,C3T ≡ a_EO,Do. This now leaves us with just four a_ij interactions to determine: EOT/W, EO/W, EOT/dodecane, EO/dodecane.

There is a distinct lack of experimental data available for parametrising the interaction of polyethene glycol beads, as discussed in existing literature,^50,51 leading to a variety of approaches from different authors. Since polyethylene glycol (PEG) is extremely soluble in water^61,62 for various PEG chain lengths, this cannot be used to define the EO/W and EOT/W interactions. Instead, we use experimental liquid–liquid equilibria data⁵² for heptane/diethylene glycol and heptane/triethylene glycol to estimate the EOT/dodecane and EO/dodecane interactions. Following this, we can use additional experimental partitioning studies of short-chain C_iE_j molecules⁵³ to define the EO/W and EOT/W interactions.

3.4.1 EOT/dodecane and EO/dodecane interactions. The infinite dilution activity coefficient is calculated as

ΔG^∞_mix is the free energy of mixing at infinite dilution and R is the universal gas constant. Alternatively, if the solubility of the solute x^sat_solute in the solvent is known, the equation becomes⁶³

where x^sat_solute is the mole fraction of the solute at saturation (solubility), ΔH_fusion is the enthalpy of fusion of the solute, and T_m is the melting temperature of the pure solute.

The difference in the free energy of mixing between molecule 1 and molecule 2

ΔG^∞_mix,1 − ΔG^∞_mix,2 = RT[ln γ^∞₁ − ln γ^∞₂] (16)

Using the expression for ln γ^∞:

Since T > T_m, we can reduce this expression to

Alternatively,

We use this expression in combination with experimental solubility values for diethylene glycol and triethylene glycol in heptane, where solute 1 is diethylene glycol and has x₁ = 0.0024⁶⁴ and solute 2 is triethylene glycol with x₂ = 0.0003.⁵² This leads us to calculate

μ^excess₁ − μ^excess₂ ≈ −1.74k_BT. (21)

Since these two molecules differ by 1 EO unit, we determine that the excess chemical potential of 1 EO unit in heptane is Δμ^EO_ex(EO → Do) ≈ 1.74k_BT. If we estimate that triethylene glycol is representative of 2 EOT units, we use the value of the solubility in heptane and eqn (17) to estimate Δμ^EOT_ex(EOT → Do) ≈ 4.1k_BT.

3.4.2 EOT/water and EO/water interactions. As already stated, it is difficult to obtain experimental IDAC values for polyethene glycols (PEG) in water from existing literature, although the reverse values have been reported (water in PEG), this would be inappropriate for our simulations.

Experimental partitioning studies of short-chain C_iE_j molecules⁵³ report that the transfer free energy of moving an ethylene oxide unit from the water phase to the dodecane phase is 0.655 kcal mol⁻¹, corresponding to Δμ^EO_ex(W → Do) = 1.1k_BT. Similarly, the cost of moving EOT from water to dodecane is Δμ^EOT_ex(W → Do) = 9.06k_BT. In combination with the μ_ex values determined in Section 3.4.1, we calculate Δμ^EO_ex(EO → W) = 0.64k_BT and Δμ^EOT_ex(EOT → W) = −4.96k_BT.

A summary of all interactions used in the simulation of C_iE_j surfactants in oil and water is given in Table 1.

Table 1 Method of defining parameters for simulating C_iE_j surfactants partitioning into oil and water phases

	Dodecane	Water	Tail	EO (head)	EOT (head)
Dodecane	Defined to reproduce single bead pressure	—	—	—	—
Water	IDAC of water in dodecane	Defined as 25.0 (single beads)	—	—	—
Tail	Set to produce Δμ^C3/C3Tex (C3/C3T → Do) = 0	IDAC of dodecane in water	Defined to reproduce single bead pressure	—	—
EO (head)	IDAC of diethylene glycol and triethylene glycol in heptane	Transfer free energy (water to dodecane) of EO from CiEj surfactants	Equivalent to EO/dodecane interaction	Defined to reproduce single bead pressure	—
EOT (head)	IDAC of triethylene glycol in heptane	Transfer free energy (water to dodecane) of EOT from CiEj surfactants	Equivalent to EOT/dodecane interaction	Self interaction of EO and EOT (which are equivalent)	Defined to reproduce single bead pressure

---

3.5 SDS and DDAO cross-interaction

We now define interactions for those related to the head groups of the SDS and DDAO molecules. It has become common in DPD simulations to define the interaction of charged beads (such as the SDS head group and counter-ion) to be the same as that of water. This is due to the fact that charged particles are typically expected to be surrounded by a number of water hydration molecules when in aqueous solution. Therefore, at short range, the ions will interact in a similar manner to interactions with pure water. This is one of the principles we adopt for our SDS parameters. Similarly, interaction of the head group with the dodecane is based on the excess chemical potential of a single water molecule in the dodecane phase, such that Δμ^SO4_ex(SO4 → Do) = 9.17k_BT.

We base the interaction of the DDAO head group with water on the solubility of trimethylamine N-oxide (TMAO) in water, which has a value of 7.28 mol kg⁻¹ (ref. 65) (mole fraction x ≈ 0.116). From this, we estimate the excess chemical potential (using eqn (15)) of the head group in water, and remove some of the contribution that comes from the additional CH₃ group, to calculate Δμ^DMAO_ex(DMAO → W) = 0.8k_BT. Once again, the interaction of the head group is defined such that Δμ^TMAO_ex(TMAO → Do) = 9.17k_BT. Extra discussion about this choice for the TMAO head group is provided in the SI.

3.6 Calculating a_ij values

In the previous sections, we have described how we determine the target excess chemical potentials μ_ex for each type of bead. We also describe in Section 2.2, how it is possible to relate this target μ_ex to the value of a_ij which should be set for the DPD simulations. However, we also demonstrated in our previous work⁵⁸ that when beads are bonded together, their overlap can lead to lower excess chemical potential values than the sum of the unbonded beads individually. To correct for the overlap, we use a method described in our previous work⁵⁸ which accounts for the overlap by calculating the overlap volume each bonded bead has with its neighbours. We use an automated script to determine the a_ij values for different molecules, which can be found in the SI, along with full details on all of the exact a_ij values used for every simulation case in this work.

3.7 Bonded interactions

The final parameters required for our DPD simulations are related to the strength of the bonding and angle control. It has been shown that properties such as the CMC¹⁶ can often be dependent on the strength of the potentials defining the angles and bonds (eqn (9) and (10)), as well as the exact choices of the equilibrium bond length l₀ and equilibrium angle θ₀.

We define the equilibrium bond length in our simulation based on the theoretical bond length expected for a molecule in the trans conformation. For dodecane, this was determined in Section 3.1 as l₀ = 0.6R_C, a value that we also use for our surfactant bonds, which have a comparable degree of coarse-graining. Following our previous works,⁵⁸ we choose to set K_B = 5 and K_A = 5 in all cases.

We define the equilibrium bond angle θ₀ based on reproducing the correct angle distribution obtained from atomistic molecular dynamics simulations. This has been performed in other DPD studies, allowing us to utilise these previous works to define that for bonds containing primarily EO groups (i.e. EO–EO–EO and C3–EO–EO bonds) as having θ₀ = 135° and those containing primarily C3 groups (i.e. C3–C3–C3 and C3–C3–EO bonds) as having θ₀ = 180°. Details of this investigation can be obtained in the SI. Note that we also investigate the sensitivity of our simulations to the choices for angle parameters, where we investigate the effect of altering K_A = 5 on the calculated IFT. However, we conclude that there is minimal effect when K_A is large enough. Further details can be obtained in the SI.

4 Simulation set-up and data analysis

4.1 Set-up

The first stage of our simulations involves constructing a simulation box consisting of coexisting oil and water phases. We construct these two phases by populating half of the simulation box with dodecane molecules and the other with water molecules, generating random initial placements within these domains. We then initialise the surfactants at each interface by arranging the surfactants with their tails primarily in the oil phase, and the head groups in the water phase. The in-plane placement of molecules is generated at random. An example of an initial configuration is illustrated in Fig. 3, where the simulation box has periodic boundary conditions in every direction.

Fig. 3 An example of an initial configuration for our simulations (C₁₂E₈ surfactants). Note that in the lower figure we have removed the oil and water molecules for clarity of observing the surfactant layers. Visualisation created using VMD.⁶⁶

Our choice of simulation box size is related to the fact that C_iE_j surfactants are expected to form relatively large micelles, ranging from diameters of around 4–7 nm⁶⁷ for the different i and j used in this work. For C₆E_j surfactants we use a box with dimensions 60R_C × 20R_C × 20R_C, while for C₁₂E_j surfactants we use 160R_C × 40R_C × 40R_C. The larger box size for C₁₂E_j surfactants is a reflection of their larger expected micelle sizes. However, we discuss the effects of the box size (both larger and smaller) in detail in Section 5.2.

4.2 Simulation analysis

In this section, we define various parameters that will be used to study the results of our simulations. These parameters include: interfacial tension, partition coefficient, surfactant monolayer thickness, and the maximum packing at the interface (area per molecule).

4.2.1 Interfacial tension. We construct our simulation box such that L_x = L_z and L_x > L_y, such that we form our interfaces in the y–z plane. The interfacial tension γ can be calculated using the following definitionwhere p_N(x) = p_xx(x) is the local normal pressure and p_T(x) = p_yy(x) = p_zz(x) is the local tangential pressure. This approach requires calculation of the local values of pressure tensor components for integration over the full box length in the x-dimension. Alternatively, for a system with two planar interfaces, the surface tension is related to the macroscopic averageswhere P_N and P_T are macroscopic normal and tangential components (average pressure component over the entire box) and L_x is the length of the simulation box in the direction perpendicular to the interfaces. Note that the factor accounts for the two interfaces in the simulation. We calculate the tangential pressure component as an average of the components in the y and z directions.

In Section 3.3 we discussed that, theoretically, one could define the oil/water interaction based on the IDAC of alkane in water or vice versa. Based on using the IDAC of water in alkane, we calculated an interaction value of a_W,Do = 198.2. In the absence of surfactants, we calculate the IFT of a pure water/dodecane system, where using a_W,Do = 198.2 leads to an interfacial tension of 50.1 mN m⁻¹. Experimentally, the interfacial tension of water/dodecane at 25° is 52.55 mN m⁻¹,⁶⁸ meaning that this choice for a_W,Do provides a good match to experimental results. Alternatively, had we used the IDAC of hydrocarbon in water, this would have led to an interaction value of a_W,Do = 60.0. For comparison, we calculate the IFT once again for a pure oil/dodecane system. In this case we find an IFT of 26.1 mN m⁻¹, which is significantly lower than the experimental value. Therefore we conclude that the oil/water interaction is best represented using a_W,Do = 198.2, and use this for the remainder of this work.

4.2.2 Partition coefficient. Experiments studying C_iE_j molecules⁵³ show that C_iE_j molecules can move into both the oil and water phases. The partition coefficient can be used to quantify the separation of the surfactant into the two phases. The partition coefficient is defined as Kwhere c^O and c^W are the concentrations in the oil and water phases, respectively. In our previous work,⁵⁸ we showed that we can express the expected partition coefficient K (at infinite dilution) in terms of the difference between the excess chemical potential of the partitioning molecule in the two phases

4.3 Surfactant layer thickness

We define the interfacial thickness of the surfactant monolayers by fitting a Gaussian distribution to the density profile. This profile takes the form

n(z) = n_i exp(−4(z − ζ)²/σ²) (26)

where σ is the full width at 1/e of the maximum height, and ζ is the position of the centre of the peak. The thickness can then be characterised using the fitted σ parameter.

4.3.1 Area per molecule. We can calculate the average area per molecule A_S at the interface by simply calculatingwhere n_S is the number of surfactants remaining at the interface after the equilibration period, and A is the area of the interface in the y–z plane.

We show that in most of our simulations, when the interface is overpacked, the interface will ‘remove’ surfactants to achieve an optimal surface coverage (optimal equilibration period, and A). Note that this calculation is only accurate for interfaces that remain relatively planar, with minimal undulation. This will be discussed in detail in our results sections.

5 Nonionic surfactants: C_iE_j

5.1 Equilibration

To assess whether simulations are fully equilibrated, we monitor various key parameters as a function of simulation time, observing: the partitioning, area per molecule at the interface, and interfacial tension (defined in Section 4.2). Once these parameters are no longer changing, indicating that equilibration is complete, we begin collecting data for analysis.

One of the most notable behaviours during the equilibration period is that the surface is reluctant to remove the longer tail (i = 12) surfactants into the bulk water. In these cases, surfactants are only removed into the water in the form of micelles, and therefore the interface must be monitored to ensure full detachment of micelles is completed (though we note that a small number of free monomers can be removed into the bulk oil). For example, Fig. 4 provides an example of how the interfacial concentration varies as a function of simulation time for C₁₂E₄ surfactants compared with C₆E₄. We observe discrete jumps in the surface concentration as the micelles are removed from each interface. By comparison, the C₆E_j surfactants remove multiple, relatively small micelles and monomers from the interface into the bulk water, and equilibration is more rapid than for the longer tail surfactants.

Fig. 4 Equilibration of C₁₂E₄ and C₆E₄ surfactants with an initial surface density of 3.2 molecules per nm². Due to the different box sizes used this corresponds to 550 molecules per interface in the C₆E₄ case and 2200 for C₁₂E₄.

Another behaviour of interest is observed at high initial surface concentrations: some of the dodecane molecules become solubilised into the water phase, due to entrapment in the micelle core. This is illustrated in Fig. 5. The uptake of oil into the micellar cores is not surprising, having been the suggested behaviour in small-angle neutron scattering experiments²⁵ of such systems.

Fig. 5 Slice of a simulation in the x–y plane for C₁₂E₆ surfactants with an initial surface concentration of 3.5 molecules per nm². The slice is chosen to cut through the centre of the micelles in the bulk to show the oil contained in the core. Beads are coloured according to type, where: tail groups (red), head groups (orange), oil (green), and water (blue). Visualisation created using VMD.⁶⁶

5.2 Interfacial tension and maximum packing

As previously discussed, we initialise our simulations by placing surfactants at the interface, and the surfactants are free to move into the bulk oil and water during the equilibration period. In particular, we observe that ‘over-packing’ the dodecane/water interface can lead to micelle formation in the water phase. This is in contrast to traditional molecular dynamics simulations,⁹ where such behaviour can rarely be observed on the length and time scales of atomistic simulations.

In this section, we determine the number of surfactants remaining at the interface after the equilibration period, with the aim of determining the maximum packing of surfactants at the CMC. We also calculate the interfacial tension as a function of the surface coverage. Fig. 6 shows the equilibrated surface coverage as a function of the initial surface concentration. For the C₆E_j surfactants, there is a clear maximum surface concentration Γ, above which the surfactants remove themselves into the bulk water. The surface consistently equilibrates to the same surface concentration, despite an increase in the initial surface concentration Γ₀. Therefore we can calculate the maximum packing of the surface Γ by averaging the plateau region of the plots. This maximum packing concentration, along with the corresponding area per surfactant, is presented in Table 2.

Fig. 6 The final surface concentration Γ of surfactants at the oil–water interface after an equilibration period, as a function of the initial, preassembled surface concentration Γ₀. Note the box size used for C₆E_j surfactants is 60R_C × 20R_C × 20R_C, while for C₁₂E_j surfactants is 160R_C × 40R_C × 40R_C.

Table 2 The final surface concentration Γ of C₁₂E_j surfactants at the oil–water interface, and the corresponding area per surfactant

Surfactant	Γ (surfactants per nm^2)	Area (Å^2 per surfactant)
C6E2	2.24	44.6
C6E4	1.84	54.4
C6E6	1.73	57.6
C12E4	1.86	53.9
C12E6	1.73	58.0
C12E8	1.63	61.2

---

For the longer tail C₁₂E_j surfactants, the behaviour is more complex. The surfactant C₁₂E₂ does not remove any molecules from the interface at all, regardless of the initial surface concentration. This results in highly warped and oscillating interfaces when the surface concentration is high, and the surface is clearly overpacked. This is illustrated in Fig. 7. This behaviour is likely to be because C₁₂E₂ surfactants are incapable of forming micelles, and therefore the surface chooses to overpack rather than remove individual molecules. This is in agreement with experimental data, where C₁₂E₂ does not form micelles⁶⁹ at room temperature. Due to our finite sized simulation box, our simulations are incapable of forming the liquid crystalline phases in the bulk as would be experimentally expected.

Fig. 7 An illustration of the distortion exhibited by overpacked C₁₂E₂ interfaces at two initial concentrations: Γ₀ = 2.3 molecules per nm² (top) and Γ₀ = 3.5 molecules per nm² (bottom). Note that we have removed the oil and water molecules for clarity. Visualisation created using VMD.⁶⁶

Surfactants with longer head portions (j = 4, 6, and 8), all remove micelles from the interface when the surface concentration is high enough. It is worth noting that we sometimes observe that one surface will remove a larger micelle than the other, and the values presented in Fig. 6 and Table 2 are averages of the two surfaces in the simulation. Uneven removal is mostly common for longer surfactants and will be discussed again in Section 5.3 when discussing the monolayer thicknesses.

There exist, however, some intermediate concentrations, which we propose are unnaturally ‘overpacked’ (e.g. C₁₂E₄ at approximately Γ = 2.6 molecules per nm²). That is there appears to be a critical concentration point, only above which micelles are removed into the bulk. Following this ‘tipping point’, we once again observe that the final surface concentration plateaus, and from this the maximum packing can be determined. At these intermediate concentrations, the surfaces are highly nonplanar showing significant oscillation, which is in agreement with the suggestion that they have surpassed their optimal packing concentration.

We suggest that surfaces overpack themselves at intermediate surface concentrations because there are not enough excess surfactants to remove a whole micelle into the bulk water. This is in agreement with our previous observation that the surfaces seem reluctant to remove individual surfactants into the water. Therefore, the simulation chooses to have an overpacked surface, as opposed to an underpacked surface and a micelle in the bulk. This hypothesis can be confirmed simply by using a larger interface, which will be discussed in Section 5.4.

The IFT and packing of the nonionic surfactants have been studied experimentally as a function of the length of the tail and head portions of the surfactant, as well as the type of alkane oil used. Sottmann et al.²⁴ showed the area per surfactant at the water–oil interface is independent of the length of the alkyl chain of the surfactant, which is in good agreement with our simulated values presented in Table 2. Sottmann et al.²⁴ also suggest that the area per surfactant can be described by experimental fitting as a_c = 29.3 + 6.2j (Ų), where j is the length of the head group. We note that this fit is obtained primarily by fitting surfactants with shorter head groups (j ≤ 5), and also note that there is quite a lot of variation between different studies^18,70–72 owing to the difficulty in determining the packing of surfactants at the interface. However, our simulated results are in agreement with the conclusion that the area per surfactant increases with increasing head group length, and the absolute values of area are roughly in agreement with experimental data.^{18,24,70–72} Finally, we note here that our results are in excellent agreement with molecular dynamics simulations reporting the IFT for C₁₂E₆ at increasing surface concentration at the dodecane/water interface.¹⁷

The minimum IFT achievable is extremely difficult to accurately measure experimentally, owing to the fact that for C_iE_j surfactants it is typically measured to be extremely close to zero.^18,70,72,73 Fig. 8 shows the interfacial tension calculated in our simulations, as a function of the initial surface concentration. We observe that there is a drop in the IFT with increasing surface concentration, until the IFT plateaus at the point it reaches its maximum surface coverage and begins to remove surfactants from the interface. This minimum IFT is lower for the C₁₂E_j surfactants than it is for the C₆E_j surfactants. In each case, the plateau of the interface is used to calculate a value for the IFT at maximum packing, and these values are given in Table 3.

Fig. 8 The final surface concentration Γ of surfactants at the oil–water interface after an equilibration period, as a function of the initial, preassembled surface concentration Γ₀. Note the box size used for C₆E_j surfactants is 60R_C × 20R_C × 20R_C, while for C₁₂E_j surfactants is 160R_C × 40R_C × 40R_C.

Table 3 Minimum IFT for various surfactants as determined from Fig. 8. Note, we see some small fluctuations in IFT within the plateau regions of the curves in Fig. 8 of the order of 1 mN m⁻¹

Surfactant	IFT (mN m^−1)
C6E2	2.79
C6E4	4.70
C6E6	2.82
C12E2	−0.29
C12E4	0.60
C12E6	−0.42
C12E8	−0.35

---

It is important to note that for some of the data points in Fig. 8 the use of eqn (23) is strictly speaking not valid, due to the fact that the interface is not always planar. However, we only use IFT data corresponding to the plateaus in Fig. 6 to calculate the IFT values presented in Table 3 (i.e. the region in which we expect the interface to be planar). An exception to this is for C₁₂E₂ where we observe no plateau, so we use IFT for data with initial surface concentrations Γ₀ ≥ 2.5 molecules per nm².

We observe in Table 3 that several of the values are below 0, and thus represent unphysical IFT values. For C₁₂E₂ we have discussed how the data points used do not correspond to planar interface values, so this is not unexpected. For C₁₂E₆ and C₁₂E₈, the actual IFT is likely to be extremely close to zero, and we attribute the fact that γ_CMC < 0 to the fact that there is a reasonable amount of uncertainty in our measurement due to small fluctuations in the predicted optimal surface coverage at CMC. As already highlighted there is a degree of fluctuation in the size of the micelle removed from the surface, which could explain the slight dip below zero in the IFT.

5.3 Surfactant layer thickness

There is limited experimental data available detailing the thickness of surfactant layers at oil–water interfaces. When surfactants adsorb at air–water interfaces, the interfacial thickness is typically measured by small-angle neutron scattering measurements.^74,75 However, this is more difficult for oil–water interfaces due to difficulties in creating a thin enough layer of either the oil or water for penetration of the neutron beam.^25,76 Therefore, simulations can be of significant use in studying this property.

Examples of density variation across the simulation box are illustrated in the SI (Fig. S5). We calculate the surfactant layer thickness in each simulation case, as defined using eqn (26), and is shown in Fig. 9. Note that when the surface is significantly overpacked (see Fig. S5b), the surface becomes non-planar and fitting a Gaussian distribution to the surfactant density at the surface becomes impossible. Data points where a Gaussian distribution could not be adequately fitted are omitted. We also exclude from this plot data points which are calculated for when we know the surface is overpacked, and thus isn’t expected to be a realistic measurement, however, these values are included in the SI for reader interest.

Fig. 9 Calculated layer thickness as a function of the initial surface concentration. Note that certain data points are omitted as explained in the main text. At the lower end of concentration, the relation is fit using a 2nd order polynomial (where the exact fitted functions are available in the SI), and the upper plateau is used to calculate the layer thickness at CMC.

We observe an increase in the thickness of the layer with an increasing number of surfactants at the interface. At lower concentrations there is a polynomal trend in the thickness vs. concentration relationship. Therefore we note that it is likely possible for one to determine the optimal maximum packing of the interfacial layer by studying the point at which the thickness of the layer deviates from the trends shown in Fig. 9.

When we increase the number of surfactants initially placed at the interface, the surface density plateaus following the removal of surfactants into the bulk, and therefore, we observe that there is also an approximate plateau in the thickness of the interface. The thickness of the interface at maximum packing is presented in Table 4.

Table 4 Average thickness of surfactant layers at the CMC as determined from Fig. 9

Surfactant	Thickness (nm)
C6E2	1.80
C6E4	1.85
C6E6	2.12
C12E4	3.22
C12E6	3.53
C12E8	3.92

---

We note that this plateau is fairly clear for the shorter surfactants (i = 6), but is less so for those that are longer (i = 12). We attribute this partially to the uneven manner in which surfactants are removed from the interface as the surfactant increases in length, as highlighted in Section 5.2. It is not always the case that the two interfaces in the system equilibrate to the same packing. The values presented in Fig. 9 are an average of the two interfaces, however the width can vary significantly between the two. For example, the final surface concentration of the two interfaces for C₁₂E₈ (when the initial Γ₀ = 3.2 molecules per nm²) are 1.58 and 1.75 molecules per nm². However, this results in corresponding thicknesses of 3.2 and 6.3, respectively. We believe this unevenness is more reflected in the plots of the interface thickness compared with the maximum packing, because of the nonlinear relationship between the concentration and the thickness.

Finally, we note that we observe the hydrocarbon tails of the surfactants are perpendicularly oriented to the interface, in agreement with other works.⁷⁷ This is in contrast to our observations of surfactants at air/water interfaces in our previous work,⁷⁸ and this orientational change is reflected in the increased thickness of the interfacial layer when compared to air/water interfaces.

5.4 Effect of box size on IFT, maximum packing, and layer thickness

It was discussed in Section 5.2 that we expect the interface size of 40R_C × 40R_C is not large enough to allow micelle removal in some cases. For our C₆E_j surfactants, we used an interface size of 20R_C × 20R_C, which proved to be sufficiently large for these surfactants given their smaller micelle sizes. For the C₁₂E_j surfactants, we chose to use a larger interface of 40R_C × 40R_C, due to our correct prediction that the micelles would be larger. However, as the simulations seem unable to remove surfactants from the interface unless it can produce a minimum micelle size, we now test the effect of a large interface on one of our simulation cases: C₁₂E₄ surfactants at a surface coverage of 2.6 molecules per nm². From our existing simulations in Fig. 6, we can estimate the minimum micelle size for this surfactant type, and this leads us to estimate that a box size 160R_C × 60R_C × 60R_C should lead to a large enough interface size for micelle removal.

Fig. 10 shows micelle removal from the interface into the bulk at this larger interface, confirming our hypothesis about the effects of box size. For comparison, the results from the original box dimensions are shown in Fig. S5, where the surface morphs instead of removing micelles into the bulk water. In this larger simulation, we calculate an IFT value of 0.07 mN m⁻¹ and a final surface coverage of 1.88 molecules per nm², which is in excellent agreement with what is predicted from our results in Table 2.

Fig. 10 C₁₂E₄ surfactants with an initial surface coverage 2.6 molecules per nm², with an increased box size of 160R_C × 60R_C × 60R_C. Visualisation created using VMD.⁶⁶

We further investigate the effects of box size on our simulations of C₁₂E_j surfactants by performing a variety of simulations at the reduced box size 60R_C × 20R_C × 20R_C. We find that we observe limited/no surfactant removal into the bulk in these smaller simulations, due to the effects discussed in this section. However, we find that the IFT calculated at a particular surface coverage is independent of the box size. This means that when there is no expected micelle removal from the interface (due to the simulations being initialised at a surface concentration which is below the expected maximum packing) there is no effect on the IFT. In summary, box size effects are only observed due to the different behaviours associated with micelle removal. More details about this study can be found in the SI.

5.5 Partitioning

Partitioning experiments on short-chain C_iE_j molecules⁵³ show that C_iE_j molecules can move into both the oil and water phases, with a partition coefficient K which depends on the length of the head (j) and tail (i) portions of the surfactant. For short-chain surfactants, there is normally a preference to partition into the water, while as the length of the hydrocarbon portion grows the preference shifts to the oil phase (see Table 5). It is important to note that the experiments highlighted are conducted in the pre-cmc regime. At higher concentrations, surfactants are capable of forming micelles (or reverse micelles), and the partition coefficient will alter due to this aggregation. Micelle formation results in the shielding of the surfactant tails, making the water phase more attractive than when the monomers are free.

Table 5 A selection of experimental partitioning values K (where K is defined in eqn (24)). Note that values labelled (a) are taken for water/hexadecane and (b) for water/hexane. All other values taken for water/dodecane systems

Surfactant	Ln K	Surfactant	Ln K
C4E2	−2.70^53	C12E2	≈8^a^18
C5E2	−1.39^53	C12E3	≈7^a^18
C4E4	−4.91^53	C12E4	≈6^a^18
C5E4	−3.56^53	C12E5	6.0^b^70
C10E4	3.47^26	C12E6	3.42^26

---

This is behaviour that we observe in our simulations. For example Fig. 11 shows the concentrations of surfactants in the bulk oil and water phases as a function of the initial surface concentration. We observe that at high enough system concentration, single surfactant molecules will remove themselves into the bulk oil, which is expected to be the preferred phase based off the partitioning values shown in Table 5. However, when the surface concentration is high enough that micelles can form in the water, there is a dramatic jump in the concentration in the water phase. To calculate an exact value for the partition coefficient to compare with experiments, our simulations should have concentrations in the bulk water and oil phases that are below the CMC concentration. However, for our surfactants with longer hydrocarbon chains, this becomes unrealistic on the length scale of our simulations. For example, surfactant C₆E₄ has an experimental CMC of 109 mM⁷⁹ while for C₁₂E₄ the CMC is several orders of magnitudes smaller, at around 0.1 mM.⁸⁰ This concentration is unobtainable, and for example, the removal of just a single surfactant into the water bulk would generally put the surfactant concentration above the CMC for all of the C₁₂E_j surfactants modelled in this work.

Fig. 11 The partitioning behaviour of C₁₂E₄ surfactants between the oil and water phases. Note that concentration is plotted using the ‘symlog’ scale in Matplotlib that combines both linear and logarithmic scaling, behaving linearly for small values (including zeros) and transitioning to a logarithmic scale for larger values.

To calculate the partition coefficient, we also need to have a suitable number of surfactants in both the water and oil phases for a precise calculation (and as stated, without having a concentration in the water phase which is above the CMC). Due to the difficulties in achieving bulk concentrations in this narrow range of concentrations (due to our finite simulation box), we approach the calculation of partitioning slightly differently from our previous simulations.

Eqn (25) highlights that the partitioning at infinite dilution is only dependent on the excess chemical potential of the surfactant molecule in the water phase and the oil phase. This means the partitioning is related to the interaction parameters between the surfactant beads and the oil, and the surfactant beads and the water, as opposed to the self-interaction parameters between the surfactant molecules themselves. One way to encourage more molecules to move into the bulk phase (and discourage micelle formation) is to increase the value of the a_ij parameters between the surfactant molecules. While this produces unrealistic behaviour, particularly if one is interested in calculating properties such as the maximum packing, IFT, or micellar behaviour, it should, importantly, have no impact on the partitioning behaviour, except making it easier to calculate. Therefore, this is the approach we take to determine the partitioning K of a couple of test surfactants as a means to check that their interaction with the oil and water phases is accurate.

We increase the self-interactions between the surfactants to a_ij = 100 (see SI for a list of interactions this applies to). We calculate the partitioning value K at infinite dilution by varying the number of surfactants placed at the interface of the box and extrapolating to infinite dilution to calculate K. We perform this calculation on surfactant types C₆E₄, and C₁₂E₆. These surfactants are selected as they represent two different partitioning behaviours (preference for water vs. oil, as illustrated by the experimental data in Table 5), while having small enough K values to be calculable in a simulation.

We find that increasing the self interaction values greatly increases the number of surfactant molecules in the bulk oil and water phases as expected. This allows us to calculate partition coefficients for C₆E₄ and C₁₂E₆ as ln K = −2.8 and ln K = 4.0, respectively. We extrapolate the experimental data in Table 5, to say that the expected values are ln K ≈ −2.1 and ln K = 3.4, meaning the partitioning behaviour in our simulations are in reasonable agreement with experimental expectations. We note that there is a high degree of uncertainty in our calculation of the partition coefficient, due to the relatively low concentration in each of the bulk phases (even with the changes made to the self-interaction values). Note that further details about this calculation, including the procedure for extrapolation to infinite dilution can be found in the SI.

6 Ionic surfactants: SDS and DDAO

In this section, we simulate anionic surfactant sodium dodecyl sulfate (SDS), and zwitterionic surfactant dodecyldimethylamine oxide (DDAO). We simulate each surfactant at the dodecane–water interface individually in their pure form, and also simulate a mixed system containing DDAO and SDS at a 50 : 50 molar ratio.

Experimental measurements report a synergistic effect between SDS and DDAO,^19,20 in which the CMC and minimum surface tension of mixed SDS/DDAO solutions is lower when compared to either pure surfactant solution. This effect is reported to be greatest at a surfactant ratio of 50 : 50.¹⁹ In this section, we also investigate whether these synergistic effects are reproducible in our DPD predictions for interfacial tension.

6.1 Initial testing

For our initial simulations, we choose to use a box with dimensions 160R_C × 40R_C × 40R_C, in line with our choice for nonionic surfactants in the previous section. A plot showing the interfacial concentration and IFT in this study is shown in Fig. 12. The relationship we observe between surface concentration Γ and IFT for pure SDS compares extremely well with results obtained from other DPD studies and MARTINI calculations.^44,49

Fig. 12 Results for charged surfactants in our initial test box size of 160R_C × 40R_C × 40R_C. Figures show the final surface concentration Γ (top) and IFT (bottom) as a function of the initial surface concentration Γ₀ for SDS, DDAO and a 50 : 50 mixture of SDS and DDAO.

At higher concentrations, similarly to our nonionic study, we observe that overpacking the interface with SDS molecules leads to micelle removal from the interface (into the bulk water) such that the concentration at the interface becomes independent of the initial surfactant concentration. The behaviour of overpacked interfaces with DDAO is somewhat different to SDS, as the DDAO molecules prefer to move into the oil phase as opposed to the water. We wish to note that we confirm that this behaviour is not because DDAO is incapable of forming micelles (see SI for more information), but rather, is because DDAO has a lower excess chemical potential in the oil phase as opposed to the water phase. Both of these partitioning behaviours are illustrated in Fig. 13.

Fig. 13 An illustration of the movement of DDAO (top) and SDS (bottom) molecules into the bulk water and oil phases. Note that we have removed the oil and water molecules for clarity, however the oil and water placement is the same as that shown in Fig. 3. Beads are coloured according to type: tail (red), DDAO head (blue), SDS head (green) and SDS counterion (purple). Visualisation created using VMD.⁶⁶

Similarly to previously, we obtain the minimum IFT and maximum packing at the interface from the plateau's observed in Fig. 13, and these values are given in Table 6. The density of surfactants at the interface is higher for DDAO than SDS, as a result of no net charges on its head group, which is the limitation on maximum packing for SDS molecules. We observe that there is very little difference between the minimum IFT of pure SDS and the mixed surfactant case, where both minimum IFT values are close to zero, making them difficult to distinguish. However, we now discuss the effects of simulation box size, which are found to be more significant when simulating these surfactants.

Table 6 Minimum IFT and maximum packing at the interface from the plateau's observed in Fig. 13, where simulations are conducted in simulations with dimensions 160R_C × 40R_C × 40R_C

Surfactant	Γ (surfactants per nm^2)	Area (Å^2 per surfactant)	IFT (mN m^−1)
SDS	2.08	48.1	1.20
DDAO	2.61	38.3	4.60
50 : 50 mix	2.45	40.9	0.883

---

6.2 Effects of box size

In our nonionic surfactant study, we observed that there were limited finite size effects related to the box size, as long as the initial density of the interface was high enough such that there were enough excess molecules allowing micelles to be removed from the interface. However, in this study, we observe that there are significant finite size effects for simulations of SDS, which result from the high concentration of ions in the bulk water. Similarly, DDAO also exhibits a box size dependence due to the relatively high concentration of DDAO in the oil phase. In an experimental system, the concentration of surfactants in the bulk would be relatively low owing to the significantly larger bulk volume. Therefore, we aim to extrapolate to a theoretically infinite volume, to find more accurate values for the maximum packing and IFT than those presented in Table 6.

To calculate the minimum IFT and maximum packing of surfactants at the interface, we perform a study in which we calculate γ and Γ as a function of increasing box length L_x, effectively lowering the ion/DDAO concentrations in the respective bulk phases while keeping the interface the same size. For this to be computationally feasible, we simulate in a box with dimensions 20R_C × 20R_C in the y–z plane, reducing the size of the interface such that we can increase the length of the box more easily. We focus our efforts on studying systems where the number of surfactants at the interface has exceeded maximum packing, due to the computational effort that would be involved in such a study for every Γ₀. Therefore we can calculate the maximum packing and minimum IFT, by fitting our data as a function of L_x and extrapolating L_x → ∞.

For our fitting, we propose that the surface concentration Γ ∝ N_I where N_I is the number of ions at the surface. We also propose that N_I ∝ 1/V_B, where V_B is the volume of the bulk. Therefore, our data is fit with a function of the form

where Γ is the final surface concentration, Γ_∞ is a fitted constant representing the surface concentration at infinite bulk size, A is a fitted constant, and L_x is the length of the box used in the simulation. We find our simulations exhibit an inverse relationship between the interfacial tension and surface concentration i.e. γ ∝ 1/Γ. Therefore, for the IFT at maximum packing, we fit a function wherewhere γ_∞ is a fitted constant representing the IFT at infinite bulk size, and B is a fitted constant.

The results of this study are shown in Fig. 14, where eqn (28) and (29) are shown to be reasonable expressions to represent Γ and γ as a function of box length. The maximum packing and minimum IFT value obtained from these fits are presented in Table 7. We see that the final estimate for the maximum packing Γ is lower and therefore the estimate for the lowest interfacial tension increases, compared to the results from the single box size shown in Table 6.

Fig. 14 Results for charged surfactants in simulation box size of L_xR_C × 20R_C × 20R_C, where we extrapolate to an infinite length box L_x → ∞. Figures show the final surface concentration Γ (top) and IFT (bottom) as a function of the box length L_x. The initial surface concentrations Γ₀ are different for each surfactant, and are chosen based on their expected maximum packing, where we use: 4.07 (DDAO), 2.32 (SDS), and 3.49 (DDAO:SDS mixture) molecules per nm².

Table 7 Minimum IFT and maximum packing at the interface from extrapolation to an infinite length box L_x → ∞ in Fig. 14, where simulations are conducted in simulations with dimensions L_xR_C × 20R_C × 20R_C

Surfactant	Γ (surfactants per nm^2)	Area (Å^2 per surfactant)	IFT (mN m^−1)
SDS	1.94	51.4	4.51
DDAO	2.57	38.9	11.61
50 : 50 mix	2.34	42.7	1.07

---

The maximum packing for SDS surfactants at water/dodecane interfaces is experimentally reported as an area per molecule ≈50–57 Ų.^81,82 There is comparatively less experimental data available for DDAO systems compared to SDS. However, DDAO at the hexane–water interface is reported as having a maximum packing corresponding to 43 Ų per molecule,⁸³ which we expect to be comparable to water/dodecane. Therefore, the results in Table 7 are in good agreement with the experimental data for pure systems. The maximum packing of the SDS:DDAO mixture falls somewhere in between that for pure SDS and pure DDAO. However, the minimum IFT achieved is much lower. This will be discussed in more detail in the following section.

6.3 DDAO/SDS synergistic effects

In experimental studies,¹⁹ the lower surface tension of SDS:DDAO mixtures is suggested to result from an enhanced adsorption of the surfactants at the interface, resulting in a smaller average area per molecule. This enhanced packing has been hypothesised to be due to an electrostatic and steric dilution effect, which results from positive interactions between the surfactants. However, we note that our DDAO:SDS mixed system (Fig. 15) exhibits the same synergistic behaviour, even though we do not predict greater interfacial packing of the mixed case compared to both pure cases. We see a greater packing compared to SDS but not compared to DDAO.

Fig. 15 A mixed 50 : 50 system of DDAO and SDS. Note that we have removed the oil, ion and water molecules for clarity. Beads are coloured according to type: tail (red), DDAO head (blue) and SDS head (green). Visualisation created using VMD.⁶⁶

The maximum ability of a surfactant to lower the interfacial tension can be considered to be a combination of two things: how well a single surfactant molecule lowers the IFT (or IFT lowered per surfactant at the interface); and also the maximum packing that can be achieved by that surfactant i.e. more surfactants also mean lower IFT. Therefore, we propose that we observe a greater reduction in IFT for the mixed system due to a combination of the more effective lowering of the interfacial tension by SDS, while a greater packing potential of DDAO.

We now consider the parameters that influence the change in the IFT per surfactant and also which influence the maximum packing. We concentrate on the effects of: the hydrophobicity of the head group; the self-repulsion between head groups; and electrostatic effects between the head group and the counterions.

One possibility is that SDS is better at lowering the IFT than DDAO is because its head group is more hydrophilic. This increased hydrophilicity is to be expected, if one considers the hydrogen bonding each head group exhibits with the water. The number of hydrogen bonds a molecule forms is closely related to its hydrophilicity (though not necessarily directly related), as hydrogen bonding is one of the key interactions facilitating a molecule's solubility in water. The number of hydrogen bonds formed by the SDS head group is expected to be much larger than that of DDAO, due to its significantly larger dipole moment. In bulk water, the average number of hydrogen bonds per water molecule is typically calculated (using molecular dynamics simulations) to be around 3.3.^13,14,84 By comparison, the number of hydrogen bonds per head group for SDS is calculated to be 5.5¹⁴ and for DDAO is 2.1.¹³ These estimates are in agreement with the head-water interactions we have assigned for our SDS and DDAO molecules.

We can easily check the impact of the hydrophilicity on the IFT, by simulating DDAO with a head-water interaction comparable to that of SDS–water (see SI for full details of these simulations). We find, however, that the decrease in IFT per molecule is relatively similar to DDAO with its original head-water interaction. However, we find that the maximum packing increases, and this leads to an enhanced reduction in the IFT. Therefore, the greater effectiveness of SDS for lowering the IFT per molecule is indicated to not be due to the hydrophilicity, but due to electrostatic effects between head groups and counterions.

Despite the increased effectiveness of the SDS head group for lowering the IFT, the packing of SDS is limited by the fact that the head groups are repulsive to each other. In contrast, the neutral DDAO head groups have a lower net self-repulsion, leading to increased packing ability at the surface. Combining the effects discussed leads to the suggestion that an optimal balance of SDS and DDAO at the interface results from a balance of the effectiveness of a surfactant at lowering the IFT (where SDS is most effective), while achieving a high maximum packing (where DDAO is most effective).

Conclusions and future research

In this article, we have developed DPD parametrisations for modelling surfactants at the water/dodecane interface. Our results show good agreement with the available experimental data, and allow us to provide insight into systems that can be difficult to study experimentally. For example, as well as making accurate predictions about the interfacial tension, we can easily study the packing and thickness of surfactant layers at the oil/water interface, and the phenomenon of the uptake of oil into micelles allowing the solubility of oil into the water phases. We can also easily investigate the origin of synergistic effects of mixed systems, which we show results from a balance of surfactant properties.

There are some aspects of our simulations that might be interesting for further investigation. For example, simulations have shown that the effective dielectric constant of water at interfaces is considerably lower than its bulk value,⁸⁵ which can be attributed to reduced dipolar fluctuations at the alkane/water interface in comparison with the bulk.^86,87 It is difficult to assign an accurate dielectric constant in DPD simulations because the effective dielectric constant varies across the water/dodecane interface. Therefore, as here, it is usual to simply assign a constant background dielectric which represents that of bulk water. However, we note this could potentially be a source of systematic error in these type of simulations. It would therefore be interesting to incorporate a polarisable water model, such as one of those developed in our previous work,⁸⁸ to see if this alters the behaviour of surfactants at interfaces in this study. In particular, the use of a polarisable water model may influence the behaviour of the counterions, which, as we demonstrate here, can influence the packing of SDS molecules at the interface.

A further issue to note is that our study shows that finite-size effects are non-negligible. Dealing fully with these, by moving to large simulations, is certainly possible. However, this adds computational cost, making this approach considerably more time-consuming for computational screening.

We also note that in this article, we have primarily focussed on IFT calculations and surfactant behaviour at interfaces, including surfactant removal into the bulk. However, a further study utilising this model could be performed to investigate micelle formation, including micelles in the presence of oil. It has been shown that the presence of oil can sometimes lead to smaller micelles for C₁₂E_j surfactants,⁸⁹ and this model would be perfect for such a study.

Finally, the methodology used in this work could be easily extended to other surfactants or oil systems, where one could obtain the chemical potentials required for defining the DPD parameters from experimental studies (as we have done in this work) or via another simulation method such as COSMO-RS, as employed in our previous work.⁷⁸

Author contributions

RH: conceptualization, data curation, formal analysis, investigation, methodology, software, validation, visualization, writing – original draft. CA: project administration, writing – review and editing. MD: writing – review and editing. MW: funding acquisition, project administration, supervision, writing – review and editing.

Conflicts of interest

There are no conflicts of interest to declare.

Data availability

Data supporting this article have been uploaded as part of the supplementary information (SI). Supplementary information contains the relationship between molecule length and self-interaction, a discussion of DDAO head group interactions with dodecane, an automated python parametrisation script, parameters for simulation models (intermolecular and intramolecular), density profiles for surfactants at different surface coverages and surfactant layer thicknesses, partition coefficient calculations, a discussion of the relationship between IFT and simulation box size, a discussion and calculations for DDAO micelle formation in water, an examination of the effect of head-water interactions on the behaviour of DDAO. See DOI: https://doi.org/10.1039/d6sm00083e.

Acknowledgements

The authors would like to acknowledge the support of Durham University for the use of its HPC facility, Hamilton. RH would like to acknowledge support through an EPSRC/P&G funded Prosperity Partnership grant EP/V056891/1.

Notes and references

1. M. Rieger, Surfactants in Cosmetics, CRC Press, 2017.
2. I. Kralova and J. Sjöblom, J. Dispersion Sci. Technol., 2009, 30, 1363–1383.
3. D. Sharma, D. Singh, G. M. Sukhbir-Singh, B. M. Karamchandani, G. K. Aseri, I. M. Banat and S. K. Satpute, Molecules, 2023, 28, 2823.
4. O. Kontogiannis, D. Selianitis, N. Lagopati, N. Pippa, S. Pispas and M. Gazouli, Pharmaceutics, 2023, 15, 501.
5. L. Rodrigues, I. M. Banat, J. Teixeira and R. Oliveira, J. Antimicrob. Chemother., 2006, 57, 609–618.
6. X. Wang, K. P. Santo and A. V. Neimark, Langmuir, 2020, 36, 14686–14698.
7. R. L. Hendrikse, C. Amador and M. R. Wilson, Soft Matter, 2024, 20, 6044–6058.
8. S. Wang, S. Yang, R. Wang, R. Tian, X. Zhang, Q. Sun and L. Liu, J. Pet. Sci. Eng., 2018, 169, 81–95.
9. J. Li, C. Amador and M. R. Wilson, Phys. Chem. Chem. Phys., 2024, 26, 12107–12120.
10. M.-T. Lee, R. Mao, A. Vishnyakov and A. V. Neimark, J. Phys. Chem. B, 2016, 120, 4980–4991.
11. A. Vishnyakov, M.-T. Lee and A. V. Neimark, J. Phys. Chem. Lett., 2013, 4, 797–802.
12. R. L. Anderson, D. J. Bray, A. Del Regno, M. A. Seaton, A. S. Ferrante and P. B. Warren, J. Chem. Theory Comput., 2018, 14, 2633–2643.
13. C. D. Lorenz, C.-M. Hsieh, C. A. Dreiss and M. J. Lawrence, Langmuir, 2011, 27, 546–553.
14. C. D. Bruce, S. Senapati, M. L. Berkowitz, L. Perera and M. D. E. Forbes, J. Phys. Chem. B, 2002, 106, 10902–10907.
15. L. Shi, N. R. Tummala and A. Striolo, Langmuir, 2010, 26, 5462–5474.
16. M.-T. Lee, A. Vishnyakov and A. V. Neimark, J. Phys. Chem. B, 2013, 117, 10304–10310.
17. M. Kanduc, C. Stubenrauch, R. Miller and E. Schneck, J. Chem. Theory Comput., 2024, 20, 1568–1578.
18. M. J. Rosen and D. S. Murphy, Langmuir, 1991, 7, 2630–2635.
19. G. Tyagi, D. Seddon, S. Khodaparast, W. N. Sharratt, E. S. Robles and J. T. Cabral, Colloids Surf., A, 2021, 618, 126414.
20. R. A. Abdel-Rahem and F. A. Al-Odail, J. Dispersion Sci. Technol., 2014, 35, 1009–1017.
21. P. Kékicheff, C. Grabielle-Madelmont and M. Ollivon, J. Colloid Interface Sci., 1989, 131, 112–132.
22. G. Montalvo, M. Valiente and E. Rodenas, Langmuir, 1996, 12, 5202–5208.
23. G. Montalvo and A. Khan, Colloid Polym. Sci., 2005, 283, 402–412.
24. T. Sottmann, R. Strey and S.-H. Chen, J. Chem. Phys., 1997, 106, 6483–6491.
25. E. Staples, J. Penfold and I. Tucker, J. Phys. Chem. B, 2000, 104, 606–614.
26. G. Catanoiu, E. Carey, S. Patil, S. Engelskirchen and C. Stubenrauch, J. Colloid Interface Sci., 2011, 355, 150–156.
27. M. Nishizawa, K. Saito and M. Sorai, J. Phys. Chem. B, 2001, 105, 2987–2992.
28. S. Müller, C. Börschig, W. Gronski, C. Schmidt and D. Roux, Langmuir, 1999, 15, 7558–7564.
29. D. J. Mitchell, G. J. T. Tiddy, L. Waring, T. Bostock and M. P. McDonald, J. Chem. Soc., Faraday Trans. 1, 1983, 79, 975–1000.
30. P. Müller, D. J. Bonthuis, R. Miller and E. Schneck, J. Phys. Chem. B, 2021, 125, 406–415.
31. M. Kanduc, J. Reed, A. Schlaich and E. Schneck, Curr. Opin. Colloid Interface Sci., 2024, 72, 101816.
32. P. Shi, H. Zhang, L. Lin, C. Song, Q. Chen and Z. Li, J. Dispersion Sci. Technol., 2018, 39, 1258–1265.
33. A. M. Luz, T. J. dos Santos, G. D. Barbosa, C. L. Camargo and F. W. Tavares, Colloids Surf., A, 2022, 651, 129627.
34. P. Yazhgur, S. Vierros, D. Hannoy, M. Sammalkorpi and A. Salonen, Langmuir, 2018, 34, 1855–1864.
35. X. He, O. Guvench, A. D. J. MacKerell and M. L. Klein, J. Phys. Chem. B, 2010, 114, 9787–9794.
36. W.-X. Shi and H.-X. Guo, J. Phys. Chem. B, 2010, 114, 6365–6376.
37. A. Maiti and S. McGrother, J. Chem. Phys., 2004, 120, 1594–1601.
38. H. Rezaei and H. Modarress, Chem. Phys. Lett., 2015, 620, 114–122.
39. R. D. Groot and P. B. Warren, J. Chem. Phys., 1997, 107, 4423–4435.
40. F. C. de Oliveira, J. M. Maia and F. W. Tavares, Colloids Surf., A, 2021, 625, 126828.
41. D. Steinmetz, B. Creton, V. Lachet, B. Rousseau and C. Nieto-Draghi, J. Chem. Theory Comput., 2018, 14, 4438–4454.
42. H. Alasiri and W. G. Chapman, J. Mol. Liq., 2017, 246, 131–139.
43. V. V. Ginzburg, K. Chang, P. K. Jog, A. B. Argenton and L. Rakesh, J. Phys. Chem. B, 2011, 115, 4654–4661.
44. M. Ndao, F. Goujon, A. Ghoufi and P. Malfreyt, Theor. Chem. Acc., 2017, 136, 1–13.
45. S. J. Marrink, H. J. Risselada, S. Yefimov, D. P. Tieleman and A. H. de Vries, J. Phys. Chem. B, 2007, 111, 7812–7824.
46. J.-C. Neyt, A. Wender, V. Lachet, A. Ghoufi and P. Malfreyt, J. Chem. Theory Comput., 2014, 10, 1887–1899.
47. M. Ndao, J. Devémy, A. Ghoufi and P. Malfreyt, J. Chem. Theory Comput., 2015, 11, 3818–3828.
48. P. Katiyar and J. K. Singh, J. Chem. Phys., 2017, 146, 204702.
49. S. Wang and R. G. Larson, Langmuir, 2015, 31, 1262–1271.
50. R. Groot and K. Rabone, Biophys. J., 2001, 81, 725–736.
51. E. Lavagnini, J. L. Cook, P. B. Warren, M. J. Williamson and C. A. Hunter, J. Phys. Chem. B, 2020, 124, 5047–5055.
52. M. A. Hughes and Y. Haoran, J. Chem. Eng. Data, 1990, 35, 467–471.
53. M. Manabe, M. Koda and K. Shirahama, Bull. Chem. Soc. Jpn., 2006, 48, 3553–3556.
54. K. Kojima, S. Zhang and T. Hiaki, Fluid Phase Equilib., 1997, 131, 145–179.
55. S. Zhang, T. Hiaki, M. Hongo and K. Kojima, Fluid Phase Equilib., 1998, 144, 97–112.
56. P. Español and P. Warren, Europhys. Lett., 1995, 30, 191–196.
57. M. González-Melchor, E. Mayoral, M. E. Velázquez and J. Alejandre, J. Chem. Phys., 2006, 125, 224107.
58. R. L. Hendrikse, C. Amador and M. R. Wilson, Phys. Chem. Chem. Phys., 2025, 27, 1554–1566.
59. Y.-C. Hung, C.-M. Hsieh, H. Machida, S.-T. Lin and Y. Shimoyama, J. Mol. Liq., 2022, 356, 118896.
60. F. H. Allen, O. Kennard, D. G. Watson, L. Brammer, A. G. Orpen and R. Taylor, J. Chem. Soc., Perkin Trans. 2, 1987, S1–S19.
61. Y.-X. Yu, J.-G. Liu and G.-H. Gao, J. Chem. Eng. Data, 2000, 45, 570–574.
62. A. Eliassi, H. Modarress and G. A. Mansoori, J. Chem. Eng. Data, 1998, 43, 719–721.
63. A. Myerson, Handbook of Industrial Crystallization, Butterworth-Heinemann, 2002, pp. 11–13.
64. I. Ashour and M. A. Abdulkarim, J. Chem. Eng. Data, 2005, 50, 924–927.
65. C. Held and G. Sadowski, Fluid Phase Equilib., 2016, 407, 224–235.
66. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38.
67. W. C. Swope, M. A. Johnston, A. I. Duff, J. L. McDonagh, R. L. Anderson, G. Alva, A. T. Tek and A. P. Maschino, J. Phys. Chem. B, 2019, 123, 1696–1707.
68. S. Zeppieri, J. Rodríguez and A. L. López de Ramos, J. Chem. Eng. Data, 2001, 46, 1086–1088.
69. M. L. Lynch, K. A. Kochvar, J. L. Burns and R. G. Laughlin, Langmuir, 2000, 16, 3537–3542.
70. R. Aveyard, B. P. Binks, S. Clark and P. D. I. Fletcher, J. Chem. Soc., Faraday Trans., 1990, 86, 3111–3115.
71. J.-W. Benjamins, K. Thuresson and T. Nylander, Langmuir, 2005, 21, 149–159.
72. R. Aveyard, N. Carr and H. Slezok, Can. J. Chem., 1985, 63, 2742–2746.
73. T. Sottmann and R. Strey, J. Chem. Phys., 1997, 106, 8606–8615.
74. J. R. Lu, T. J. Su, Z. X. Li, R. K. Thomas, E. J. Staples, I. Tucker and J. Penfold, J. Phys. Chem. B, 1997, 101, 10332–10339.
75. J. R. Lu, Z. X. Li, T. J. Su, R. K. Thomas and J. Penfold, Langmuir, 1993, 9, 2408–2416.
76. D. Petsev, Emulsions: Structure, Stability and Interactions, Academic Press, 2004.
77. J. Chanda and S. Bandyopadhyay, J. Phys. Chem. B, 2006, 110, 23482–23488.
78. R. L. Hendrikse, C. Amador and M. R. Wilson, Soft Matter, 2023, 19, 3590–3604.
79. L. Ambrosone, L. Costantino, G. D'errico and V. Vitagliano, J. Colloid Interface Sci., 1997, 190, 286–293.
80. K.-V. Schubert, R. Strey and M. Kahlweit, J. Colloid Interface Sci., 1991, 141, 21–29.
81. S. J. Rehfeld, J. Phys. Chem., 1967, 71, 738–745.
82. J. Bergfreund, S. Siegenthaler, V. Lutz-Bueno, P. Bertsch and P. Fischer, Langmuir, 2021, 37, 6722–6727.
83. S. D. Stoyanov, V. N. Paunov, H. Rehage and H. Kuhn, Phys. Chem. Chem. Phys., 2004, 6, 596–603.
84. J. Faraudo and A. Travesset, Biophys. J., 2007, 92, 2806–2818.
85. M. Dinpajooh and D. V. Matyushov, J. Mol. Liq., 2023, 376, 121400.
86. C. Qi, Z. Zhu, C. Wang and Y. Zheng, J. Phys. Chem. Lett., 2021, 12, 931–937.
87. C. D. Wick, T.-M. Chang, J. A. Slocum and O. T. Cummings, J. Phys. Chem. C, 2012, 116, 783–790.
88. R. L. Hendrikse, C. Amador and M. R. Wilson, Soft Matter, 2024, 20, 7521–7534.
89. Y. Einaga, Y. Totake and H. Matsuyama, Polym. J., 2004, 36, 971–978.

---