Nanoscale surfactant transport: bridging molecular and continuum models

Abstract

Surfactant transport is central to a diverse range of natural phenomena with numerous practical applications in physics and engineering. Surprisingly, this process remains relatively poorly understood at the molecular scale. Here, we use non-equilibrium molecular dynamics (NEMD) simulations to study the spreading of sodium dodecyl sulphate on a thin film of liquid water. The molecular form of the control volume is extended to a coordinate system moving with the liquid–vapour interface to track surfactant spreading. We use this to compare the NEMD results to the continuum description of surfactant transport on an interface. By including the molecular details in the continuum model, we establish that the transport equation preserves substantial accuracy in capturing the underlying physics. Moreover, the relative importance of the different mechanisms involved in the transport process is identified. Consequently, we derive a novel exact molecular equation for surfactant transport along a deforming surface. Close agreement between the two conceptually different approaches, i.e. NEMD simulations and the numerical solution of the continuum equation, is found as measured by the surfactant concentration profiles, and the time dependence of the so-called spreading length. The current study focuses on a relatively simple specific solvent–surfactant system, and the observed agreement with the continuum model may not arise for more complicated industrially relevant surfactants and anti-foaming agents. In such cases, the continuum approach may fail to predict accompanying phase transitions, which can still be captured through the NEMD framework.

1. Introduction

The dynamics of surfactants at the liquid–vapour interface represents a critical area of study with widespread implications across various disciplines, including engineering and medicine. This is underscored by a significant body of research highlighting its importance in processes ranging from respiration (Stetten et al. 2018; Possmayer et al. 2023) to breaking waves (Erinin et al. 2023), and to solar energy applications (Morciano et al. 2020). Surfactants, by virtue of their amphiphilic structure, play a pivotal role in modulating interfacial dynamics. These molecules align at the interface, orientating their hydrophilic heads towards, and hydrophobic tails away from, the aqueous phase. The surfactant layer adsorbed at the interface reduces the net cohesion or surface tension of the surface region. The spread of surfactants, and their subsequent reorganisation and reorientation, leads to gradients in concentration, causing differential surface tension within the surface plane. Such gradients result in the development of Marangoni stresses (Scriven & Sternling 1960). In addition to their crucial role in film thinning, which eventually gives rise to spontaneous rupture (Neél & Villermaux 2018; Ruckenstein & Jain 1974; Rahman et al. 2024), Marangoni flows are known to have influences on interfacial dynamics in complex ways; see Basu & Gianchandani (2007), Roché et al. (2014), Liu et al. (2017), Kitahata & Yoshinaga (2018), Trittel et al. (2019) and Benouaguef et al. (2021).

The migration of surfactant molecules across the liquid–vapour interface involves the formation of a precursor film that alters surface tension in its path (Afsar-Siddiqui, Luckham & Matar 2003). The extent of this migration, or spreading, follows a power law of the form $l(t) \sim kt^n$ , where $l(t)$ denotes the length of the spreading region as a function of time $t$ . The temporal variation of the exponent represents a transition between the dominant spreading mechanisms: for instance, $n\sim 1/10$ for capillarity-driven spreading, while for gravity-driven spreading $n\sim 1/8$ (Afsar-Siddiqui et al. 2003). The implicit role played by the collective behaviour of the surfactant molecules has led to them being referred to as the ‘hidden variables’ controlling fluid flows (Manikantan & Squires 2020). Empirical observations, such as the findings of Benouaguef et al. (2021) using particle image velocimetry, reveal the complexity of such flows. This underscores the limitations of analytical approaches in capturing the intricacies of surfactant-driven fluid motion.

Our understanding of surfactant transport has been advanced significantly through the contributions of several researchers, including Stone (1990), Gaver & Grotberg (1990), Crowdy (2021), Bickel & Detcheverry (2022) and Crowdy, Curran & Papageorgiou (2023). In particular, Stone (1990) derived a simple yet versatile convection–diffusion equation, building on the form of Scriven (1960), to describe the transport of surfactant molecules along a deforming liquid–vapour interface. Gaver & Grotberg (1990) showed that the initial distribution of surface tension, which depends directly on the surfactant concentration, dictates the system dynamics. Later, Crowdy (2021) showed that the complex Burgers equation can be employed to describe surfactant transport, which was further developed through exact solutions (Bickel & Detcheverry 2022; Crowdy et al. 2023). However, the case-specific nature of these solutions hinders efforts to establish a broader understanding of the underlying physical principles (Temprano-Coleto & Stone 2024).

In addition to the inherent complexities that lie in describing the transport process, the focus of this present study is to investigate the mechanism at the smallest scales where local fluctuations are crucial. Microscale and nanoscale processes inherently involve complexity and heterogeneity of molecular-level interactions, with the consequence of highly localised and dynamic surface properties evolving over space and time. Conventional continuum models simplify these details into spatio-temporal averaged and static parameters, leading to inaccuracies in predicting real-world behaviour at the nanoscale (Ilgen et al. 2024). For example, continuum models use an average value for the transport coefficients, such as the diffusion coefficient, which can be obtained independently from equilibrium molecular dynamics (MD). However on the microscopic length scale, their spatial, temporal and concentration dependence is of crucial importance which non-equilibrium MD (NEMD) implicitly includes, and notable advances on these fronts have been made (Rideg et al. 2012; Tomilov et al. 2012).

The question of whether a full three-dimensional hydrodynamic description holds at the molecular dimension has been explored for several decades (Davidovitch, Moro & Stone 2005). The continuum framework is known often to fall short in addressing the increasing importance of spatio-temporal variations and fluctuations of liquid properties on molecularly small scales (Moseler & Landman 2000). Growing success in employing MD simulations to describe dynamics below the ‘hydrodynamic’ threshold has been found. These have validated and extended continuum theories to the nanoscale (Bocquet & Charlaix 2010; Lohse & Zhang 2015; Koplik & Maldarelli 2017; Maldarelli et al. 2022). For example, the assumption of the continuum no-slip boundary condition for hydrophobic surfaces has been revisited using MD, which has demonstrated a strong link between surface atom interactions and water slippage (Huang et al. 2008). The NEMD simulations in the interfacial region have shown that the Stokes solutions fail to capture velocity profiles accurately (Bocquet & Barrat 2007). Also, MD simulations have provided confirmation of the strong dependence of electro-osmotic flow on surface wettability (Joly et al. 2004, 2006).

For the molecular picture of an interface at the smallest scale, the pioneering work of Chacón & Tarazona (2003) fitted an intrinsic interface down to the molecular spacing that was then used in later works to understand the origins of capillary forces (Delgado-Buscalioni, Chacon & Tarazona 2008), pick apart the interface layer by layer (Sega, Fábián & Jedlovszky 2015), or explore the molecular stress holding the interface together (Braga et al. 2018; Rahman et al. 2022). Furthermore, MD has been very useful in exploring systems involving chemical agents, such as surfactants or anti-foams, where continuum descriptions are insufficient. Phenomena like adsorption, micellization (Marrink, Tieleman & Mark 2000; Kanduc et al. 2024) and gas enrichment (Dammer & Lohse 2006) are inherently molecular processes, and are beyond the descriptive capabilities of conventional continuum models.

A number of studies added a stochastic stress contribution to the continuum model in order to capture the effects of thermal fluctuations on the nanoscale (Moseler & Landman 2000; Davidovitch et al. 2005; Zhang, Sprittles & Lockerby 2021; Sprittles et al. 2023). Davidovitch et al. (2005), for example, studied the spreading of viscous drops on solid substrates, and observed significantly different behaviour compared to Tanner’s law classical theory of spreading. Moseler & Landman (2000) investigated the break-up of jets of molecular diameters, where continuum assumptions cease to be valid. The authors showed that just before the break-up event, the dynamics of the system is dominated by thermal fluctuations. Since the scale in question uniquely suits the strengths of molecular simulations (Evans 1979), these works compared and validated their stochastic models against MD results. Perumanath et al. (2023) identified a novel wetting mode in their MD study of droplet spreading on a solid surface. This observation led to the subsequent development of a theoretical model. A slip boundary condition derived from MD simulations captures nanoscale dynamics typically overlooked by conventional approaches.

Quite surprisingly, the nanoscale processes that govern surfactant transport along the interface have not been comprehensively understood (Maldarelli et al. 2022), and the dynamic nature of the surfactant transport process makes NEMD a particularly appropriate tool to scrutinise these problems. The MD community has focused largely on understanding the structural and thermodynamic properties of surfactant systems (Chanda & Bandyopadhyay 2005; Rideg et al. 2012; Shi et al. 2019; Bui et al. 2021), surface tension calculations (Sresht et al. 2017; Peng et al. 2021), surface–bulk surfactant partitions (Kanduc et al. 2024; Yang & Sun 2014), micelle formation, adsorption kinetics and self-assembly (Marrink et al. 2000; Shinoda, DeVane & Klein 2008; Sangwai & Sureshkumar 2011; Wang & Larson 2015; Weiand et al. 2023; Kanduc et al. 2024). There are molecular simulations that have explored transport across the interface (Ahn et al. 2011; Xu et al. 2017), coalescence of surfactant-laden droplets (Arbabi et al. 2023, 2024), and the spreading mechanism of surfactant-laden droplets on solid substrates (Kim, Qin & Fichthorn 2006; Theodorakis, Smith & Müller 2019). Laradji & Mouritsen (2000) showed that a variation of the surfactant chain length alters the bending rigidity of the interface in a non-monotonic way. In addition to reducing surface tension, the presence of surfactant increases the thermal roughness of the surface, and thereby contributes to its instability, which is a non-intuitive outcome from a recent MD investigation by Zhang & Ding (2023). Regarding the dynamical behaviour of interfacial surfactant molecules, Johansson, Galliéro & Legendre (2022) studied a system of two immiscible Lennard-Jones fluids with surfactants at the interface. Their study highlighted the importance of incorporating local surfactant concentration and interfacial viscosity variation in modelling these flows accurately.

