Memory Effects in Contact Line Friction

Abstract

When a drop of liquid comes into contact with a solid surface, it relaxes towards an equilibrium configuration, either wetting the surface or remaining in a droplet-like shape with a finite contact angle. The force driving the process towards equilibrium is the corresponding out-of-balance Young’s force. However, the speed with which the liquid front advances depends strongly on an opposing friction force arising from dissipative processes due to the moving solid-liquid-gas contact line. In analogy to the treatment of hydrodynamic friction we present an exact method, based on the Mori–Zwanzig formalism, to extract this friction from equilibrium data. We find that the contact line exhibits long-lasting memory with a characteristic power-law decay due to coupling to the systems hydrodynamic modes. Within linear response regime, we obtain the frequency-dependent dissipative and elastic response of the contact line to an external perturbation, including a frequency-dependent friction coefficient. Similar to hydrodynamic friction in liquids, we find that the friction decreases beyond a characteristic frequency and the system exhibits predominantly elastic behavior.

I Introduction

Theories that describe dynamic wetting phenomena generally consider the microscopic motion of the contact line and the associated contact line friction to be Markovian. An overlap between the relaxation of the liquid front and other microscopic relaxation processes is not considered. For theories that employ continuum mechanics to describe dynamic wetting phenomena,[hoffman_study_1975, voinov_hydrodynamics_1976, hocking_moving_1977, cox_dynamics_1986, pismen_mesoscopic_2002], this is a natural choice. However, molecular kinetic theory (MKT)[blake_kinetics_1969, blake_influence_2002] (and related approaches,[brochard-wyart_dynamics_1992, petrov_combined_1992]), which attribute the contact-line motion to molecular adsorption at the liquid front, similarly rely on a Markovian assumption. From a comparison to hydrodynamic friction in fluids it could be argued that a Markovian description of the contact line starts to break down only at time scales comparable to the relaxation of molecular degrees of freedom.[mason_optical_1995, gittes_microscopic_1997, squires_fluid_2010, straube_rapid_2020]

Herein, we demonstrate that the Markovian assumption breaks down at least one order of magnitude earlier, as the movement of the contact line and the contact line friction are dominated by persistent and long-lasting memory effects. We introduce a framework to obtain the full frequency-dependent contact-line friction, utilizing information from non-critical equilibrium fluctuations of the contact line, reducing to the approach by Toledano et al.[fernandez-toledano_contact-line_2019, fernandez-toledano_moving_2020] in the Markovian limit. Within linear response theory, we then obtain the response of the contact line to an external force. The dynamics of the contact line at short times are dominated by low-frequency molecular vibrations, with no modes that could be attributed to MKT absorption modes. On long time scales, the contact-line correlation functions are governed by a power-law decay, due to coupling to the system’s hydrodynamic modes. Notably, the majority of the contact-line friction emerges due to this coupling.

II Theoretical Background

At equilibrium, the equation of motion of the contact line may be written as a generalized Langevin equation (GLE)[mori_transport_1965] in terms of the contact line position $X$ and velocity $V$

$$m\frac{\differential V_{t}}{\differential t}=F_{t}=-kX_{t}-\int_{0}^{\infty}\differential s\,K(s)V_{t-s}+F^{R}_{t},$$ (1)

where $K$ is the memory kernel, $F^{R}$ is a random force, $F$ is the force acting on the contact line, $k$ is a harmonic force constant and $m$ is the mass of the contact line.

Figure 1: (a) A simulation snapshot of the water–graphene liquid bridge viewed in the $xz$-plane. (b) A schematic illustration of how the density field is calculated to obtain a continuously changing field in time. The black dots indicate the center of each grid cell with cell length $l$ and $\Delta x$ and $\Delta z$ the distance of the particle to the lower left cell along the $x$- and $z$-axis respectively. The contribution of a particle (blue circle) to the density field is distributed to the four nearest grid cells via bilinear interpolation, e.g., the contribution to the bottom left cell is given by $(l-\Delta x)(l-\Delta z)/l^{2}$. (c) The instantaneous density profile $\rho(x,z)$ in the $xz$-plane and (d) its ensemble average $\expectationvalue{\rho(x,z)}$. (e) The density profile along $x$ at $z_{\mathrm{Cl}}$ for the instantaneous density (dotted blue line) and its time average (dashed blue line). Both profiles are approximated with Eq. 7 (solid red lines) to determine the contact line position from the interfacial position. (f) The resulting trajectory of the contact line position $X$.

