Interplay between surfactant self-assembly and adsorption at hydrophobic surfaces: insights from dissipative particle dynamics

ABSTRACT We study the self-organisation of aqueous surfactants in bulk phase, and the adsorption and self-organisation of the aqueous surfactants at a planar hydrophobic surface by dissipative particle dynamics. Nonionic surfactants, n-alkyl poly(ethylene oxide) C E , and water are coarse-grained into mesoscopic beads comprising 1–3 heavy atoms and two water molecules, respectively. The size of the mesoscopic beads is related to the molar volume of the underlying molecular fragments while the bead–bead interaction parameters are calibrated against the water-octanol partition coefficients. We focus on the C E , C E , C E , and C E surfactants in water that form spherical micelles in the bulk. The bulk micellization is primarily affected by the alkyl tail length which is demonstrated by an order of magnitude decrease in the critical micelle concentration when going from the aqueous C E to aqueous C E solutions. Surfactants strongly adsorb on the hydrophobic surface, adopting lying-down configurations and forming hemispheres which are in equilibrium with the spherical micelles in the bulk. In contrast to the bulk phase, the surfactant adsorption behaviour is influenced by both the alkyl tail and head chain lengths. GRAPHICAL ABSTRACT


Introduction
Surfactants, or surface-active agents, are amphiphilic molecules which tend to aggregate. Surfactant selfassembly in bulk solution results in a variety of distinct morphologies, such as spherical and cylindrical (wormlike) micelles, hexagonal and lamellar phases, or multilamellar vesicles. The surfactant self-assembly is affected by the system temperature, solvophobicity of surfactant blocks, and surfactant architecture [1]. CONTACT  Surfactants also have the ability to undergo selfassembly and adsorption at surfaces which can dramatically change the surface properties [2,3]. The aggregation and adsorption of amphiphilic molecules on surfaces has been shown to be useful in microelectronics, lithography, sensors, optics and photonics, see e.g. Refs. [4][5][6]. The type of application dictates the desired shapes and sizes of the self-assembly domains. Surfactant self-assembly at a surface can be different from that in bulk solution because interactions between the surfactants and the solid surface play a key role. Surfactants can adsorb onto the surface in various morphologies, such as hemispheres, hemicylinders, planar micelles, or self-assembly monolayers [7][8][9][10]. The morphology of adsorbed surfactant films is then a key determinant of the interfacial properties [11,12].
Nonionic surfactants, such as n-alkyl poly(ethylene oxide) C n E m , have been found to adsorb strongly on carbon-type, hydrophobic surfaces where surfactant alkyl tail groups interact primarily with the hydrophobic surfaces via van der Waals interactions. In addition, the hydrophobic surfaces exert nonspecific hydrophobic interactions between water and the surface which drives a minimisation of the water-surface interfacial area [13]. The type of self-assembly structures of C n E m surfactants at hydrophobic surface-liquid interfaces strongly depends on both the surfactant tail-group and headgroup sizes. For example, at graphite-liquid interfaces, C 9 E 3 and C 12 E 3 surfactants, which self-assemble into bilayers as a bulk lyotropic phase, form monolayers perpendicular to the graphite substrate [14,15]. On the other hand, C 12 E 5 , C 12 E 8 and C 12 E 10 surfactants, which aggregate into wormlike micelles in the bulk, form hemicylinders at the graphite-liquid interface [15]. Nonionic surfactants in monolayers adopt standing-up configurations which allow more surfactants to adsorb on the substrate. On the other hand, nonionic surfactants in hemispheres or hemicylinders tend to adopt lying-down configurations to maximise their interactions with the surface but the surfactants also occupy a larger surface area [16].
Molecular modelling has become a complementary tool to experiments to explore ways of controlling the self-assembly and adsorption of surfactants at surfaces by tuning the surfactant-surface and surfactant-surfactant interactions via the choice of the surfactant, solvent, and surface chemistry as well as the surfactant architecture [17]. However, the time and length scales for surfactant aggregation and adsorption at a surface from a bulk solution are too long to be conveniently studied by all-atom molecular dynamics [13]. Instead, researchers rely on mesoscopic modelling, such as coarse-grain molecular dynamics (CGMD) [18] or dissipative particle dynamics (DPD) [19,20]. These methods employ coarse-grain models where molecules or molecular fragments are grouped into mesoscopic beads [21]. The DPD method has been proved to be an effective and flexible modelling tool and has been used to address a wide range of mesoscale problems in and out of equilibrium [22][23][24][25][26][27][28].
Modelling of solid-liquid interfaces within DPD is not straightforward, see e.g. Refs. [29][30][31]. The soft repulsion between DPD beads cannot prevent liquid beads from penetrating the surface boundaries, and thus extra effort is needed to impose the no-slip (or partial slip) wall boundary conditions. A well-accepted approach represents a solid surface by frozen beads with a density equal to the liquid density [32]. Liquid and solid beads then preserve their separation through the use of proper reflections when liquid beads are about to cross the given position of the solid surface [33]. The reflection scheme suppresses unphysical fluctuations in the liquid density close to the surface and imposes no-slip boundary conditions, i.e. at a solid boundary, the liquid beads have zero velocity relative to the boundary [29,33].
Besides GCMD and DPD simulations, the surfactant adsorption can also be studied by theoretical methods such as density functional theory (DFT) or lattice gas theory (LGT). For example, Rosenthal and Klapp [34] employed DFT to investigate the structure formation of amphiphilic molecules at planar walls, considering neutral, hydrophilic, and hydrophobic surfaces. They focused on the competition between the surface field and the interaction-induced ordering phenomena. A similar DFT study was performed by Borówko at al. [35] who considered planar walls modified by tethered chain molecules. Yet in another theoretical study, Bock and Gubbins [36] used LGT to explore the surfactant selfassembly from aqueous solutions on solid surfaces and the role of hydrogen bonding between water and the surfactant headgroups.
In this work, we deal with the nonionic surfactants C 6 E 3 , C 6 E 4 , C 8 E 3 , and C 8 E 4 , which form spherical micelles in bulk water, and explore their adsorption and self-organisation at carbon-type, hydrophobic surfaces by DPD simulations over a wide range of surfactant concentrations. We employ a recently developed generic DPD force field, which is calibrated against the molar volume of the underlying molecular fragments and water-octanol partition coefficients, and which provides a faithful prediction of the micellar behaviour of both aqueous ionic and nonionic surfactant solutions [37,38]. The remainder of the paper is organised as follows. The methodology, including the DPD method, coarse-grain models, and simulation details together with the definition of the system observables, is outlined in Section 2. In Section 3, we then present and discuss our results of surfactant self-assembly in bulk water, and the adsorption and self-organisation at the hydrophobic surface. Finally, we present our conclusions in Section 4.