Despite these previous works, the literature still lacks a comprehensive molecular-level investigation of the transport processes of industrially relevant surfactants along the deforming surfaces. This study aims to fill this gap by modelling the transport of a prototypical anionic surfactant, sodium dodecyl sulphate (SDS) on the liquid–vapour surface of a thin water film by employing NEMD simulations, and simultaneously solving the transport equation. An objective of this work is to establish the extent to which the continuum transport equation can model the redistribution of surfactant molecules, and what aspects of the redistribution process require a molecular description.

In what follows, § 2 provides the necessary details of the modelled system. In § 2.1, the specifics of the coarse-grained model and the nature of the MD simulations employed in this work are presented. Section 2.2 describes the continuum framework for surfactant transport, elaborating on the macroscopic view of the transport dynamics. In §§ 2.3 and 2.4, the molecular and continuum descriptions are connected by deriving a molecular-level transport equation that provides a physical interpretation of each term relevant to the transport mechanism. Details of the numerical methods used to solve the continuum model are given in § 2.5. Section 3 presents the results and discussion. The exact balance of the proposed molecular-scale transport equation is shown in § 3.1. Particular emphasis is placed on the time-dependent evolution of the surfactant concentration (§ 3.2), and the so-called ‘spreading length’ (§ 3.3). A summary of the main conclusions arising from this work is presented in § 4.

2. System description and methodology

The primary objective of this study is to achieve a molecular-level understanding of surfactant spreading over thin liquid water films. For illustration of the transport phenomena, figure 1(a) shows macroscale experimental images of the dispersion of surfactant molecules across a thin layer of soapy water. The procedure for this thin film experiment and the visualisation techniques involved are described in an earlier study by one of the authors, Shen et al. (2020). Initially deposited over a small region on the left-hand side of the film, the surfactant layer eventually spreads to cover the entire film. As the colour maps show, regions of variable thickness emerge as the transport processes progresses. Figure 1(b) shows the MD geometry for the present study, where a thin water film of length $2L=100\,\mathrm {nm}$ , width $w=20\,\mathrm {nm}$ and thickness $h=10\,\mathrm {nm}$ is considered. A soluble surfactant monolayer composed of SDS molecules was placed in the central region of both the top and bottom surfaces of the film, which spanned area $20\times 20\, \mathrm {nm^2}$ ; see the top panel of figure 1(b). A droplet simulation (shown in figure 1(c), $R=17\,\mathrm {nm}$ , $w=20\,\mathrm {nm}$ , effective length $2L=2\pi R \approx 106.8\, \mathrm {nm}$ with similar volume to the film) is carried out with a surfactant monolayer added to a point on the surface. As is well known, the presence of the surfactant molecules reduces the local surface tension when compared with the surfactant-free regions, which is responsible for the ensuing evolution of the liquid film.

Figure 1. (a) Illustration of the surfactant spreading over a film of soapy water obtained from experiments. The colour maps show the thinning and thickening of the film. (b) The MD arrangement of a monolayer of model SDS molecules deposited on the central area of a thin water film of dimensions $100\, \mathrm {nm} \times 20\, \mathrm {nm} \times 10\, \mathrm {nm}$ , where the third dimension denotes film thickness $h_0$ . The top panel illustrates surfactant initially concentrated over area $20\, \mathrm {nm} \times 20\, \mathrm {nm}$ . The bottom panel illustrates the spreading of the monolayer across the surface, and subsequent deformation of the film. The local curvature of the surfactant–water–air interface is shown in the zoomed inset; the surface element ${\rm d}\boldsymbol{s}$ is resolved into normal $\boldsymbol{e}_n$ and tangential $\boldsymbol{e}_{t_{1}}$ and $\boldsymbol{e}_{t_{2}}$ components, or in 2D when $y$ is averaged as $\boldsymbol{e}_t$ . The solid line shows $y$ -averaged spline fitting $\zeta (x,t)$ to the surface. The mapping from $(x, y, z)$ space to $(\chi , \psi , \omega )$ space through Jacobian $J$ is illustrated schematically. (c) Surfactant transport along the surface of a (cylindrical) droplet of radius $R\sim 17\,\mathrm {nm}$ and width $w\sim 20$ nm. Initially, the surfactant molecules occupied an area $\sim 12\,\text {nm}\times w$ at both the top and bottom surfaces. In both (b) and (c), only the upper halves of the film/drop are shown.

2.1. Molecular description of the system

A coarse-grained modelling approach is employed. The MARTINI polarisable water model (Yesylevskyy et al. 2010) used here is composed of three particles: a charge-neutral central particle $W$ , accompanied by positively charged $W^+$ and negatively charged $W^-$ particles. The central particle $W$ engages in Lennard-Jones interactions with surrounding particles, whereas $W^+$ and $W^-$ interact via Coulombic forces, with no interaction occurring within the same water bead. The mass of each coarse-grained bead is set at $72$ amu, reflecting the aggregate mass of four actual water molecules. The SDS molecule is described by a five-bead coarse-grained model (Weiand et al. 2023); see figure 2. Bonded interactions were described using harmonic potentials, i.e. $V_b (r_{ij}) = 1/2 k_{ij}^b (r_{ij} - b_{ij})^2$ , where, $V_b (r_{ij})$ is the potential energy associated with the bond between atoms $i$ and $j$ at distance $r_{ij}$ , $b_{ij}$ is the equilibrium bond length between atoms $i$ and $j$ (i.e. the natural length of the bond at which the potential energy is minimised), and $k_{ij}^b$ is the force constant that determines the stiffness of a bond.

To model the energy associated with the bending of angles between three consecutive bonded atoms, the cosine-squared potential was used. This is generally expressed as $V_a (\theta ) = {1}/{2} k^\theta ( \mbox {cos} \ \theta - \mbox {cos} \ \theta _0 )^2$ , where $V_a$ is the potential energy as a function of angle $\theta$ formed by three bonded atoms $i,j,k$ , $\theta _0$ is the equilibrium angle, and $k^\theta$ is the force constant in terms of the angle. Non-bonded interactions were accounted for through a combination of Lennard-Jones and Coulombic potential terms, employing a hybrid overlay approach to capture the interplay of forces on different length scales. This is generally expressed as $V_{\textit {non-bonded}} (r) = V_{ {LJ}} (r) + V_{ {Coul.}} (r)$ , where $ V_{{LJ}} (r)$ is the Lennard-Jones potential, i.e. $V_{{LJ}} (r) = 4\epsilon [ ({\sigma }/{r} )^{12} - ({\sigma }/{r} )^6 ]$ . The Coulombic interactions are given by $V_{{Coul.}} (r) = q_i q_j/ (4\pi \epsilon _0 \epsilon _r {r})$ , where $q_i$ is the charge of index $i$ , $\epsilon _0$ is the permittivity of free space, and $\epsilon _r$ is the relative permittivity.

The initial configurations of film surfactant and drop surfactant were prepared using the open source codes Packmol (Martínez et al. 2009) and Moltemplate (Jewett et al. 2021). A velocity-Verlet integrator was used, with time step 5 fs. The long-range electrostatic forces were calculated using the particle–particle particle-mesh (PPPM) method, which enabled the electrostatic component of the system interactions to be computed at reasonable cost. The SHAKE algorithm was used to maintain rigid O–H bonds within water molecules. The initial states were energy-minimised separately, and equilibrated (figures 2a–c) in the NVT (canonical) ensemble to equilibrate the system at the target temperature 300 K. The separate equilibration was to restrict any surfactant transport before the system had time to stabilise before the production phase was started. Following this, the system underwent an NVT production phase during which data were collected for analysis. Figures 2(d,e) show, respectively, the coarse-grain descriptions, as well as the chemical structures, of water and SDS. All simulations for this study were carried out using the large-scale atomic/molecular massively parallel simulator (LAMMPS) package (Thompson et al. 2022). All properties from the MD simulations were averaged over the width of the film ( $y$ -axis) to improve the statistics.

Figure 2. (a,b) Initial configuration of an SDS layer and water film, with hydrophilic head group of SDS close to the water surface, and hydrophobic tail groups pointing away from the surface. These were equilibrated separately: equilibration of (a) SDS, (b) water. (c) Zoomed-in view of a portion of the surfactant–water interface after equilibration. Schematic of the chemical structure and the coarse-grained (CG) descriptions of (d) the MARTINI polarisable water model consisting of three beads, and (e) the SDS molecule.

Figure 3. Schematic diagram showing the mapping of a deforming surface from $(x, y, z)$ Cartesian space to $(\chi , \psi , \omega )$ space.

2.2. Continuum description of surfactant transport