We review a derivation of Eq. 1 in supplementary information (SI) S1. We denote functions of time with brackets, e.g., $K(t)$, and the time dependence of observables along a microscopic trajectory by subscript $t$, e.g., $V_{t}$. We include a harmonic force in Eq. 1, because if the droplet is at rest, the contact line will fluctuate around some average position $\expectationvalue{X}=0$. While the memory kernel can be obtained from projected force correlation functions,[carof_two_2014] this is usually not done. Multiplying the Eq. 1 with $V_{0}$ yields a Volterra equation

$$m\dot{C}(t)=-k\expectationvalue{V_{0},X_{t}}-\int_{0}^{\infty}\differential s\,K(s)C(t-s),$$ (2)

where $\expectationvalue{\dots}$ denotes an ensemble average and $C(t)=\expectationvalue{V_{0},V_{t}}$ is the velocity auto-correlation function (VACF). This equation can be inverted[straube_rapid_2020, kowalik_memory-kernel_2019] to extract the memory kernel from the corresponding VACF. We want to point out that Eq. 1 is exact if the probability distribution of the contact line position around its average is Gaussian, i.e., the potential of mean force is harmonic.

For a stationary non-equilibrium ensemble under an external perturbation $P$, it follows that[kubo1991]

$$m\frac{\differential}{\differential t}\expectationvalue{V}_{1}(t)=P(t)-\int^{\infty}_{-\infty}\differential s\,K_{+}(t-s)\expectationvalue{V}_{1}(s),$$ (3)

where the subscript $+$ indicates multiplication with the unity step function and $\expectationvalue{\dots}_{1}$ denotes the average over the non-equilibrium ensemble. Fourier transforming this expression and rearranging for the velocity yields

$$\expectationvalue{\hat{V}}_{1}(\omega)=\frac{1}{-i\omega m+\hat{K}_{+}(\omega)}\hat{P}(\omega)=\hat{\Upsilon}(\omega)\hat{P}(\omega),$$ (4)

where $\hat{\Upsilon}(\omega)$ is the (linear) response function of the contact line. Here $\real(\hat{\Upsilon}(\omega))$ is the dissipative response of the contact line and $\imaginary(\hat{\Upsilon}(\omega))$ is the elastic response. The frequency-dependent contact line friction $\zeta$ follows from the dissipative response

$$\zeta(\omega)=\frac{1}{L}\real(\hat{\Upsilon}^{-1}(\omega))=\frac{1}{L}\real(\hat{K}_{+}(\omega)),$$ (5)

where $L$ is the length of the contact line. Computing the memory kernel via Eq. 2 usually requires correlation functions with a high time resolution, but the zero-frequency contact line friction coefficient can be obtained from the position correlation function at a much lower resolution with a Green-Kubo-like relation[kowalik_memory-kernel_2019]

$$\zeta(0)=\frac{k_{\mathrm{B}}T}{L}\int^{\infty}_{0}\differential t\,\frac{\expectationvalue{X_{0},X_{t}}}{\expectationvalue{X_{0},X_{0}}^{2}},$$ (6)

where $k_{\mathrm{B}}$ is the Boltzmann constant and $T$ is the temperature. We review a derivation of Eq. 6 in SI S2.

III Extracting the Contact Line Position

Figure 2: (a) the VACF and (c) the corresponding memory kernel $K$ for an LJ-octamer fluid confined between generic solid plates, and (b) the VACF and (d) the corresponding memory kernel for water confined between graphene sheets. For both systems, the VACF and, to a lesser extent, the memory kernel are well captured by a short-time expansion of the form $1-(t\omega_{0})^{2}$ (dotted black lines) with the slope approaching zero for $t\rightarrow 0$. This indicates that the contact line velocity is a well-defined observable. At long times, both the VACF and memory kernel decay according to a power law (dashed black lines) with an exponent of $2/3$.