Dissipative particle dynamics
Dissipative Particle Dynamics (DPD) is a mesoscopic model that treats groups of molecules or molecular fragments as single beads; in standard DPD, as beads of the same size [39][40][41]. In this work following Anderson et al. [37], we employ coarse-grain models with different bead sizes, enhancing realism against the standard DPD model. In our DPD, bead-bead interactions are described by a soft-sphere repulsive potential where superscript * denotes DPD units. The DPD units are introduced using m 0 , kT 0 and r c as units of mass, energy, and length, respectively, where m 0 is the mass of a reference bead (typically solvent bead), k is the Boltzmann constant, T 0 is the reference temperature, and r c is the cut-off distance. In Equation (1), a * ij = a ij r c /(kT 0 ) is the maximum repulsion between the interacting beads, r * ij = r ij /r c is their separation, r ij = |r ij |, r ij = r i − r j , and R * ij = R ij /r c is the bead size parameter. The repulsion parameter a * ij is proportional to the Flory-Huggins interaction parameter, χ ij : where C FH−DPD is the proportionality constant; e.g. in the standard DPD, C FH−DPD 0.3 [41]. Equation (2) indicates that attraction between the interacting beads is treated as weaker repulsion than the repulsions between the like beads [41].
To model chain molecules such as surfactants, beads are connected into chains, typically with a simple harmonic potential which represents bonds between connected beads. In is the bond constant, r * = r/r c is the distance between the connected beads and r * 0 = r 0 /r c is the nominal bond length. Chain rigidity can be introduced using a harmonic angular potential between pairs of bonds. In Equation (4), K * a = K a /(kT 0 ) is the bending constant, θ and θ 0 are, respectively, the angle and nominal angle between the pairs of bonds. Conservative forces, F C , F bond and F bend , are then defined as a negative derivative of corresponding interaction potentials, and for example, the bead-bead conservative force becomes where e ij = r ij /r ij is the unit vector. Generally, coarse-graining (removing degrees of freedom from the system) causes a loss of friction and 'smoothing' of the free-energy landscape, leading to faster dynamics at the mesoscale level [42]. In DPD, the removed degrees of freedom can be effectively returned to the system using pair-wise dissipative and random forces respectively. In Equations (6) and (7), /m 0 and p * α are, respectively, the mass and momentum of a bead α, ζ ij is a Gaussian random number with zero mean and unit variance, t * = t kT 0 /m 0 /r c is the time step, γ * ij = γ ij r c / √ m 0 kT 0 and σ * ij are, respectively, the friction coefficient and noise amplitude for the interacting beads, and ω(r ij ) is a weighting function.
The friction coefficient and noise amplitude are related via the fluctuation-dissipation theorem [40] σ * ij 2 = 2 T * γ * ij (8) where T * = T/T 0 and T is the system temperature. The weighting function is typically given by where R * RD = R RD /r c is a cut-off distance of the dissipative and random forces. By adjusting γ ij to experimental values of self-diffusivity or viscosity, DPD can capture realistic system dynamics at the mesoscale level [43]; not relevant for this work. The dissipative and random forces act together as a thermostat, and they along with F C also guarantee local momentum conservation and ensure correct hydrodynamic behaviour.