The spreading of the surfactant across the liquid–air interface resulted in a locally dilute surfactant–water layer that extended over a larger surface area. This behaviour of the surfactant monolayer is described by the convection–diffusion equation (Gaver & Grotberg 1990; Stone 1990), which is a continuum model that captures the separate underlying physical mechanisms of the transport phenomena on the surface:

(2.1) \begin{equation} \underbrace {\partial _t \Gamma }_{\mathrm {unsteady}} + \underbrace { \boldsymbol{\nabla} _s \boldsymbol{\cdot} (\Gamma \boldsymbol{u}_s) }_{\mathrm {advection}} + \underbrace {\Gamma (\boldsymbol{\nabla} _s \boldsymbol{\cdot} \boldsymbol{e}_n)(\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{e}_n)}_{\mathrm {geometric\,effect}} = \underbrace {D_s\, \boldsymbol{\nabla} _s^2 \Gamma }_{\mathrm {diffusion}} \,\, +\underbrace {j_n}_{\mathrm {source\, term}}. \end{equation}

Here, $\Gamma$ denotes the local surfactant concentration, $\boldsymbol{u}_s$ represents the velocity vector along the surface, $\boldsymbol{e}_n$ is the surface normal unit vector, $\boldsymbol{\nabla} _s = (\boldsymbol{I-e}_n \boldsymbol{e}_n)$ is the surface gradient operator, and $D_s$ characterises the surfactant diffusion rate across the interface. The nature of the unsteady time derivative in (2.1) had been ambiguous, which led to potential misinterpretations that several later studies clarified (Wong, Rumschitzki & Maldarelli 1996; Cermelli, Fried & Gurtin 2005; Pereira et al. 2007). Among these studies, Wong et al. (1996) addressed this ambiguity by deriving the surfactant mass balance equation. They showed that the time derivative must follow the interface in the surface-normal direction, which ensures that the time evolution of the surfactant concentration is captured correctly as the interface deforms. Our procedure aligns with that in Wong et al. (1996), which will be clarified in § 2.3. The second term in (2.1) describes surfactant advection along the surface. The third term, which is the product of local curvature ( $\kappa =-\boldsymbol{\nabla} _s \boldsymbol{\cdot} \boldsymbol{e}_n$ ) and the surface-normal velocity ( $v=\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{e}_n$ ), represents the geometric effect of the locally deforming interface on the transport process. Together, the second and third terms capture the effects of surface compression and expansion due to the non-uniformity of advective flows and curvature, accounting for the redistribution of the surfactant molecules driven by the various dynamical processes, including the Marangoni effect. The subscript $s$ indicates that the quantities are considered along the tangential direction to the surface. The $D_s\, \boldsymbol{\nabla} _s^2 \Gamma$ term builds in the diffusive component of the surfactant transport. The local adsorption or desorption of surfactants is denoted by $j_n$ , which acts as a source or sink in the system. Through this term, surfactant molecules are introduced to or eliminated from the interface, influencing the overall dynamics of the transport process.

2.3. Derivation of a molecular-scale transport equation

Any attempt to trace the surfactant particles in experimental measurements inevitably alters their intrinsic properties, thereby affecting the dynamics of the system (Manikantan & Squires 2020). In contrast, the molecular system under investigation in this study provides direct access to individual molecules and their instantaneous properties in a non-intrusive way, providing an ideal platform to verify the efficacy of (2.1) in describing the transport phenomena.

The link between the continuum description and a discrete molecular paradigm is formalised in the seminal work of Irving & Kirkwood (1950). This starts from the postulate that the continuum concept of density $\rho$ at a point in a molecular system can be expressed by

(2.2) \begin{align} \rho (\boldsymbol {r}, t) = \sum _{i=1}^{N} \langle m_i\, \delta (\boldsymbol {r}- \boldsymbol {r}_i) ; f \rangle, \end{align}

where the sum is over all $N$ molecules in the space, and the angular brackets $\langle \cdots ;f \rangle$ denote an inner product with the probability distribution function $f$ over the $6N$ -dimensional phase space (3 position and 3 momentum coordinates). Formally, $f$ is the Gibbs ensemble, a probability drawn from a representative ensemble of similar systems, although in practice this inner product is replaced by a time average when this can be justified by the ergodic hypothesis. Key to this approach is the Dirac delta function $\delta (\boldsymbol{\cdot} )$ , which is the manifestation of the continuum approximation in a discrete system, an infinitely high and thin peak that is non-zero only when a molecule position ${\boldsymbol {r}}_i= [x_i, y_i, z_i]$ is located at $\boldsymbol {r} = [x, y, z]$ . When a particle is at point $\boldsymbol {r}$ , its mass $m_i$ is counted towards the density at that point, with the Dirac delta function having units of inverse volume. By analogy with Irving & Kirkwood (1950), we define the density of surfactant as

(2.3) \begin{align} \rho _s(\boldsymbol {r},t) = \sum _{i=1}^{N_S} \langle m_i\, \delta (\boldsymbol {r}- \boldsymbol {r}_i) ;f \rangle . \end{align}

Here, the sum is limited to only surfactant molecules, denoted by the sum from $i=1$ to $N_S$ , which is shorthand for the more formal notation $\sum _{i \in \mathcal {S}}$ summing over the set of surfactant molecules $\mathcal {S}$ with $|\mathcal {S}| = N_S$ . This set of surfactant molecules does not change during the simulation as we simulate an NVE or NVT ensemble.

In a molecular system, all quantities must also be averaged over a region in space, so the control volume is the most natural form to work in. The integral over a three-dimensional control volume $V_s$ , which follows an interface $z=\zeta (x,y,t)$ with the corners $x^{\pm }$ , $y^{\pm }$ and $z^{\pm }$ being functions of $\zeta$ , is then

(2.4) \begin{align} \int _{V_s} \rho _s(x,y,z,t) \,{\rm d}V & = \int _{V_s} \sum _{i=1}^{N_S} \langle m_i\, \delta (\boldsymbol {r}- \boldsymbol {r}_i) ; f \rangle\, {\rm d}V \nonumber \\ & = \int _{x^-}^{x^+} \int _{y^-}^{y^+} \int _{z^-}^{z^+} \sum _{i=1}^{N_S} \langle m_i\, \delta (x - x_i)\, \delta (y - y_i)\, \delta (z - z_i) ;f \rangle \, {\rm d}z \, {\rm d}y \, {\rm d}x \nonumber \\ & = \int _{\chi ^-}^{\chi ^+} \int _{\psi ^-}^{\psi ^+} \int _{\omega ^-}^{\omega ^+} \sum _{i=1}^{N_S} \langle m_i\, \delta (\chi - \chi _i)\, \delta (\psi - \psi _i)\, \delta (\omega - \omega _i)\, J ;f \rangle \, {\rm d}\chi \, \nonumber \\ & \quad\times {\rm d}\psi \, {\rm d}\omega \nonumber \\ & = \sum _{i=1}^{N_S} \langle m_i \vartheta _i J_i ;f \rangle, \end{align}

where we have used the property that the Dirac delta function at any point in vector space $\delta (\boldsymbol {r} - \boldsymbol {r}_i)$ is the product of the three orthogonal vectors $\boldsymbol {r} = x \boldsymbol{e}_x + y \boldsymbol{e}_y + z \boldsymbol{e}_z$ , using the three Cartesian coordinates, or equivalently it could be written in terms of coordinates aligned with the two tangential directions to surface $\zeta$ at every point, denoted by $\chi$ and $\psi$ , and the surface normal $\omega$ , so $\boldsymbol {r} = \chi \boldsymbol{e}_{\chi } + \psi \boldsymbol{e}_{\psi } + \omega \boldsymbol{e_{\omega }}$ . The integral over a curved shape that deforms with the interface in real space is mapped, using a Jacobian $J=J(\chi , \psi , \omega ,t)$ , to a constant cuboidal control volume of length $\Delta \chi$ , width $\Delta \psi$ and height $\Delta \omega$ . This is centred on the interface $\zeta$ with half its total volume $\Delta \omega /2$ either side, and constant length/width values $ \Delta \chi$ and $\Delta \psi$ , respectively. The notation for the cuboid integral limits is $\chi ^{\pm } = \chi \pm \Delta \chi /2$ , and similarly for $\psi$ and $\omega$ . The sifting property of the Dirac delta function $\int \delta (x-a)\, f(x) = f(a)$ then results in a Jacobian $J_i = J(\chi _i, \psi _i, \omega _i,t)$ evaluated at the location of each molecule in the sum. The integral between constant limits is called the control volume function $\vartheta _i$ , i.e. the product of three boxcar functions with $\vartheta _i \equiv \Lambda _{\chi } \Lambda _{\psi } \Lambda _{\omega }$ , each the difference of two Heaviside functions $H$ , $\Lambda _{\chi } \equiv (H(\chi ^+ - \chi _i)-H(\chi ^- - \chi _i))$ , and similarly for $\psi$ and $\omega$ , so the product is a cuboid in three-dimensional $(\chi ,\psi ,\omega )$ space. This has the property that the derivative of $\vartheta _i$ in a given direction is a function that selects molecules crossing a square face of the cuboid in $[\chi ,\psi ,\omega ]$ space, i.e. $\partial \vartheta _i / \partial \chi = [\delta (\chi ^+ - \chi _i)-\delta (\chi ^- - \chi _i) ]\Lambda _{\psi } \Lambda _{\omega } \equiv {\rm d}S_{\chi i}$ , which is the mass flux over the $\chi$ face. We introduce the notation ${\rm d}S_{\chi i}$ , the flux of molecules in the $\chi _i$ direction crossing surfaces at $\chi ^{\pm }$ , and again similarly for $\psi$ and $\omega$ .