We study two liquid bridge systems: A generic Lennard-Jones (LJ) fluid, where each molecule consists of eight beads, with negligible vapor pressure confined between two simple cubic plates with a contact angle of 113.0°$\pm$0.2°, and water confined between graphene sheets with a contact angle of 80.3°$\pm$0.6°, shown in Fig.1 (a). Details on the molecular dynamics (MD) simulations are provided in the Appendix and details on the contact angle extraction in SI S3. Projecting the system onto the $xz$-plane reveals four contact lines aligned along the $y$-axis. To determine their positions, we compute the density field $\rho(x,z)$ on a grid with cell length $l$. We then follow common approaches used in, e.g., hybrid particle–field theory [milano_hybrid_2009] to distribute the contribution of each particle to the surrounding cells by bilinear interpolation, schematically shown in 1 (b). This ensures that the contact line position changes continuously in time (see SI S4). The resulting density profile and its ensemble average are shown in Fig. 1 (c) and (d), respectively. We assume the $z$-position of the contact line $z_{\mathrm{Cl}}$ to be at the first density maximum normal to the plate and extract the position along $x$ from the density profile $\rho(x,z_{\mathrm{Cl}})$ (dotted blue line) and its ensemble average (dashed blue line) shown in Fig. 1 (e). Using that both systems have a negligible vapor pressure, we approximate the density profiles (solid red lines) with,

$$\frac{\rho(x,z_{\mathrm{Cl}})}{\langle\rho(x_{\mathrm{L}},z_{\mathrm{Cl}})\rangle}=\frac{1}{2}+\frac{1}{2}\tanh(2\frac{x-x_{\mathrm{Cl}}}{w}),$$ (7)

where $w$ is the interfacial width and $x_{\mathrm{L}}$ denotes a position where the density at the solid-liquid interface is not influenced by the liquid-gas interface. Finally, we obtain the contact line position as the difference between the instantaneous and the average interfacial position. This procedure yields a smooth contact line trajectory shown in Fig. 1 (f), with a Gaussian probability density for both the LJ- and the water–graphene system (see SI S5), from which the contact line velocity is obtained by numerical differentiation. Finite size effects with respect to the length of the contact line do not play a major role (see SI S6).

IV Non-Markovian Effects and Frequency Dependent Response

Figure 3: The dissipative $\real(\hat{\Upsilon})$ (solid red lines) and elastic $\imaginary(\hat{\Upsilon})$ (solid purple lines) response of the contact line of a) an LJ-octamer fluid confined between generic solid plates, and b) water confined between graphene sheets. The real part of the half-sided Fourier transform of the memory kernel $\real(\hat{K}_{+})$, corresponding to the frequency-dependent contact line friction $\zeta(\omega)$ (solid yellow lines), and the imaginary part of the transform $\imaginary(\hat{K}_{+})$ of c) an LJ-octamer fluid confined between generic solid plates, and d) water confined between graphene sheets.

Figure 2 shows the VACF $C(t)$ and the corresponding memory kernel $K(t)$, obtained by inverting Eq. 2, on a logarithmic scale for (a,c) a generic LJ liquid bridge and (b,d) a water–graphene liquid bridge. We can distinguish three regimes labeled I-III. At short times (regime I), the VACF -and to a lesser extent the memory kernel- follow a quadratic expansion $1-\nu^{2}t^{2}+\mathcal{O}(t^{4})$ (black dotted lines). Odd-order derivatives vanish due to the skew symmetry of time derivatives. The excellent agreement of the VACF with the quadratic expansion confirms that the contact line velocity is a well-defined observable. By contrast, $K(t)$ deviates more strongly, which we attribute to numerical limitations. Following the initial decay, both the VACF and the memory kernel show damped oscillations (regime II). We shown in SI S7, that these oscillations coincide with low–frequency molecular vibrations arising from fast, molecular-scale motions. For the water–graphene contact line, in particular, we observe a strong peak at $1.2$ THz in the Fourier transform of the VACF coinciding with SPC/Es low-frequency librations.[carey_measurement_1998, lin_two-phase_2010] At large times (regime III), both the VACF and the memory kernel follow a power law decay with an exponent of $-2/3$. The shared asymptotic behavior follows from the Laplace transform of Eq. 2.[corngold_behavior_1972] Power-law decays in correlation functions are generally a result of coupling to slow modes of the system. Note that if the movement of the contact line were determined by fast processes, the decay of correlations would be exponential. We show in SI S8 that all four contact lines couple weakly to each other, indicating that the contact lines couple to hydrodynamic modes. This coupling results in corresponding finite size effects (see SI S8).