Coarse-Grain models
To model the aqueous solutions of nonionic surfactants poly(ethylene oxide) alkyl ethers (C n E m ) in the bulk and at hydrophobic surfaces, we employ a generic force field (FF) proposed by Anderson et al. [37]. The FF was developed by a systematic top-down thermodynamic parametrisation using water-octanol partition coefficients, supplemented by water-octanol phase equilibria and pure liquid phase density data. Beads in the FF represent either two water molecules, a solvent bead, or molecular fragments comprising 1-3 heavy atoms, i.e. C and O in the case of C n E m surfactants.
First following Groot and Warren [41], the FF asserts that the density of liquid water corresponds to the bead number density ρ * = ρ r 3 c = 3 or, equivalently, the pressure P * = P r 3 c /(kT 0 ) = 23.7 ± 0.1. Then, the use of the molecular volume of liquid water (V m = 30 Å 3 [44]), number density (ρ * = 3), and two water molecules per solvent bead (N m = 2) results in r c 5.65 Å. Second, the FF constructs C n E m surfactants with a general formula CH 3 (CH 2 ) n−1 (CH 2 OCH 2 ) m CH 2 OH from connected beads comprising (i) alkyl molecular fragments [CH 2 Figure 1 shows the coarse-grain models of water and C 8 E 4 surfactant as examples. Third, the surfaces are modelled as walls of frozen CH 3 beads with ρ * = 3 which aim to mimic carbon-type, hydrophobic surfaces at the mesoscale level.
The values of R * ii for the surfactant beads are determined from their molar volume, V m , which are listed in Table 1. The FF considers that R * ii 3 ∝ V m,ii with the proportionality constant set by the water bead mapping; therefore, Then, the values of R * ij are given by the arithmetic mixing rule as R * ij = (R * ii + R * jj )/2 which is equivalent to assigning an effective radius R * i = R * ii /2 to individual beads and to assert that The values of a * ii for the water and surfactant beads are adjusted, respectively, to match the experimental compressibility of liquid water and experimental densities of various simple molecular liquids containing the surfactant groups. Finally, starting from a * ij = (a * ii + a * jj )/2, the FF adjusts the values of a * ij to fit the experimental water-octanol partition coefficients [45]. The values of the non-bonded interaction parameters, Table 1. The values of the size parameter, R * ij = R ij /r c , molar volume, V m , maximum repulsion, a * ij = a ij r c /(kT 0 ), and Flory-Huggins parameter, χ ij , for water and surfactant beads; r c 5.65 Å is the cut-off distance, T 0 = 300 K is the reference temperature, and k is the Boltzmann constant.  3 ]), 0.1 is subtracted from r * 0 = 0.39. The nominal bond length r * 0 = 0.39 corresponds to 2.52 Å which is approximately equivalent to the distance between the centres of two alkane beads, each with two heavy atoms. For bending interactions, the FF employs the bending constant K * a = 5 and the nominal angle between the pairs of bonds θ 0 = 180 • .