To express the conservation of surfactants in terms of these fluxes, the control volume function can be given the same treatment as a cuboidal control volume (Smith et al. 2012), linking time evolution in a volume to molecular surface flux terms. The Jacobian has all the time dependence of the moving volume $J_i =J(\chi _i, \psi _i, \omega _i, t)$ , so only the molecular positions are functions of time, i.e. ${\rm d} H [\chi - \chi _i(t) ]/{\rm d}t = -\dot { \chi }_i\, \delta [\chi - \chi _i(t) ]$ , where $\chi$ has no time dependence (similarly for $\psi$ and $\omega$ ). In control volume form, these equations are valid without ensemble averages (Smith et al. 2012), simply representing an instantaneous density in a volume in space, so we no longer need to take the ensemble averages $\langle \cdots ; f \rangle$ . The time evolution of a control volume on the interface is then

(2.5) \begin{align} \frac {{\rm d}}{{\rm d}t}\int _V \rho _s(x,y,z,t) \,{\rm d}V &= \frac {{\rm d}}{{\rm d}t}\sum _{i=1}^{N_S} m_i J_i \vartheta _i \nonumber \\ & =\sum _{i=1}^{N_S}m_i \left [ J_i\dot {\Lambda }_{\chi } \Lambda _{\psi } \Lambda _{\omega } + J_i\Lambda _{\chi } \dot {\Lambda }_{\psi } \Lambda _{\omega } + J_i\Lambda _{\chi } \Lambda _{\psi } \dot {\Lambda }_{\omega }\right. \nonumber \\ & \quad\left. + \dot {J}_i\Lambda _{\chi } \Lambda _{\psi } \Lambda _{\omega } \right ] \nonumber \\ & = \sum _{i=1}^{N_S}m_i \left [ - J_i \dot {\chi }_i\, {\rm d}S_{\chi i} - J_i \dot {\psi }_i\, {\rm d}S_{\psi i} - J_i \dot {\omega }_i\, {\rm d}S_{\omega i} + \dot{J_i} \vartheta _i \right ]. \end{align}

Here, the flux terms ${\rm d}S_{\chi i}$ in mapped space checks if the mapped molecule position $\chi _i$ crosses a flat surface $\chi ^+$ or $\chi ^-$ , which in real space represents flow along the tangent to the surface $\xi$ . When transformed back into real space, this indicates that there is a molecule crossing the control volume face as it moves with the interface, for $\chi$ and $\psi$ along the interface, and for $\omega$ a flux over the interface.

Rearranging the equation, the equivalent of each term in the continuum equation (2.1) is annotated as follows:

(2.6) \begin{equation} \sum _{i=1}^{N_S} \left [ \underbrace { \frac {{\rm d}}{{\rm d}t} \left ( m_i\vartheta _i J_i \right )}_{\mathrm {unsteady}} + \underbrace { m_i J_i \dot {\boldsymbol {s}}_i \boldsymbol{\cdot} {\rm d}\boldsymbol{S}_{i} \vphantom {\frac {1}{J_i}}}_{\mathrm {advection}} - \underbrace {m_i\vartheta _i J_i \frac {\dot {J_i}}{J_i}}_{\mathrm {geometric\; effect}} \right] = \sum _{i=1}^{N_S} \underbrace {m_i J_i \dot {w_i}\, {\rm d}S_{wi}^-\vphantom {\frac {1}{J_i}}}_{\mathrm {source\, term}}, \end{equation}

where the velocities along the surface are denoted by $\dot {\boldsymbol {s}}_i = [\dot {\chi }_i, \dot {\psi }_i]$ with ${\rm d}\textbf {S}_{i} = [{\rm d}S_{\chi i}, {\rm d}S_{\psi i}]$ . Equation (2.6) is tested in § 3.1 and shown to be exactly satisfied. This is not surprising since this simply counts the surfactant molecules as they move between boxes, with a Jacobian weight due to the surface deformation. More important is how these terms correspond to the continuum terms in (2.1). The justification for the equivalence of the various terms to the continuum equation denoted by the underbraces will be shown in § 2.4, where the limit of zero volume will be taken to give the differential form of the equations.

2.4. Zero-volume limit

The control volume form is the most natural for MD, consistent with binning or chunking (Thompson et al. 2022) where a grid of cells is overlayed on the MD system to obtain fields such as density, momentum and energy. However, the continuum equations such as (2.1) are recovered in the limit that cell volumes tend to zero. In this limit, each side length (e.g. $\Delta \chi$ ) tends to zero, and the limits approach each other so the surface flux term $ \lim _{\Delta \chi \to 0} {\rm d}S_{\chi i} = ({\partial }/{\partial \chi })\, \delta (\chi - \chi _i)\, \Lambda _{\omega }$ becomes the derivative of a delta function (Smith et al. 2012). The similar limit of a boxcar function is simply the delta function, $\lim _{\Delta \psi \to 0}\Lambda _{\psi } = \delta (\psi - \psi _i)$ . As we are in this limiting case of zero volume, we reintroduce the ensemble averaging $\langle \cdots ;f \rangle$ to ensure equivalence of continuum and molecular terms. Starting with the first term in (2.6), we have

(2.7) \begin{align} \lim _{\substack {\Delta \chi \to 0 \\ \Delta \psi \to 0 \\ \Delta \omega \to 0}} \sum _{i=1}^{N_S} \left \langle \frac {{\rm d}}{{\rm d}t} (m_i \vartheta _i J_i) ;f \right \rangle & = \sum _{i=1}^{N_S} \left\langle \frac {{\rm d}}{{\rm d}t} m_i \left (\delta (\chi - \chi _i)\, \delta (\psi - \psi _i)\,\delta (\omega - \omega _i)\,J_i \right ) ;f \right\rangle \nonumber \\ & = \frac {{\rm d}}{{\rm d}t} \sum _{i=1}^{N_S} \bigg \langle m_i\, \delta (\tilde {\boldsymbol {r}} - \tilde {\boldsymbol {r}}_i)\, J_i ;f \bigg \rangle = \partial _t \Gamma. \end{align}

Here, $\tilde {\boldsymbol {r}} = [\chi , \psi , \omega ]$ , with $\chi$ , $\psi$ and $\omega$ the mapped coordinates following, respectively, the two tangential and normal directions to the interface $\zeta (x,y,t)$ . This expression is non-zero only on the surface as the Delta functions ensure that we select only molecules with position $\omega _i = \omega$ . This can therefore be used as the definition of the surface concentration $\Gamma$ on the interface,

(2.8) \begin{align} \Gamma (\boldsymbol {s},t) = \sum _{i=1}^{N_S} \bigg \langle m_i\, \delta (\tilde {\boldsymbol {r}} - \tilde {\boldsymbol {r}}_i)\, J_i ;f \bigg \rangle , \end{align}

which is consistent with Wong et al. (1996) in that $\tilde {\boldsymbol {r}}$ is fixed to a coordinate following the surface, and the surface deformation is included through the $J_i$ term.

The second term of (2.6) can be seen to be the advection term in the continuum equations, as the derivative in the tangential direction of the velocity times mass of surfactant along the tangent:

(2.9) \begin{align} \lim _{\substack {\Delta \chi \to 0 \\ \Delta \psi \to 0 \\ \Delta \omega \to 0}}\sum _{i=1}^{N_S} \bigg \langle m_i J_i \left [\dot {\chi }_i\, {\rm d}S_{\chi i} + \dot {\psi }_i\, {\rm d}S_{\psi i} \right ] ;f \bigg \rangle = \frac {\partial }{\partial \boldsymbol {s}}\boldsymbol{\cdot} \sum _{i=1}^{N_S} \bigg \langle m_i J_i \dot {\boldsymbol {s}}_i\, \delta (\tilde {\boldsymbol {r}} - \tilde {\boldsymbol {r}}_i) ;f \bigg \rangle = \boldsymbol {\boldsymbol{\nabla} }_s \boldsymbol{\cdot} {\Gamma \boldsymbol {u}_s}, \end{align}

where $\boldsymbol {s} = [\chi , \psi ]$ , and velocities along the surface are $\dot {\boldsymbol {s}}_i = [\dot {\chi }_i, \dot {\psi }_i]$ . Equation (2.9) is analogous to one in Irving & Kirkwood (1950), which defined the gradient of momentum as $\boldsymbol{\nabla} \boldsymbol{\cdot} \rho \boldsymbol {u} = ({{\rm d}}/{{\rm d}\boldsymbol {r}}) \sum _{i=1}^{N} \langle m_i \dot {\boldsymbol {r}}_i\, \delta (\boldsymbol {r} - \boldsymbol {r}_i) ;f \rangle$ for a fixed Eulerian reference frame, which is extended here to a surface tracking Lagrangian one.

The third term in (2.6), the deformation of the surface in time, denoted as geometric effect in (2.6), is given by the time derivative of the Jacobian $\dot {J}$ ,