Fig 3 shows the dissipative response $\real(\hat{\Upsilon})$ (solid red line) and the elastic response $\imaginary(\hat{\Upsilon})$ (solid purple line) given by Eq. 4 for (a) the LJ-fluid and (b) the water-graphene liquid bridge. For further details on the computation of the Fourier transform, see SI S9. In close analogy to hydrodynamic friction, at low frequencies, the response is static and dissipative. In this case, the response is well described by a Markovian contact line friction coefficient $L_{y}\hat{\Upsilon}(\omega\to 0)=1/\zeta(0)$ (upper gray dashed line), and there is no elastic response. At high frequencies beyond the characteristic frequency of the contact line, corresponding to a maximum in the elastic response, dissipation decreases exponentially to zero (lower gray dashed line), and the contact line behaves elastically. In this limit the response follows the asymptotic $L_{y}\hat{\Upsilon}(\omega\to\infty)=iL_{y}/(m\omega)+\mathcal{O}(\omega^{-3})$ (black dotted line). Fig. 3 shows the frequency-dependent contact line friction $\zeta(\omega)=\real(\hat{K}_{+}(\omega))$ (yellow solid line) and the corresponding imaginary part $\imaginary(\hat{K}_{+})$ (solid green line) for (c) the LJ-fluid and (d) water-graphene liquid bridge. At low frequencies, the contact line friction approaches the above-mentioned Markovian contact line friction coefficient, followed by a decrease in friction. The time scales in this frequency range fall into regime III of the VACF and memory kernel, meaning the majority of the contact-line friction emerges on the same time scales where the contact line couples to slow modes. At high frequencies, beyond the characteristic frequency of the system, the contact line friction decreases since the real part of $\hat{\Upsilon}$ vanishes, with the imaginary part following an asymptotic $\imaginary(\hat{K}_{+}(\omega\to\infty))=iK(0)/\omega+\mathcal{O}(\omega^{-3})$ (dotted black line).

While MKT correctly captures the relationship between the work of adhesion and contact-line friction,[blake_influence_2002, de_ruijter_dynamic_1999, duvivier_experimental_2011, nold_hydrodynamic_2024] the picture that contact-line motion is driven by fast absorption processes involving only single molecules is not supported by our findings. While at intermediate times (Regime II in Fig. 2), the contact-line motion is dominated by molecular vibrations, these oscillations are not followed by an exponential decay of correlations at large times (Regime III in Fig. 2). In contrast, we find that at large times the dynamics are governed by coupling to the hydrodynamics of the system with the majority of the contact-line friction emerging due to this coupling. This observation is in conflict with descriptions of dynamic wetting where the movement of the contact line is attributed to rapid processes involving only a few molecules at the contact line, like MKT.

V Conclusion

We have derived a generalized description of contact-line friction, in close analogy to the description of hydrodynamic friction. This approach yields, within the linear response regime, a frequency-dependent description of contact-line friction obtained directly from equilibrium data. We introduce a robust procedure to extract the contact-line coordinate and its correlation functions from MD data, and show that the corresponding correlation functions obey the expected short-time expansion of well-defined variables. The VACF and memory kernel feature damped oscillations caused by low-frequency molecular vibrations, and on large time scales a power law decay arising from coupling to the systems hydrodynamic modes. From the frequency-dependent friction, we infer that friction mostly emerges due to this coupling. As a result the microscopic contact line dynamics are neither Markovian nor determined by the movement of a few molecules at the contact line. On the contrary the dynamics are strongly non-Markovian and influenced by hydrodynamic modes in the bulk of the droplet.

Acknowledgments

The authors thank Vishal Dadich, Viktor Klippenstein, Florian Müller-Plathe, and Burkhard Dünweg for helpful discussions. This project was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) in the framework to the collaborative research center “Multiscale Simulation Methods for Soft Matter Systems” (TRR 146) under Project No. 233630050.

Data Availability Statement

The data that supports the findings of this study, input files for simulations, and the raw data used to generate the figures shown is freely available.

Appendix A MD Simulation Details

All MD simulations are performed with LAMMPS[thompson_lammps_2022] (stable release 2. Apr. 2025). Integration of the liquid particles is performed in the $NVT$ ensemble. The positions of the solid particles are not updated. The center of mass of the liquid is fixed by removing the center of mass momentum of the liquid at every step. Initial configurations are generated using the molecular builder Moltemplate[jewett_moltemplate_2021], and snapshots visualized with VMD.[humphrey_vmd_1996] During production runs, the density field is computed with a custom LAMMPS compute. From the density, we obtain the contact line position as outlined in the main text. Correlation functions are obtained via fast Fourier transforms as implemented in NumPy[harris_array_2020] and then averaged over $10$ production runs and the $4$ contact lines.