Simulation details
The DPD simulations of aqueous surfactant solutions in the bulk and at hydrophobic surfaces were performed at ρ * = 3 and T * = 1, considering T 0 = 300 K. For the bulk simulations, we used a rectangular simulation box 20 r c × 20 r c × 40 r c with periodic boundary conditions in all three directions [24]. For the simulations with hydrophobic surfaces, we employed a rectangular simulation box 20 r c × 20 r c × 42 r c with confining walls of the width equal to 1 r c placed perpendicular to the z-direction and with periodic boundary conditions in the x-and y-directions only [27]. The box size 42 r c in the z-direction and wall width 1 r c guarantee bulklike behaviour in the middle of the simulation box and no mutual influence of the opposing walls. The walls were modelled as an equilibrium layer of frozen beads with the same density as in the bulk solution which also allows us to employ the interaction parameters a * ij and R * ij for modelling the solid-liquid interactions. Because the soft repulsion between DPD beads cannot prevent the solution from penetrating the walls, we separated the solution and wall domains by a reflecting surface and imposed mirror reflections when solution beads are about to cross the solid surface. The reflection scheme suppressed unphysical fluctuations in the solution density close the wall [46].
We further employed t * = 0.01, and γ * = 4.5 which may affect the estimate of time in real units based on t = t * r c / kT 0 /m 0 . Instead, one may follow Groot and Rabone [47] and gauge the physical unit of time τ by matching the experimental self-diffusivity of water, D exp w , with that of DPD model water, D DPD * w , which leads to where N m = 2 in this work, D exp w = 2.30 · 10 −9 m 2 /s [48], and D DPD * w 0.3 [49]. Equation (10) gives τ 80 ps and in turn, t 0.9 ps. The simulation started from random configurations and production runs spanned from 20 · 10 6 to 50 · 10 6 time steps, following an equilibration period. The DPD trajectories were generated using the LAMMPS simulation package [50], followed by post-processing using in-house codes to evaluate the system observables.

System Observables
The post-processing of the saved trajectories began with the identification of the surfactant associates (micelles) where we considered that two surfactants belong to the same associate if the distance between any pair of their hydrophobic alkane beads was less than the cut-off distance, r c [26]. It was followed by the evaluation of the distributions of system observables, O, such as the association number, A S . We computed the number-, weight-, and z-distributions as respectively, along with the corresponding number-, weight-, and z-averages respectively. In Equations (11) and (12) The surfactant self-assembly takes place above the critical micelle concentration (CMC), and its value for the studied systems was determined from the number distribution function of association numbers, F n (A S ). For a surfactant system with a concentration, c, above the CMC, the CMC is equal to the concentration of 'free' surfactants, c f , corresponding to the number of surfactants in submicellar clusters with A S 's below a cut-off association number, A S,m , i.e. CMC ≡ c f for A S < A S,m [51].
The gyration tensor is a convenient measure of the shape and size of surfactant associates [30,52]. In Equation (13),Ñ is the number of surfactant beads in an associate, r m i stands for the mth Cartesian coordinate of the position vector of bead i, r i , and r m mean is the m-th coordinate of the average position vector of all beads in the surfactant associate, r mean . The gyration tensor is symmetrical with eigenvalues g 2 a ≤ g 2 b ≤ g 2 c which are proportional to the principal half-axes of an equivalent ellipsoid and therefore, their ratio can be used to estimate the shape of the surfactant associate. The radius of gyration, R G , which can be obtained experimentally using static light scattering or small-angle X-ray and/or neutron scattering, is then determined via For the adsorption of surfactants at hydrophobic surfaces, we considered that a surfactant is adsorbed if at least one of its alkane-tail beads is within a distance r w from the wall which in turn determines the number of adsorbed surfactants, n adsorb . We employed r w = 1 r c and expressed the amount of the adsorbed surfactants using its adsorbed surface density where L x and L y are the box size in the x-and ydirections, respectively [46]. The orientation of an adsorbed surfactant chain was described by the orientational-order parameter where θ is the angle between a specified chain axis and the surface normal. Besides S z for the adsorbed surfactant chains, we also evaluated S z separately for the surfactant head E m and tail C n chains where the specified chain axis corresponds to a vector between the first and last beads of the chains. The values of S z vary between −0.5 and 1; values close to −0.5 and 1 indicate that the chains have become parallel and perpendicular, respectively, with respect to the surface normal. A random orientation of chains corresponds to the value of S z equal to 0.25, i.e. cos 2 θ ranging uniformly between 0 and 1 [53]. For aqueous surfactant solutions at hydrophobic surfaces, we also calculated the number density profiles, ρ α (z), for water and surfactant beads as where i δ(z − z α i ) gives the number of α beads inside a slice of the simulation box of the thickness z along the z -axis. The density profiles are plotted relative to the density of liquid water.