(2.10) \begin{align} \lim _{\substack {\Delta \chi \to 0 \\ \Delta \psi \to 0 \\ \Delta \omega \to 0}}\sum _{i=1}^{N_S} \left \langle m_i J_i \vartheta _i \frac {\dot {J_i}}{J_i} ;f \right\rangle =\sum _{i=1}^{N_S} \left\langle m_i J_i\, \delta (\tilde {\boldsymbol {r}} - \tilde {\boldsymbol {r}}_i)\, \frac {\dot {J_i}}{J_i} ;f \right\rangle = \Gamma \frac { \dot {a}}{2a}, \end{align}

which is simply the definition of $\Gamma$ from (2.8) multiplied by the time-evolving Jacobian divided by its original value $\dot {J}_i/J_i$ at the position of each molecule. Interestingly, this shows the effect of surface geometry in a molecular system, and is therefore weighted by the deformation at the location of each molecule.

Finally, for the right-hand side of (2.6), $\dot {\omega }_i\, {\rm d}S_{\omega i}$ , the limit of zero volume would give a derivative in the surface-normal direction:

(2.11) \begin{align} \lim _{\substack {\Delta \chi \to 0 \\ \Delta \psi \to 0 \\ \Delta \omega \to 0}}\sum _{i=1}^{N_S} \bigg \langle m_i J_i \dot {\omega }_i\, {\rm d}S_{\omega i} ;f \bigg \rangle = \frac {\partial }{\partial \omega }\sum _{i=1}^{N_S} \bigg \langle m_i J_i \dot {\omega }_i\, \delta (\tilde {\boldsymbol {r}} - \tilde {\boldsymbol {r}}_i) ;f \bigg \rangle = j_n. \end{align}

To understand the link to the surfactant source/sink term $j_n$ , it is clearer to interpret this in terms of fluxes. The surfactant flux over $\omega ^+$ given by $m_i \dot {\omega }_i\, {\rm d}S_{\omega i}^+$ can be assumed to be zero, as this surface is above the interface so a non-zero value would mean evaporation of surfactant. The flow of surfactant over the bottom surface, $m_i \dot {\omega }_i\, {\rm d}S_{\omega i}^-$ , is therefore the surfactant movement from bulk to the surface, a source/sink term in the surface equations.

Collecting all these terms, we obtain the equivalent form to the continuum differential equation:

(2.12) \begin{align} \overbrace {\sum _{i=1}^{N_S} \left \langle \frac {{\rm d}}{{\rm d}t} \bigg [ m_i J_i\, \delta (\tilde {\boldsymbol {r}} - \tilde {\boldsymbol {r}}_i) \bigg ] ;f \right \rangle }^{\mathrm {unsteady}} + \overbrace { \frac {\partial }{\partial \boldsymbol {s}} \boldsymbol{\cdot} \sum _{i=1}^{N_S} \bigg \langle m_i J_i \dot {\boldsymbol {s}}_i\, \delta (\tilde {\boldsymbol {r}} - \tilde {\boldsymbol {r}}_i) \vphantom {\frac {1}{J_i}} ;f \bigg \rangle }^{\mathrm {advection}} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \nonumber \\ {}- \underbrace { \sum _{i=1}^{N_S} \bigg \langle m_i \dot {J_i}\, \delta (\tilde {\boldsymbol {r}} - \tilde {\boldsymbol {r}}_i) ;f \bigg \rangle }_{\mathrm {geometric\, effect}} = \underbrace { \frac {\partial }{\partial \omega }\sum _{i=1}^{N_S} \bigg \langle m_i J_i \dot {\omega }_i\, \delta (\tilde {\boldsymbol {r}} - \tilde {\boldsymbol {r}}_i) ;f \bigg \rangle \vphantom {\sum _{i=1}^{N_S}}}_{\mathrm {source\, term}}. \end{align}

The over/underbraces identify the equivalence between the terms in (2.12) and those in the continuum equation (2.1). The unsteady term represents the temporal change of molecules at a point. The advection term corresponds to the flux of molecules along the tangential direction of the deforming surface. The geometric effect term, characterised by $\dot {J}_i/J_i$ , accounts for the impact of the deforming volume over time, but expressed as a weighted sum based on the location of the molecules. The source term in (2.6) describes the flux of surfactant from the bulk to the surface, under the assumption that no surfactant molecules evaporate from the surface.

Unlike the form introduced by Stone (1990), there is no obvious diffusion contribution obtained in the molecular system. This is consistent with the work of Scriven (1960), who presented a surface transport equation without diffusion. It might be possible to justify a diffusion-like effect through (i) consideration of Fick’s law style flux, (ii) interactions with the non-surfactant molecules, or (iii) fluctuations, which remain when the molecular equations are averaged over time; formally, $D_s\, \boldsymbol{\nabla} _s^2 \Gamma = \boldsymbol{\nabla} _s \boldsymbol{\cdot} \sum _{i=1}^{N_S} \langle m_i J_i ( \dot{\boldsymbol {s}}_i - \boldsymbol {u}_s)\, \delta (\tilde {\boldsymbol {r}} - \tilde {\boldsymbol {r}_i}) ; f\rangle$ . The presented equation (2.6) is formally exact, and simply counts the number of molecules moving into and out of a surface tracking (Lagrangian) volume, with a term $\dot {J}_i/J_i$ for the deformation of the volume in time. It is worth noting that the absence of an explicit diffusion term in (2.6) is analogous to kinetic theory, where balance equations derived from the Boltzmann equation lack diffusivity until small perturbations around a local equilibrium are introduced (Chapman & Cowling 1990).

Similarly, an MD framework provides an exact account of the molecular trajectories and does not explicitly include diffusion as an additional term. Diffusional behaviour becomes apparent only through analysis of these trajectories . Small-scale fluctuations have been proven to play a critical role in emergent macroscopic behaviour (Pototsky & Suslov 2024; Yang et al. 2024). Averaging over small-scale fluctuations near local equilibrium could also provide a means to recover diffusive-like contributions.

2.5. Numerical methods

In this subsection, it is shown how the control volume equations can be collected directly in an MD simulation to validate the equations in the previous section. Then the method to solve the continuum equation (2.1) using a simple finite differences scheme is discussed, with boundary, initial and other coupled values obtained from the MD.

2.5.1. Control volume fluxes

The control volume form can be used directly in a molecular simulation, having the advantage of exact conservation because it simply counts how many molecules are in each box and how many cross surfaces, with the Jacobian term scaling by the volume change. Taking the time integral of both sides of (2.6) over one time step of a molecular simulation $\Delta t$ , a form of equation is obtained that is usable directly in MD:

(2.13) \begin{align} \sum _{i=1}^{N_S} \frac {m_i\,\vartheta _i (t+\Delta t)\, J_i (t+\Delta t) - m_i\,\vartheta _i (t)\, J_i (t) }{\Delta t} + \sum _{i=1}^{N_S} \sum _{k}^{N_{{crossing}}} m_i J_i \dot {\boldsymbol {s}}_i \boldsymbol{\cdot} {\rm d}\textbf {S}_{i k } \;\;\;\;\;\;\; \;\;\;\;\;\;\; \nonumber \\ - \sum _{i=1}^{N_S} m_i\,\vartheta _i (t)\, \frac {J_i(t+\Delta t)- J_i(t)}{\Delta t} = \sum _{i=1}^{N_S} \sum _{k}^{N_{{crossing}}} m_i J_i \dot {w_i}\, {\rm d}S_{wi k}^-, \end{align}

where the $k$ subscript on ${\rm d}S_{\alpha i k }$ denotes all crossings in the time step $\Delta t$ , formally obtained by expressing the delta function of molecular position in terms of the sum of its roots in time $\delta (r - r_i(t)) = \sum _k \delta (t - t_k)/|\dot {r}_i|$ , and taking the time integral.

This complex-looking equation (2.13) has a simple algorithmic interpretation. (i) Sum the surfactant molecules in a volume at time step $t$ and at $t+\Delta t$ . The difference gives the first term, the unsteady change in surfactant concentration. (ii) Count any surfactant molecules crossing the left or right surfaces of the box; these give the second term for the advection along the interface. (iii) Calculate the Jacobian at times $t$ and $t+\Delta t$ , which can be done assuming a simple linear finite element in the volume (see SM), and then evaluating $J$ at the location of all molecules in the volume. The sum of the time-evolving Jacobian, the geometric effect, is the third term of (2.13). (iv) Finally, any molecules crossing the bottom or top surface give the fourth term for evaporation/absorption, typically zero below the Critical Micelle Concentration (CMC) for a volume wide enough to encompass the interface.

Equation (2.13) is tested directly in molecular simulations, and is shown to give an exact balance in § 3.1. It is demonstrated that (2.13) includes all relevant terms, including the geometric effect due to surface deformation, and the source term due to surfactant absorption/de-absorption into/from the bulk. Note that the surface is obtained by a fitting process to the molecular positions or average density field, so it is decoupled from the time evolution of the molecules themselves. The evolution of the molecular system is not changed by this fitting, and the surface is purely defined. Therefore, the flux along this surface can be obtained, with fluxes over the next time step obtained on a surface assumed to be constant for that time in (2.13). The surface evolution itself is calculated from the simple forward Euler approximation to $\dot {J}_i$ . In this way, it corresponds to the assumption of a system with fixed surface coordinates (Wong et al. 1996), and so is consistent with Scriven (1960).

