2.1. Molecular description of the system

A coarse-grained modelling approach is employed. The MARTINI polarisable water model (Yesylevskyy et al. Reference Yesylevskyy, Schäfer, Sengupta, Marrink and Levitt2010) 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. Reference Weiand, Rodriguez-Ropero, Roiter, Koenig, Angioletti-Uberti, Dini and Ewen2023); 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. Reference Martínez, Andrade, Birgin and Martínez2009) and Moltemplate (Jewett et al. Reference Jewett2021). 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. Reference Thompson2022). 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 Reference Gaver and Grotberg1990; Stone Reference Stone1990), 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 Reference Wong, Rumschitzki and Maldarelli1996; Cermelli, Fried & Gurtin Reference Cermelli, Fried and Gurtin2005; Pereira et al. Reference Pereira, Trevelyan, Thiele and Kalliadasis2007). Among these studies, Wong et al. (Reference Wong, Rumschitzki and Maldarelli1996) 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. (Reference Wong, Rumschitzki and Maldarelli1996), 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 Reference Manikantan and Squires2020). 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 (Reference Irving and Kirkwood1950). 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 (Reference Irving and Kirkwood1950), 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. Reference Smith, Heyes, Dini and Zaki2012), 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. Reference Smith, Heyes, Dini and Zaki2012), 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. Reference Thompson2022) 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. Reference Smith, Heyes, Dini and Zaki2012). 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. (Reference Wong, Rumschitzki and Maldarelli1996) 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 (Reference Irving and Kirkwood1950), 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 (Reference Stone1990), there is no obvious diffusion contribution obtained in the molecular system. This is consistent with the work of Scriven (Reference Scriven1960), 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 Reference Chapman and Cowling1990).

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 Reference Pototsky and Suslov2024; Yang et al. Reference Yang, Du, Xiao, Wang, Zhao and Xiong2024). 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. Reference Wong, Rumschitzki and Maldarelli1996), and so is consistent with Scriven (Reference Scriven1960).

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 Reference Todd, Evans and Daivis1995). However, the statistics are known to be poor for surface flux measurements (Smith, Heyes & Dini Reference Smith, Heyes and Dini2017), 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. Reference Heyes, Smith, Dini and Zaki2011), 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 Reference Hirsch2007) 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 Reference Shu and Gong2001; Ramírez-Garza et al. Reference Ramírez-Garza, Méndez-Alcaraz and González-Mozuelos2021). 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. Reference Rideg, Darvas, Varga and Jedlovszky2012); 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. Reference Rideg, Darvas, Varga and Jedlovszky2012; Tomilov et al. Reference Tomilov, Videcoq, Chartier, Ala-Nissilä and Vattulainen2012). 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. Reference Ribeiro, Lobo, Azevedo, Miguel and Burrows2003), 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. Reference Vögele, Holm and Smiatek2015; Schmitt et al. Reference Schmitt, Fleckenstein, Hasse and Stephan2023). 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 Reference Marrink and Tieleman2013). 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 (Reference Savitzky and Golay1964) 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.