Results and discussion
We performed DPD simulations of aqueous surfactant solutions in the bulk and at hydrophobic surfaces, varying the surfactant concentration in a wide range. We considered nonionic surfactants C n E m with different lengths Table 2. Surfactant concentrations in terms of the mass fraction, w, and molar concentration, c, along with the number of surfactant chains, n, in simulation boxes. The size of the simulation box is 20r c × 20r c × 40r c for the bulk systems and 20r c × 20r c × 42r c for the systems with hydrophobic surfaces with the bead number density ρr 3 c = 3; r c 5.65 Å is the cut-off distance. of both the hydrophobic C n tail and hydrophilic E m head chains; specifically, C 6 E 3 , C 6 E 4 , C 8 E 3 , and C 8 E 4 . Table 2 summarises surfactant concentrations studied in terms of the mass fraction, w, and molar concentration, c, along with the corresponding number of the surfactant chains, n. The simulations of the bulk phase provide necessary information about the values of the CMC and the distribution, size, and shape of the surfactant associates while the simulations at hydrophobic surfaces focus on molecular-level insights into the surfactant adsorption and self-organisation, and their interplay with the bulk surfactant self-assembly. Figure 2(a,d) shows typical simulation snapshots for the aqueous C 6 E 3 and C 8 E 3 solutions, respectively, at a surfactant concentration above the CMC. For C 6 E 3 systems (Figure 2(a)), smaller micelles coexist with many submicellar clusters and free chains. For C 8 E 3 systems (Figure 2(d)), there are larger micelles and only a few free chains. A similar behaviour was observed for the aqueous C 6 E 4 and C 8 E 4 solutions (see Figure S1(a,d) in the Supplemental Material), indicating that a change in the hydrophobic C n tail length has more pronounced effects on the C n E m micellization behaviour than a change in the hydrophilic E m head length. Figure 3 then quantitatively elaborates on the micellization behaviour by plotting the weight distribution function of association numbers, F w (A S ), for different surfactant concentrations. The C 6 E 3 and C 6 E 4 concentrations range from fully dissolved surfactants to surfactant concentrations above the CMC, as evident, respectively, in Figure 3(a,b), by the transition from monotonously decreasing F w (A S ) to the bimodal distribution with increasing surfactant concentration. On the other hand, all C 8 E 3 and C 8 E 4 concentrations are above the CMC as evidenced by the bimodal shapes of their F w (A S ) ( Figure 3(c,d)). Similarly to Figure 2(a,d), and Figure S1(a,d) in the Supplemental Material, lengthening the hydrophilic E m head does not significantly affect the micellar behaviour, whereas lengthening the hydrophobic C n tail greatly reduces surfactant solubility and increases its association number.