2.5.2. Volume average forms

The surface fluxes in (2.13) are exact, directly counting molecular crossings using a mapped version of the velocity method of planes (Todd, Evans & Daivis 1995). However, the statistics are known to be poor for surface flux measurements (Smith, Heyes & Dini 2017), especially when compared to averages of quantities inside a volume, the so-called volume average forms. Since surface fluxes are equivalent to the volume measurements in the limit that the volumes are small (Heyes et al. 2011), it is often preferable to work with these. Here, the volume $\tilde {V}$ is a cuboidal control volume in mapped $(\chi , \psi , \omega)$ space:

(2.14) \begin{align} \int _{\tilde {V}} \Gamma (\boldsymbol {s},t)\, {\rm d}\tilde {V} = \sum _{i=1}^{N_S} m_i \vartheta _i J_i. \end{align}

The volume average forms make the assumption that there is a single constant continuum value in this control volume, i.e. $\int _{\tilde {V}} \Gamma\, {\rm d}\tilde {V} \approx \Gamma\, \Delta \tilde {V}$ , so that the MD surfactant concentration at location $\boldsymbol {s}$ and time step $t$ is

(2.15) \begin{equation} \overset {MD}{\Gamma }(\boldsymbol {s},t) = \frac {1}{\Delta \tilde {V}} \sum _{i=1}^{N_S} m_i J_i \vartheta _i[\boldsymbol {s}, t]. \end{equation}

Molecules are assigned to bins of width $\Delta s = s^+ - s^-$ centred on $s$ by the control volume (CV) function $\vartheta _i [s,t] = [H(s^+ - s_i(t)) - H(s^- - s_i(t))]\Lambda _{\omega }$ , giving the average surfactant concentration in that volume, assuming an average over the $y$ or $\psi$ direction. Similarly, the momentum along the surface can be obtained from the volume integral form of (2.9),

(2.16) \begin{equation} \overset {MD}{\Gamma \boldsymbol {u}_s} = \frac {1}{\Delta \tilde {V}} \sum _{i=1}^{N_S} m_i J_i \dot{\boldsymbol {s}}_i \vartheta _i[\boldsymbol {s}, t], \end{equation}

so the volume averaged velocity is (2.16) divided by (2.15):

(2.17) \begin{equation} \overset {MD}{\boldsymbol {u}_s} (\boldsymbol {s},t) = \frac {1}{\Delta \tilde {V} \overset {MD}{\Gamma } } \sum _{i=1}^{N_S} m_i J_i \dot{\boldsymbol {s}}_i \vartheta _i[\boldsymbol {s}, t]. \end{equation}

These equations are used to give boundary and initial conditions for the surfactant concentration equations, as well as to define the surface velocity.

2.5.3. Finite difference continuum solver

To evaluate the applicability of the continuum model in describing transport phenomena at the nanoscale, (2.1) is solved numerically using direct inputs from the MD simulations. Specifically, the concentration data at the left-hand boundary of the domain are taken directly from MD and used as input for the numerical model. Additionally, the velocity field (decomposed into its normal and tangential components) and the surface deformations are provided by the MD data. This formulation results in a closed system of equations, which we solve to assess how closely the convection–diffusion model (2.1) agrees with the MD system, without requiring a continuum equation to approximate the surface evolution and momentum balance.

The evolution of $\Gamma$ in time is discretised using a simple backwards Euler scheme,

(2.18) \begin{equation} \frac {\partial \Gamma }{\partial t} \approx \frac {\Gamma _{I}^{n} - \Gamma _{I}^{n-1}}{\Delta t}, \end{equation}

where superscripts denote time steps, and subscripts denote spatial cells. The advection term is one-dimensional in $s$ , hence $\boldsymbol{\nabla} _s \boldsymbol{\cdot} (\Gamma \boldsymbol{u}_s) = \partial \Gamma u / \partial s$ with the subscript $s$ dropped for notational simplicity:

(2.19) \begin{equation} \frac {\partial \Gamma u }{\partial s}= u \frac {\partial \Gamma }{\partial s} + \Gamma \frac {\partial u }{\partial s} \approx u_I \frac {\Gamma _{I+1} - \Gamma _{I} }{\Delta s} + \Gamma _I \frac {u_{I+1} - u_{I}}{\Delta s}. \end{equation}

Forward difference is used to include upwinding (Hirsch 2007) as $u\gt 0$ with surfactant moving from left to right, and $\Gamma \geqslant 0$ . The sum of (2.19) over all $N_{ {cells}}$ gives the total advection $\eta _{ {adv.}}$ at each time. For the diffusion term, a second-order discretisation is applied:

(2.20) \begin{equation} D_s \frac {\partial ^2 \Gamma }{\partial s^2} \approx D_s \frac {\Gamma _{I+1} - 2\Gamma _{I} + \Gamma _{I-1}}{(\Delta s)^2}. \end{equation}

The sum of (2.20) over all $N_{ {cells}}$ gives the total diffusion $\eta _{ \textit{diff.}}$ at each time. Together, the discretised form of the transport equation is

(2.21) \begin{equation} \frac {\Gamma _{I}^{n} - \Gamma _{I}^{n-1}}{\Delta t} + u_I^n \frac {\Gamma _{I+1}^n - \Gamma _{I}^n}{\Delta s} + \Gamma _I^n \frac {u_{I+1}^n - u_{I}^n}{\Delta s} + \Gamma _{I}^n \kappa _I^n v_I^n = D_s \frac {\Gamma _{I+1}^n - 2\Gamma _{I}^n + \Gamma _{I-1}^n}{(\Delta s)^2}. \end{equation}

Assuming a zero source term (i.e. $j_n = 0$ ), the discretised form of the transport equation, with like terms grouped together in the form of a matrix $a\; \Gamma = b$ , is given by

(2.22) \begin{align} \Gamma _{I+1}^n \overbrace {\left [ \frac {\Delta t\, u_I^n}{\Delta S} - \frac {D\, \Delta t}{(\Delta S)^2} \right ]}^{{a_{I+1}^n}} &+ \Gamma _I \overbrace {\left [1 + \Delta t \left ( \frac {u_{I+1}^n - 2u_I^n}{\Delta s} + \kappa _I^nv_I^n + \frac {2D}{(\Delta s)^2} \right ) \right ]}^{{a_I^n}} \nonumber \\ &{}+ \Gamma _{I-1}^n \overbrace {\left [ -\frac {D\, \Delta t}{(\Delta S)^2} \right ]}^{{a_{I-1}}} = \Gamma _I^{n-1}. \end{align}

The system of equations is solved using the inbuilt linear algebra solver in Numpy (v.1.26.0). In solving the transport model, the domain assumes symmetry about $x=0$ , and hence is discretised from $x=0$ to $x=l$ along the positive $x$ -axis. Note that the continuum equation is solved in $s$ coordinates following the surface; $x$ and $s$ are equivalent, and the only difference arises in quantifying the length travelled by the molecules, i.e. $\Delta s \neq \Delta x$ . As such, the distance along the surface is worked out from $(\Delta s)^2 = (\Delta x)^2 + (\Delta z)^2$ . The following set of initial and boundary conditions and the surface fit provides a closed set of equations to solve (2.22):

(2.23) \begin{align} & \Gamma (s, t=0) =\overset {MD}{\Gamma } (s, t=0), \end{align}

(2.24) \begin{align} & \Gamma (s=0,t) = \alpha\, {\rm e}^{-\lambda t} + \alpha _0, \end{align}

(2.25) \begin{align} & \Gamma (s=l,t) =0, \end{align}

(2.26) \begin{align} & \boldsymbol {e}_n(s,t) = - \frac {\boldsymbol{\nabla} _s \zeta (s,t)}{\|\boldsymbol{\nabla} _s \zeta (s,t)\|}, \end{align}

(2.27) \begin{align} & u_s (s,t) = \overset {MD}{{u}_s} (s,t), \end{align}

(2.28) \begin{align} & D_s = D_0\,\mbox {e}^{-\beta \Gamma }. \end{align}

The initial concentration field is set to the MD data at $t=0$ in (2.23) obtained using (2.15); the left-hand boundary condition (2.24) is modelled based on a fit to the data from MD at $s=0$ . As surfactant molecules do not reach the right-hand end of the film during the simulation time, zero concentration is assumed at the far right (2.25). A careful tracking of the spatially and temporarily evolving surface to get the local properties is necessary to obtain an accurate molecular description of the transport process.