For the generic LJ liquid, the mass $m_{\mathrm{LJ}}$, the depth of the potential well $\epsilon_{\mathrm{LL}}$ and the distance $\sigma$ at which the potential is $0$ are set to $1$. Each molecule is modeled as an eight-bead oligomer with adjacent beads connected by finitely extensible non-linear elastic bonds[kremer_dynamics_1990] using a bond constant of $k=35\,\epsilon\,\sigma^{-2}$ and an extensibility of $R_{0}=1.5\,\sigma$ and all interactions are cutoff at $2.5\,\sigma$. The thermal energy $k_{\mathrm{B}}T$ is equal $\epsilon_{\mathrm{LL}}$ and the integration time step is $0.005\,\tau_{LJ}$ for all simulations. Both solid plates are modeled as single-layer primitive cubic sheets with a lattice constant of $0.5\,\sigma$. Liquid-solid interactions are modeled with LJ-parameters of $\epsilon_{\mathrm{LS}}=0.125\,k_{\mathrm{B}}T$ and $\sigma_{\mathrm{LS}}=1$. The liquid phase contains $4,000$ molecules, for a total of $32,000$ liquid beads, and the two solid sheets contain $16,000$ particles each, for a total of $32,000$ solid particles in a simulation box with dimensions $L_{x}=100\,\sigma$, $L_{y}=40\,\sigma$, and $L_{z}=21.001\,\sigma$, with a distance of $21.0\,\sigma$ between the solid plates. The boundaries in the $x$ and $y$ directions are periodic, and a non-periodic boundary in the $z$ direction. To relax the liquid density, we first integrate the system for $2,500\,\tau_{\mathrm{LJ}}$ with a Langevin thermostat[schneider_molecular-dynamics_1978, dunweg_brownian_1991] constant of $0.5\,\tau_{\mathrm{LJ}}$. Then the contact angle is equilibrated for $25,000\,\tau_{\mathrm{LJ}}$ using a Nosé–Hoover-thermostat[tuckerman_liouville-operator_2006] with a thermostat constant of $1\,\tau_{\mathrm{LJ}}$. To obtain $10$ independent configurations, we randomize the velocity of all liquid particles and integrate the system using a Nosé–Hoover thermostat with a thermostat constant of $5\,\tau_{\mathrm{LJ}}$ for $50,000\,\tau_{\mathrm{LJ}}$. We then perform a production run with all $10$ systems using a Nosé–Hoover thermostat with a thermostat constant of $5\,\tau_{\mathrm{LJ}}$ for $50,000\,\tau_{\mathrm{LJ}}$, outputting the density field at every step. The density field is calculated from the bead positions (not the center of mass).

Water molecules are modeled with the SPC/E water model.[berendsen_missing_1987] Water–graphene interactions are described by a Lennard-Jones potential with LJ-parameters $\epsilon_{\mathrm{WC}}=0.4389\,\mathrm{kJ}\,\mathrm{mol}^{-1}$ and $\sigma_{\mathrm{WC}}=0.3367\,\mathrm{nm}$.[carlson_modeling_2024] Bond lengths and angles in water are constrained to their equilibrium values using the SHAKE algorithm with a tolerance of $0.0001$ and a maximum of $10$ iterations.[ryckaert_numerical_1977] The system consists of $14,850$ water molecules and $20,000$ carbon atoms in two graphene sheets with a $CC$ bond length of $0.142\,\mathrm{nm}$. The dimensions of simulation box are $L_{x}=24.595\mathrm{nm}$, $L_{x}=10.65\mathrm{nm}$, and $L_{z}=5.0001\mathrm{nm}$ with a distance of $5.0\,\mathrm{nm}$ between the graphene sheets. The boundary is periodic in the $x$- and $y$-directions and non-periodic in the $z$-direction. Long-range electrostatics are treated with the slab particle–particle particle–mesh method using a real-space cutoff of $1.0\,\mathrm{nm}$, an accuracy of $10^{-4}$, and a slab factor of $3$. We first minimize the energy of the system and then equilibrate for $4\,\mathrm{ns}$ using a Nosé–Hoover thermostat with a thermostat constant of $1\,\mathrm{ps}$. To obtain $10$ independent configurations, we randomize the velocity of all liquid particles and then integrate the system for $2\,\mathrm{ns}$ using a Nosé–Hoover thermostat with a thermostat constant of $1\,\mathrm{ps}$. Finally, we perform production runs with all $10$ systems for $10\,\mathrm{ns}$ using a Nosé–Hoover thermostat with a thermostat constant of $1\,\mathrm{ps}$, outputting the density field every step. The density field is calculated from the position of the oxygen atoms.