Bulk phase
Determination of the CMC is not straightforward because micellization is not a true first-order phase transition. In an idealised behaviour of micellar systems, the concentration of free surfactants is equal to the surfactant concentration below the CMC while above the CMC, the concentration of free surfactants remains constant with increasing surfactant concentration because all added surfactants form micelles. The transition is therefore sharp unlike in real cases where the transition is gradual, and the concentration of free surfactants may depend on the surfactant concentration even above the CMC. We determined the values of the CMC from a submicellar portion of F n (A S ), detailed in insets of Figure 3(a,b), as the concentration of free surfactants, c f [51]. We counted the free surfactants up to the cut-off association number, A S,m , equal to 3 or 4 for the aqueous C 6 E 3 and C 6 E 4 The calculated values of the CMC for the aqueous C n E m solutions are summarised in Table 3 where they are compared with experimental data [54,55]. Note that experimental measurements of the CMC are based on abrupt changes in properties such as the surface tension or osmotic pressure, leading to CMC values depending on the experimental techniques employed. Therefore, Table 3 reports the experimental CMCs as a range of values from different experimental techniques. The calculated values of the CMC are in excellent agreement with experimental values [54,55]. The CMCs for the aqueous C 6 E 3 and C 6 E 4 solutions are almost an order of magnitude higher than those for the aqueous C 8 E 3 and C 8 E 4 solutions, whereas the CMCs for the aqueous C n E m solutions differing in the length of the hydrophilic head (C 6 E 3 vs. C 6 E 4 and C 8 E 3 vs. C 8 E 4 ) are roughly the same.
In addition, we also evaluated the eigenvalues of the gyration tensor for C n E m surfactant associates along with the average association number, A S , and average radius of gyration, R G , which are also listed in Table 3. All the three eigenvalues of the gyration tensor are roughly equal, indicating that the surfactant associates are spherical. In Table 3, we report ranges of A S and R G values because generally, the values of the observables A S and R G depend on the surfactant concentration. The ranges are for wellequilibrated aqueous C n E m solutions above the CMC; namely, c = 87−347 mM for C 6 E n and c = 58−216 mM for C 8 E n . The observables A S and R G show similar trends as the CMCs, i.e. lengthening hydrophobic C n tail affects A S and R G values more profoundly than lengthening the hydrophilic E m head. For example, the values of A S for the aqueous C 6 E 3 and C 6 E 4 solutions are about half of those for the aqueous C 8 E 3 and C 8 E 4 solutions. Although A S slightly decreases with the length of the hydrophilic E m head, R G increases since the longer E m heads spread further from the associate cores to the solvent.