To obtain the local curvature of the surface, as used in (2.26), the deforming liquid–vapour surface is identified from the MD simulation. This uses the fluid density field $\rho$ of water molecules, with a threshold to identify the higher-density part as liquid before implementing edge detection. The surface is then fitted to a spline using the UnivariateSpline function from the scipy library to obtain a functional $\zeta (s,t)$ ; see § S1 in the supplementary material. Using the derivatives of the spline, the distance ${\rm d}s$ as well as the normal and tangential vectors are obtained at each time step. These derivatives of the surface at each point are used to get $K_I^n$ from the definition $\kappa (x) = {\zeta ''}/{ [1 + \zeta '^2 ]^{3/2}}$ , where $\zeta '$ and $\zeta ''$ represent the first and second derivatives of the spline fit with respect to $x$ , respectively. The local velocities obtained from MD simulations are resolved along the surface, as in (2.27). Traditionally, the mass transport equation is solved alongside a thin film equation to capture surface deformations and their effect on the transport dynamics. In this study, this is bypassed by directly extracting the evolving surface shape and surfactant velocities from MD simulations, which inherently captures this deformation at the molecular scale without the need for continuum approximations.

Surfactant molecules on a curved surface can occupy a larger area compared to a flat surface, potentially slowing down the evolution of the concentration gradient. Additionally, the paths available for diffusion are geometrically constrained by curvature, leading to anisotropic diffusion (Shu & Gong 2001; Ramírez-Garza et al. 2021). In a dedicated series of simulations focusing on uniformly distributed surfactants on a water film, the surface diffusion of SDS molecules was estimated. To distinguish surface diffusion from bulk diffusion, only the surfactant molecules that remained on the film’s surface were tracked. The trajectories of these molecules were analysed to calculate the mean squared displacement (MSD) over the simulation period, i.e. $\mathrm {MSD} = \langle |s(t_0-t) - s(t_0)|^2 \rangle$ , where $s(t)$ is the position of a molecule at time $t$ , and $s(t_0)$ is its initial position. Subsequently, we employed Einstein’s equation for surface diffusion, $D_s = \mathrm {MSD}/t$ (Rideg et al. 2012); the methodology employed here to obtain the diffusion coefficient from molecular simulations has been detailed in § S2 in the supplementary material. The diffusion coefficient obtained shows a concentration dependence as in (2.28), i.e. $D_s = D_0\,\mbox {e}^{-\beta \Gamma }$ , where $D_0 = 3.22\, \mathrm {nm}^2\,\mathrm{ns}^{-1}$ , $\Gamma$ is the local number of surfactant molecules per unit area, and $\beta =4/3$ . For dimensional consistency, $\beta$ should have units of area.

The reduced mobility of the surface molecules at higher concentration leads to a lower rate of diffusion. Similar reduction in diffusion rate due to particle crowding was observed in previous studies. A longer tail length to achieve the diffusive limit was also observed (Rideg et al. 2012; Tomilov et al. 2012). This functional relation employed in solving (2.1), which originally considered a constant diffusion coefficient, accounts for variations in local diffusivity caused by the inhomogeneous local surfactant concentration. The effective diffusion coefficient, which updates ‘on the fly’ based on the local concentration, reflects an inherent nonlinearity in the governing transport equations. The inclusion of these features is crucial to model accurate transport in systems with spatial- and time-dependent ‘driving’ properties. This concentration-dependent diffusion coefficient suggests that the diffusive term in (2.1) is essentially $\boldsymbol{\nabla} _s \boldsymbol{\cdot} ( D_s(\Gamma )\, \boldsymbol{\nabla} _s \Gamma )$ . Note that the diffusion coefficients obtained in our simulations are higher than the experimentally measured bulk diffusion of SDS in aqueous solutions (Ribeiro et al. 2003), by a factor of $2$ – $3$ . This can be attributed to the speed-up process of the kinetics in coarse-grained models that contributes to the over-estimation of diffusion (Vögele et al. 2015; Schmitt et al. 2023). In Martini models, the best estimated kinetic speed-up factor is approximately $4$ ; however, the polarised water model has been reported to show approximately a $2.5$ times speed-up (Marrink & Tieleman 2013). To overcome any potential mass imbalance arising from numerical diffusion or discretisation errors, a conservation correction mechanism is integrated to ensure that the total quantity of the transported scalar is preserved throughout the simulation, i.e. $\int _0^l \Gamma\, {\rm d}s = \mathrm {constant}$ .

Since any derivative property on the molecular scale is noisy, the data were smoothed without significantly distorting the higher moments of the signal. This was achieved by utilising a Savitzky & Golay (1964) filter, which smooths input signals by fitting a polynomial using least squares techniques. These facilitate accurate calculation of the spatio-temporal surface vectors, i.e. $\boldsymbol{e}_t$ and $\boldsymbol{e}_n$ , which are instrumental in resolving the velocities of surfactant molecules parallel and normal to the interface. Note that all properties are averaged along the $y$ direction.

3. Results and discussions

3.1. Control volume conservation

The equations derived in § 2.3 track exactly the molecules in a control volume on the surface. In this subsection, the control volume (2.6), as discretised in (2.13), was tested for a moving interface, showing exact mass conservation. Figure 4(a) illustrates the process of establishing the system of control volumes by fitting a second-order spline to the outermost surfactant head group in each volume. This is shown in figure 4(a) at a single arbitrary time, with the corresponding mapping from ( $x,z$ ) space to ( $\chi , \omega$ ) space shown in figure 4(b). Figure 4(c) shows the control volume balance for volume 14 (which is second from the right) as it contains only a single molecule and clearly highlights the various terms in (2.13). Once the molecular positions in $(\chi , \omega)$ space are obtained, the flux terms are calculated by simply counting if a molecule leaves one volume and joins another. These fluxes are shown in figure 4(c) as circles. These are almost equal to the corresponding change in the number of molecules in a volume over time, as shown by the black line. Any difference between the fluxes and change in time is due to the stretching of the volume as the surface moves. This change due to surface movement is shown as a dashed red line in figure 4(c), while the inset shows a period of time when there are no surface fluxes, highlighting the effect on measured surface concentration of the deforming volume. The contour in figure 4(b) shows the Jacobian at each point in a given volume, due to the shape of the surface. With the convex shape of the surface seen in figure 4(a), the result is a stretching at the top of the volumes, and a squeezing at the bottom (it would be reversed for a concave surface), where these values are normalised by the ratio of unmapped to mapped volumes $\Delta x\, \Delta y / (\Delta \chi\, \Delta \omega )$ , to highlight the small effects of the surface curvature.

Figure 4. (a) Surfactant head groups at the liquid–vapour interface corresponding to (non-dimensional) time $t=17$ , fitted by a spline (dashed black line) that is shifted vertically by $\pm \Delta y=25$ to define the top and bottom surfaces of the control volumes, to represent the surface along which the surfactant molecules move. The surface is divided into 15 control volumes, each of width $\Delta x \approx 25$ , separated by the line normal to the middle spline at intervals of $\Delta x$ . The surfactant molecules are colour-coded according to the Jacobian used for mapping, except cell 8, which is highlighted in green. (b) The mapped coordinate system in terms of $\chi , \omega$ , with the Jacobian as a contour in the background. Each volume has the same width and height $\Delta \chi = \Delta \omega = 2$ , so the Jacobian contour is shown normalised by the ratio of unmapped to mapped volume for the flat surface case, $(\Delta x\, \Delta y) /(\Delta \chi\, \Delta \omega )$ , so shows squeezing at the volume top $J\lt 1$ and stretching at the bottom $J\gt 1$ . (c) The control volume balance equation (2.6) is shown for volume 14. The grey circles are surface fluxes $\sum _{i=1}^N m_i J_i \dot{\boldsymbol {s}}_i \boldsymbol{\cdot} {\rm d}\textbf {S}_i$ , the black lines are the control volume time evolution ${{\rm d}}/{{\rm d}t} \sum _{i=1}^N m_i J_i \vartheta _i$ , and the red dashed line is the sum of the molecular Jacobian terms $\sum _{i=1}^N m_i \dot {J}_i \vartheta _i$ . The sum of these three terms is zero at all times and is shown as the thin black line, which serves as the abscissa. No absorption $j_n$ is observed for this volume, but may occasionally occur in other volumes.

One final noteworthy point is that at time $t=21$ , when a molecule is seen to simultaneously leave and enter, the net change is zero except for the difference due to the Jacobians of the leaving and entering molecules. For all volumes, including those containing a single molecule, the control volume equations derived in § 2.3 are seen to be satisfied exactly. The horizontal black line in figure 4(c) shows the sum of all terms, which is observed to be zero to near machine precision. This approach was tested for all 15 control volumes over all 50 time units shown, and (2.13) is respected for all cases. In some volumes, reabsorption $j_n$ to the bulk is observed; for example, near cell 10 in figure 4(a), a few molecules can be seen below the bottom surface of the control volumes. Given that surfactant molecules tend to stay at the surface, these events were followed by an absorption event at a later time. Nevertheless, the framework would be able to capture those events in cases where significant absorption would be expected (e.g. past the CMC when micelles would form). The size of the control volume $\Delta y$ could be chosen to ensure that these events do not occur unless they represent a systematic concerted process. The full procedural details for obtaining the spline and defining these mappings from a system of bilinear finite elements determined using the spline fit have been elaborated in § S4 in the supplementary material, with the link to open-source code provided in ‘Code availability’.

3.2. Surfactant transport

Figure 5 examines the time evolution of the local surfactant concentration $\Gamma$ (normalised by the initial concentration at $x=0$ , i.e. $\Gamma _0$ ) at different times during the spreading process. This considers the case when there are $0.35$ nanogram of surfactant molecules deposited on a planar film. The finite difference solution of (2.1) is depicted by the solid line, and the symbols denote corresponding data obtained directly by MD using (2.15) to average. The close agreement between these two quite different procedural routes affirms the validity of the continuum model in capturing the nanoscale transport mechanism, at least for the simple SDS surfactant case. The dotted lines represent the spline fit $\zeta$ to the instantaneous liquid–vapour surface obtained from the MD systen. Section S3 of the supplementary material includes similar plots for different amounts of deposited surfactants, illustrating the robust effectiveness of the transport equation employed. It is observed that results of MD simulations with an initially denser monolayer align more closely with the numerical solution of the transport equation. The reasons for this are not obvious, but could be due in part to uncertainties associated with estimating the diffusion coefficient in these surfactant-deficient systems, and a greater role of fluctuation-related cooperative effects. For less dense monolayers, local surfactant concentration fluctuations are likely, which would not be captured by the mean-field continuum model.

Figure 5. Spatio-temporal variation of the normalised surfactant concentration along the length of the film at distinct times, where $x=0$ denotes the centre of the film. The solid line represents surfactant evolution predicted by the continuum transport equation (2.1), while the symbols denote data from the molecular simulations. The dotted line corresponds to spline fitting to the evolving surface. The arrows show the local tangential ( $\boldsymbol{e}_t$ ) and normal ( $\boldsymbol{e}_n$ ) surface vectors (not to scale), and $l$ denotes the spreading length.

In addition to the spreading of surfactant across a planar film, the analysis was extended to an initial globally curved geometry, i.e. a droplet as in figure 1(b). Although a droplet has a higher global curvature than the planar films, figure 6 shows that the transport mechanism remains largely unaffected. For planar films, the furthest right end of the film remained devoid of surfactant, resulting in a persistent surface tension gradient that favoured continuous transport throughout the MD simulation. In contrast, in the droplet case, surfactant molecules from both the upper and lower regions of the droplet readily redistributed, reducing the fractional area of surfactant-poor regions. Consequently, the transport process quickly approached a state where surfactants from both halves of the droplet came into close proximity, which significantly diminished the driving force for further transport, and nearly halted any directional flow of the surfactants. Note that the solid lines in the figure, representing the solution of the transport model, show some noticeable non-smoothness, which can be attributed to the inherent fluctuations in the MD simulations, injected into the finite difference solver through the MD derived boundary conditions, surface curvature and velocity field. Fluctuations in the MD data are a natural feature of liquids when viewed on molecular dimensions, particularly in the velocity fields and surface deformations. Additionally, the evolving surface curvature introduces fluctuations in the transport process. While the above results demonstrate that the advection–diffusion model provides an effective description of surfactant transport within molecularly thin films for the considered initial conditions and geometries, a more comprehensive investigation would be necessary to assess the robustness of the continuum model under a broader range of surfactant systems of differing chemistry. In some cases, the chemistry of the surfactant molecules may induce changes in their distribution during spreading (such as phase or ordering changes) that are beyond the scope of the continuum description to capture (at least as currently formulated). Future studies could explore the transport dynamics of, for example, insoluble surfactants, employing different coarse-grained representations of water, and investigate a variety of solvent–surfactant hydrophobic–hydrophilic interactions. Nevertheless, the present work is still a necessary and useful starting point for future research to facilitate theoretical advances in this field.

Figure 6. Temporal variation of the surfactant concentration on a droplet. Data from MD simulations are represented by symbols, while the solid lines indicate the numerical solutions to the transport equation.

The diffusion coefficients of the surfactant molecules are higher for low initial concentrations. At higher concentrations, the diffusion rate of the surfactant molecules diminishes and becomes a highly collective process resulting from the coupled interactions between surfactant molecules in close proximity. This suggests that for initially dilute surfactant layers, the self-diffusional contribution to the spreading event, $\eta _{\textit{diff.}} (t)$ , might have a relatively greater impact on the overall spreading event. The remaining component, i.e. contributions from a ‘geometric’ term, $\eta _{\textit {geom.}} (t)$ might be consistently insignificant in all cases, possibly because both $\kappa _I$ and $v_I$ fluctuate around zero, and are small numbers, rendering their product an order of magnitude smaller than the advective and diffusive terms. They might gain more importance for surfaces with higher local curvatures, or when the surface is on a growing/shrinking droplet, film or bubble. These contribution terms are plotted in figure 7(a).

Figure 7. (a) Relative contributions of advection, diffusion and geometrically induced transport. (b) Spreading of the monolayer. Symbols denote MD data, with empty and solid symbols denoting independent simulations. Solid lines show $l\sim t^{1/2}$ fitting to the data. The shaded region shows the $90\,\%$ confidence interval of the fitting. (c) Zoomed-in view of the initial rapid transient period (boxed in (b)) showing the early time (linear) spreading.

3.3. Spreading length

A useful parameter that characterises the surfactant transport across the film is the ‘spreading rate’, which denotes the rate at which the front of the initial surfactant layer covers the liquid–vapour interface. The surfactant monolayer manifests a spatial difference in the surface tension across the interface, the magnitude at every point in space being a function of the local surfactant concentration. This gradient in tension along the surface is a contributing factor in driving the spreading of the monolayer. For spreading, surfactants must overcome the viscous drag from the water sub-phase and other surfactant molecules. By balancing these two forces, Ahmad & Hansen (1972) derived a scaling law for the spreading length $l$ , which was found to have the simple analytic form $l\sim t^{1/2}$ . Despite such simplicity, the authors found that this expression represented the experimental data very well. Borgas & Grotberg (1988) derived a similar power-law dependency between $l$ and $t$ , using a more formal analysis route that involved solving the Navier–Stokes equations, including the lubrication approximations, taking into account the surface tension gradients, and the viscous dissipation within the thin liquid film underlying the monolayer.

Figure 7(b) shows that the spreading of the surfactant monolayer in the present simulations follows the characteristic square root dependency on time for a range of initial concentrations. This observation agrees with past experimental studies; for example, Hanyak, Sinz & Darhuber (2012) studied the spreading of SDS on deformable liquid surfaces, and obtained an average exponent $0.48$ for a range of concentrations. Similar power-law dependency was observed in several other studies (Borgas & Grotberg 1988; Hamraoui et al. 2004; Howell & Stone 2005; Hernández-Sánchez et al. 2015). The scaling exponent may assume a range of values, depending on the geometry of the problem, and the physicochemical properties of the liquids involved, such as viscosity, surface tension, miscibility and surface contamination (Jensen 1995; Santiago-Rosanne, Vignes-Adler & Velarde 2001; Rafaı & Bonn 2005; Lee & Starov 2007). While the $t^{1/2}$ rate dependence persists over the whole simulation time window in the present study, the initial spreading (for approximately $t\lt 0.5\, \mathrm {ns}$ ) appears to be faster with the exponent $n\to 1$ ; see figure 7(c) for a zoomed-in view of the initial spreading period. Similar accelerated initial spreading was also observed in several previous experimental studies of monolayer spreading (Lee & Starov 2007; Mitra & Mitra 2016; Motaghian et al. 2022). This relatively rapid spreading, or more generally diffusion at short times, is also observed for simple bulk liquids, where it is usually attributed to initial ‘ballistic’ molecule trajectories. It is an intrinsic feature at molecules on the molecular scale (Alder & Wainwright 1967), although whether the origin is still the same for these molecular systems is not entirely clear, and would be worthy of further investigation. The initial rapid spreading may, for example, be linked to the development of shear flow, where starting from a quiescent field, the surfactant might spread quickly in the beginning, driven by advection. As the shear flow starts to develop, spreading slows down. It has been found that a ballistic to diffusive crossover is expected as interactions and scattering of the particles attain greater importance (Huang et al. 2011).

4. Conclusions

This study uses empirical data from NEMD simulations to parametrise a continuum model of surfactant transport, confirming its efficacy at the molecular scale. Spreading is influenced not only by the boundary conditions, but also by the chemical nature, and the molecular structure and solubility of the surfactant. The inherently molecular nature of the surfactant spreading process on a deforming liquid–vapour surface makes NEMD simulations with molecular resolution an ideal choice for investigating these phenomena. This can compliment continuum methods, where approximations such as diffusion, surface tension and viscosity have to be made to ‘incorporate’ surfactant chemistry.

An NEMD methodology for modelling surfactant spreading on a thin liquid film was developed in this work. The NEMD simulations were carried out for an SDS monolayer patch deposited on a thin film of liquid water. The spreading of the monolayer across the surface and subsequent curved deformation of the film were followed. Local surfactant concentration profiles as a function of time were computed. The spreading length of the surfactant patch was followed and interpreted in terms of a power law. The MD results of the surfactant transport process were compared with the coupled convection–diffusion transport equation. This transport equation describes the dynamical processes that arise naturally during the spreading event. The comparison was facilitated by developing a novel coarse-graining of the surfactant molecules using a Jacobian transform of the control volume formulation. The continuum-level calculations were started from the initial coarse-grained density profile of the NEMD system, sharing the NEMD boundary conditions and velocity field to test the evolution of surfactant density with diffusion coefficients calculated by separate equilibrium bulk water–surfactant MD simulations.

The NEMD-informed continuum model accurately describes the nanoscale transport process for the particular surfactant modelled here. Although the continuum transport equation is reasonably accurate for this prototypical SDS surfactant on thin water films, it may require further modifications for other surfactant–solvent systems. Using the proposed control volume framework, NEMD can now be employed as a valuable tool to complement continuum methods by capturing the molecular details of important transport processes.