Hydrophobic surfaces
Typical simulation snapshots for the aqueous C 6 E 3 and C 8 E 3 solutions at the hydrophobic surface ( Figure  2(b,c,e,f), respectively) at two surfactant concentrations illustrate that the surfactants strongly adsorb onto the surface. At low surfactant concentration (Figure 2(b,e)), the adsorbed surfactants are in equilibrium with a bulk solution away from the walls. Surfactant concentration in the bulk solution is around the CMC, and the surfactants form submicellar clusters. At high surfactant Table 3. Values of the critical micelle concentration from our simulations, CMC sim , and experiments [54,55], CMC exp , average number-, weight-, and z-association numbers, A S n , A S w , and A S z , respectively, and average number-, weight-, and z-radii of gyration, R G n , R G w , and R G z , respectively. The simulation uncertainties of CMC sim are given as subscripts. The ranges of values for A S and R G correspond to concentration ranges of c = 58−216 mM for aqueous C 8 E n solutions and c = 87−347 mM for aqueous C 6 E n solutions.
concentration (Figure 2(c,f)), the adsorbed surfactants are in equilibrium with a bulk micellar phase. A similar behaviour is observed for the C 6 E 4 and C 8 E 4 systems; see Figure S1(b,c,e,f) in the Supplemental Material. The interplay between surfactant adsorption at the hydrophobic surface and surfactant micellization in the bulk solution is quantified by plotting the fraction of the adsorbed surfactants, n adsorb /n, as a function of the number of surfactants, n, in Figure 4(a). At the lowest n studied, most of the surfactants (more than 80 % of C 6 E m surfactants and about 98 % of C 8 E m surfactants) adsorb onto the surface, forming disordered domains and leaving the concentration of surfactants in the bulk water below or around the CMC. Upon increasing n, adsorbed surfactants self-assemble into hemispheres and the surface becomes increasingly unavailable for adsorption. Simultaneously, surfactants start self-assembling into micelles away from the surface, i.e. in the bulk water. The surfactants in the bulk solution feel less attracted towards the surface due to the screening of their interactions with the surface. These lead to the decrease in n adsorb /n with increasing n. Figure 4(a) further indicates that the C n E m surfactant adsorption on the surface is enhanced by lengthening the hydrophobic C n tail and conversely, weakened by lengthening the hydrophilic E m head. Figure 4(b) then confirms that the surface does not affect micellization in the bulk solution. Here, we compare the weight average aggregation number, A S w , as a function of the surfactant bulk concentration, c bulk , from simulations of the bulk phase (see Subsection 3.1) with A S w vs. c bulk from simulations with the hydrophobic surfaces, considering only surfactants away from the surfaces. Figure 4(b) shows excellent agreement between the two types of simulations. Note that in the latter simulations, surfactants belonging to the bulk solution as well as the size of the bulk domain, which in turn defines c bulk , were determined based on the water density profile; see below. Figure 5 shows adsorption isotherms for the aqueous C n E m surfactant solutions expressed as the surfactant adsorbed surface density, ρ adsorb , against c bulk with vertical dashed lines denoting the calculated values of the CMC (see Table 3). For less readily adsorbing and more soluble C 6 E m surfactants, ρ adsorb increases relatively quickly with increasing c bulk below the CMC. Above the CMC, the increase in ρ adsorb with increasing c bulk becomes less pronounced, and ρ adsorb may eventually reach its saturation value at high c bulk . For more readily adsorbing and less soluble C 8 E m surfactants, ρ adsorb increases rapidly below the CMC. Any further increase in ρ adsorb when c bulk increases above the CMC is small, indicating that ρ adsorb is close to its saturation value around the CMC. The adsorption isotherms further elaborate on the interplay between the hydrophobic C n tail length and the hydrophilic E m head length on the enhancement/weakening of the C n E m surfactant adsorption on the hydrophobic surface, as already seen in Figure 4(a). Specifically, ρ adsorb (C 8 E 3 ) > ρ adsorb (C 6 E 3 ) and ρ adsorb (C 8 E 4 ) > ρ adsorb (C 6 E 4 ) while ρ adsorb (C 6 E 3 ) > ρ adsorb (C 6 E 4 ) and ρ adsorb (C 8 E 3 ) > ρ adsorb (C 8 E 4 ). Figures 6 and 7 present number density profiles, ρ(z), for the aqueous C 6 E 3 and C 8 E 3 surfactant solutions, respectively, in the interfacial region of the hydrophobic surface at low and high surfactant concentrations. Figures S2 and S3 in the Supplemental Material then show analogous ρ(z) for the aqueous C 6 E 4 and C 8 E 4 surfactant solutions, respectively. We evaluated ρ(z) separately for the water beads, surfactant hydrophobic C n tail, and surfactant hydrophilic E m head, ρ w (z), ρ t (z) and ρ h (z), respectively, as well as for the individual surfactant beads, ρ α (z). The density profiles are normalised by the density of pure water, ρ water , and z is the distance of a bead from the wall. In the snapshots, the spheres forming chains represent E 3 head and C 6 tail of the surfactants while packed spheres depict the walls; water beads are suppressed for the sake of clarity. The density profiles are normalised by the density of pure water, ρ water , and z is the distance of a bead from the wall. In the snapshots, the spheres forming chains represent E 3 head and C 8 tail of the surfactants while packed spheres depict the walls; water beads are suppressed for the sake of clarity.
Density profiles ρ t (z) and ρ h (z) show two distinct maxima, indicating the presence of two interfacial layers. The surfactants lying directly on the surface are characterised by strong density peaks at the solid-liquid interface, while the second layer is characterised by the minor density peaks at the distance (0.5-0.6) nm from the surface. By lying on the surface, the first-layer surfactants occupy a larger surface area [16] and maximise the interactions of their hydrophobic C n tails with the surface, but due to the chain rigidity, the hydrophilic E m heads cannot fully point away from the surface, leaving the surfactant tails exposed to water. Nevertheless, the tendency of the hydrophilic head beads to move away from the hydrophobic surface is evidenced by ρ t > ρ h and ρ occ 1 > ρ occ 2 > ρ occ 3 > ρ o at z = 0. The secondlayer surfactants are also parallel to the surface to maximise their hydrophobic-hydrophobic interactions with the first layer as indicated by the similar positions of the minor ρ t , ρ h , and ρ α peaks. However, being further from the surface gives the chains more freedom to angle their hydrophilic chains into water as suggested by the minor ρ h peak being slightly further from the wall than the minor ρ t peak. At higher surfactant concentrations, this behaviour leads to the formation of hemispheres at the surface visible in the snapshots in Figures 6 and 7 rather than the formation of an adsorbed monolayer. The above-described tendencies become more pronounced with increasing surfactant concentration.
Water number density profiles, ρ w (z), exhibit a depletion layer at the surface due to surface hydrophobicity which is mediated by the presence of the surfactant interfacial layers. Water beads can solvate the surfactant hydrophilic E m heads in the interfacial layers, leading to the presence of two water density peaks in ρ w (z) inside the depletion layer. The water depletion layer becomes more pronounced and the water density peaks inside the depletion layer lower upon an increase of the surfactant concentration. The latter is due to increasing hydrophilic-hydrophobic interactions. Finally, the adsorption behaviour of the C 6 E m surfactants in Figures 6 and 7 and that of C 8 E m in Figures S2 and S3 of the Supplemental Material differs only quantitatively but not qualitatively. Figure 8. The average of the orientational-order parameter, S z , as a function of the surfactant adsorbed surface density, ρ adsorb , for the C n E m surfactants adsorbed on the hydrophobic surface. Besides S z for the whole C n E m chains, S z was also separately evaluated for the C n tail and E m head of the adsorbed surfactants. Figure 8 elaborates on the lying-down configurations of the adsorbed C n E m surfactants by plotting the average of the orientational-order parameter, S z , as a function of ρ adsorb . We evaluated S z for the whole C n E m chains (Figure 8(a)) as well as for the surfactant C n tail chains (Figure 8(b)) and E m head chains (Figure 8(c)). The negative values of S z confirm the tendency of the adsorbed surfactants to orient themselves parallel to the hydrophobic surface, and the higher values for heads than for tails also confirm the hydrophilic heads are oriented more towards the bulk solution, while the hydrophobic tails maximise their interactions with the surface. The S z dependence on the size of surfactant blocks reflects the interplay between the length of the hydrophobic C n tail and the length of the hydrophilic E m head, already seen in Figures 4(a) and 5. The values of S z increase with increasing ρ adsorb which is associated with self-assembly of the adsorbed surfactants into the hemispheres.

Conclusion
We employed dissipative particle dynamics [41] to study the interplay between the self-organisation of aqueous surfactants solutions in the bulk, and adsorption and self-organisation of the aqueous surfactant solutions at a hydrophobic surface. We considered the nonionic surfactants, n-alkyl poly(ethylene oxide) C n E m , in water which were modelled by a generic mesoscopic force field [37]. The force field coarse-grained the surfactants and water into mesoscopic beads with 1-3 heavy atoms and two water molecules, respectively. The hydrophobic surface was represented by frozen methyl beads which mimicked carbon-type, hydrophobic surfaces at the mesoscale level. The force field relates the size of the mesoscopic beads to the molar volume of underlying molecular groups and it calibrates the bead-bead interaction parameters against the water-octanol partition coefficients.
For the aqueous C 6 E 3 , C 6 E 4 , C 8 E 3 , and C 8 E 4 solutions, we first explored their self-organisation in the bulk. Our mesoscopic simulations showed the formation of spherical micelles above the critical micelle concentration (CMC) which were further characterised by distribution functions of association numbers, average association numbers, and radii of gyration. The simulated values of the CMC were in excellent agreement with experiments [54,55], which provided evidence about the predictive power of the employed mesoscopic force field. The bulk surfactant micellization behaviour is primarily affected by the alkyl chain length which is for example demonstrated by an order of magnitude decrease in the CMC when moving from the aqueous C 6 E m to aqueous C 8 E m solutions.
The surfactants C 6 E 3 , C 6 E 4 , C 8 E 3 , and C 8 E 4 strongly adsorb on the hydrophobic surface. The surfactant adsorption behaviour was quantified by the adsorption isotherm, density profiles and chain orientational-order parameters and in contrast to the bulk phase, was influenced by both the alkyl tail length and head length. At low surfactant concentrations, the adsorbed surfactants formed disordered domains at the surface which were in equilibrium with submicellar clusters in the bulk solution. At high surfactant concentrations, the adsorbed surfactants self-organised into hemispheres which were in equilibrium with spherical micelles in the bulk. Surfactants formed two interfacial layers: a strong layer on the surface and a weak layer at (0.5-0.6) nm away from the surface; in both the layers, the surfactants tend to adopt lying down/parallel configurations. In the former layer, the surfactants maximised their interactions with the surface and occupied a larger surface area. The latter layer was caused by the hydrophobic-hydrophobic interactions between the surfactants in this layer and the adsorbed surfactants. Water was depleted from the surface as a result of surface hydrophobicity which was mediated by water molecules that solvate the surfactant hydrophilic E m heads in the interfacial layers.
Although our simulations on the surfactant adsorption were performed in the slit-pore geometry, the wall distance used was large enough to suppress confinement effects. Surfactant systems confined in nanopores are also a topic of great importance. When confined in planar slits, adsorption interactions, symmetry breaking, structural frustration, and confinement-induced entropy loss influence the morphology of the surfactant selfassembled domains, where transitions into surfactant self-assembled structures not found in the bulk are possible [30,31,34].