2.1. Stochastic thin film equation

By assuming a small aspect ratio of the viscous channel, $\hat {h}/L\ll 1$ , one can apply the lubrication approximation to describe the flow of a Newtonian fluid with viscosity $\mu$ in the channel, leading to a parabolic velocity profile (Batchelor Reference Batchelor2000). A random stress tensor is introduced in the momentum equations to account for the thermal fluctuations in the fluid, which are significant on the relevant nanometric scale (Landau & Lifshitz Reference Landau and Lifshitz1987). By imposing no-slip and kinematic boundary conditions at the membrane, one arrives at the following stochastic thin film (Davidovitch, Moro & Stone Reference Davidovitch, Moro and Stone2005; Mecke & Rauscher Reference Mecke and Rauscher2005; Grün et al. Reference Grün, Mecke and Rauscher2006; Durán-Olivencia et al. Reference Durán-Olivencia, Gvalani, Kalliadasis and Pavliotis2019) that can be used to describe $\hat {h}(\hat {x},\hat {y},\hat {t})$ :

(2.1) \begin{equation} \frac {\partial \hat {h}(\hat {x},\hat {y},\hat {t})}{\partial \hat {t}} = \hat {\boldsymbol{\nabla }}\boldsymbol{\cdot }\left (\frac {\hat {h}^3(\hat {x},\hat {y},\hat {t})}{12\mu }\hat {\boldsymbol{\nabla }} \hat {p}(\hat {x},\hat {y},\hat {t}) \right) + \sqrt {\frac {k_B T}{6\mu }}\hat {\boldsymbol{\nabla }} \boldsymbol{\cdot }\left (\hat {h}^{3/2}(\hat {x},\hat {y},\hat {t})\boldsymbol{\hat {\eta }}(\hat {x},\hat {y},\hat {t})\right) \end{equation}

where $\hat {\boldsymbol{\nabla }}$ represents the two-dimensional (2-D) gradient operator $(\partial /\partial \hat {x},\partial /\partial \hat {y})$ and the dependent variable $\hat h$ represents height in the third spatial direction.

The first term on the right-hand side of (2.1) represents the change in film height due to a lateral flux driven by the horizontal pressure gradient $\hat {\boldsymbol{\nabla }} \hat {p}$ . The second term on the right-hand side incorporates the effect of thermal fluctuations at temperature $T$ , by introducing a random flux in accordance with the fluctuation–dissipation theorem and averaged across the channel (Davidovitch et al. Reference Davidovitch, Moro and Stone2005; Mecke & Rauscher Reference Mecke and Rauscher2005; Grün et al. Reference Grün, Mecke and Rauscher2006). Here, $\boldsymbol{\hat {\eta }} (\hat {x},\hat {y},\hat {t})$ , is a random vector with Gaussian white noise uncorrelated in both time and space, i.e. $\langle \hat {\eta }_i(\hat {x},\hat {y},\hat {t})\rangle =0$ and $\langle \hat {\eta }_i(\hat {x},\hat {y},\hat {t})\hat {\eta }_j(\hat {x}^{\prime},\hat {y}^{\prime},\hat {t}^{\prime})\rangle =\delta _{ij}\delta (\hat {x}-\hat {x}^{\prime})\delta (\hat {y}-\hat {y}^{\prime})\delta (\hat {t}-\hat {t}^{\prime})$ , where $\langle \; \rangle$ is the ensemble average, $\delta _{ij}$ is the Kronecker symbol, $\delta$ is the Dirac distribution and $k_B$ is the Boltzmann constant.

Here, we note that the nonlinear prefactor $\hat {h}^3/(12\mu)$ arises from viscous resistance to the flow, describing the mobility of the fluid. In fact, (2.1) can be re-written as a gradient flow (Durán-Olivencia et al. Reference Durán-Olivencia, Gvalani, Kalliadasis and Pavliotis2019; Cates Reference Cates2022; Sprittles et al. Reference Sprittles, Liu, Lockerby and Grafke2023) of the form

(2.2) \begin{equation} \frac {\partial \hat {h}(\hat {x},\hat {y},\hat {t})}{\partial \hat {t}} = \hat {\boldsymbol{\nabla }}\boldsymbol{\cdot }\left (M(\hat {h})\hat {\boldsymbol{\nabla }} \frac {\delta \hat {F}}{\delta \hat {h}} + \sqrt {2 k_B TM(\hat {h})}\boldsymbol{\hat {\eta }}\right)\!, \end{equation}

where $M(\hat {h})$ is the mobility and $\hat {F}[\hat {h}(\hat {x},\hat {y})]$ is an energy functional that gives rise to the pressure. Depending on the problem and assumptions, the mobility can take other forms (Glasner Reference Glasner2008). For example, in a crowded, narrow section of cytoplasm, one could use Darcy’s law to describe the flow through a porous medium, which reduces the mobility to ${\sim} \hat {h}$ (Gnann & Petrache Reference Gnann and Petrache2018). A slip boundary, however, would enhance mobility by introducing an additional ${\sim} \hat {h}^2$ term (Zhang, Sprittles & Lockerby Reference Zhang, Sprittles and Lockerby2020). The ${\sim} M^{1/2}$ prefactor of the fluctuation term ensures that detailed balance is satisfied (Durán-Olivencia et al. Reference Durán-Olivencia, Gvalani, Kalliadasis and Pavliotis2019; Liu et al. Reference Liu, Sprittles and Grafke2024). The free energy $\hat {F}$ depends on what forces drive the flux in the channel, as described in the following sections.

2.2. Membrane dynamics

Deformation of the membrane, as shown in figure 1 $(a)$ , results in tension and bending forces, which lead to a change in the pressure $\hat p(\hat {x},\hat {y},\hat {t})$ of the fluid in the channel. As the membrane changes shape, so does the fluid flux in accordance with (2.1). The membrane is idealised as an isotropic linearly elastic solid with Young’s modulus $E$ , Poisson ratio $\nu$ and thickness $d$ (Landau & Lifshitz Reference Landau and Lifshitz1986). In the limit of small deflections, i.e. small spatial gradients in $\hat {h}(\hat {x},\hat {y},\hat {t})$ , which is natural from the scale separation, and small strains, the coupling between the pressure in the channel and the deflection of the membrane can be represented by the Föppl–von Kármán (FvK) equations (Audoly & Pomeau Reference Audoly and Pomeau2008; Howell, Kozyreff & Ockendon Reference Howell, Kozyreff and Ockendon2008). The FvK equations incorporate the effects of both bending and tension, the latter of which is highly coupled to the deformations in the membrane. We consider a membrane subject to constant tension $\gamma$ , which could arise from some external stretching of the membrane at its edges (or a pressure inside an adhered cell), which is much greater than the tension caused by local membrane deflections. Under this assumption, the relationship between membrane deflections and channel pressure becomes (Evans Reference Evans1985; Peng & Lister Reference Peng and Lister2019)

(2.3) \begin{equation} \hat {p}_{\textit{elastic}}(\hat {x},\hat {y},\hat {t}) =B \hat {{\nabla }} ^4\hat {h}(\hat {x},\hat {y},\hat {t}) - \gamma \hat {{\nabla }} ^2\hat {h}(\hat {x},\hat {y},\hat {t}), \end{equation}

where $B=Ed^3/12(1-\nu ^2)$ is membrane bending rigidity and $\gamma$ is the membrane tension (which is analogous to the surface tension at a fluid–fluid interface). Both of these components are generally present, but their relative strength will depend on the ratio of the characteristic horizontal length scale of the system, $L$ , to the bendocapillary length $l_{BC}= \sqrt {B/\gamma }$ (Deserno Reference Deserno2015). For $L/l_{BC}\ll 1$ , bending should be the prominent driving force, whereas for $L/l_{BC}\gg 1$ , we expect membrane tension to dominate. Note that the tension in a membrane can also be described by an integral constraint for the membrane length (Kodio, Griffiths & Vella Reference Kodio, Griffiths and Vella2017), or by solving the full FvK to give a more complete description of membrane deformation (Pihler-Puzović et al. Reference Pihler-Puzović, Périllat, Russell, Juel and Heil2013).

The pressure terms in (2.3) can be formulated as a gradient flow of (2.2) as the two contributions can be rewritten as functional derivatives of free energies. The bending energy is $\hat {F}_{{bend}}[\hat {h}(\hat {x},\hat {y})] = \int ({B}/{2}) |\hat \nabla^2 \hat h |^2\, {\rm d}\hat {A}$ and the surface energy is $\hat {F}_{\gamma }[\hat {h}(\hat {x},\hat {y})] = \int ( {\gamma }/{2}) \lvert \hat {\boldsymbol{\nabla }} \hat {h}\rvert ^2 \,{\rm d}\hat {A}$ .

2.3. Dynamics of adhesive molecules

When proteins or other adhesive molecules bind across the gap between the two membranes, additional forces are generated that contribute to the pressure $\hat p(\hat {x},\hat {y},\hat {t})$ (Carlson & Mahadevan Reference Carlson and Mahadevan2015a , Reference Carlson and Mahadevanb ). Here, we use a Hookean elastic spring model to describe how bound proteins resist stretching and compression. Given a concentration $\hat {c}(\hat {x},\hat {y},\hat {t})$ of the bound adhesive molecules, an additional contribution to the pressure $\hat p_{\textit{adh}}$ takes the form

(2.4) \begin{equation} \hat p_{\textit{adh}}(\hat {x},\hat {y},\hat {t})=\kappa (\hat {h}(\hat {x},\hat {y},\hat {t})-l)\hat {c}(\hat {x},\hat {y},\hat {t}), \end{equation}

where $l$ is the equilibrium bond length and $\kappa$ is the spring coefficient.

To determine $\hat {c}(\hat {x},\hat {y},\hat {t})$ , we use a minimal description of the chemical kinetics of a single protein species, a full description of this model is given in Appendix A (Bell, Dembo & Bongrand Reference Bell, Dembo and Bongrand1984; Carlson & Mahadevan Reference Carlson and Mahadevan2015a ). The adhesive molecules are assumed to be uniformly distributed on both surfaces at a constant concentration $c_0$ . The concentration of pairs of unbound proteins is thus $c_0-\hat {c}(\hat {x},\hat {y},\hat {t})$ . The molecules can bind or unbind at a rate of $(c_0-\hat c)K_{{on}}$ or $\hat cK_{\textit{off}}$ , respectively. The rate constant for binding, $K_{{on}}$ , is a Gaussian distribution about the equilibrium length of the bond, $l$ , with a standard deviation $\sigma$ ; whereas the rate constant for unbinding, $K_{\textit{off}}$ , is constant (Bell et al. Reference Bell, Dembo and Bongrand1984; Qi et al. Reference Qi, Groves and Chakraborty2001). Furthermore, we assume that the molecular dynamics of binding/unbinding are much faster than the time scale of viscosity-mediated membrane deformation, so that the adhesive molecule concentrations immediately adjust to the equilibrium values for a specific membrane height. Balancing the binding and unbinding rates then gives us the following expression for the concentration of a single species of bound molecules (Carlson & Mahadevan Reference Carlson and Mahadevan2015b ):

(2.5) \begin{equation} \hat {c}(\hat {x},\hat {y},\hat {t})=c_0\frac {\exp \left (- ((\hat {h}(\hat {x},\hat {y},\hat {t})-l)/\sigma )^2\right)}{\exp \left (-((\hat {h}(\hat {x},\hat {y},\hat {t})-l)/\sigma)^2\right)+\tau _{{on}}/\tau _{\textit{off}}}, \end{equation}

where a value of $\sigma /l=0.2$ is used in the results of this study and the ratio of kinetic times for binding to unbinding is set to $\tau _{on}/\tau _{\textit{off}}=1/3$ Carlson & Mahadevan (Reference Carlson and Mahadevan2015a ,Reference Carlson and Mahadevan b ). We note that with $\hat {c}(\hat {x},\hat {y},\hat {t})$ given by (2.5), the molecule kinetics can be reduced to a single function of $\hat {h}(\hat {x},\hat {y},\hat {t})$ , meaning that we do not need to solve for $\hat {c}(\hat {x},\hat {y},\hat {t})$ as a variable in our system, which would be the case if the time scale for molecules to bind or unbind were of the order of the time scale for viscous fluid transport in the channel. The pressure given by (2.4) and (2.5) can conveniently be written as a free energy $\hat {F}_{\textit{adh}}[h]$ , as is described in Appendix A.

2.4. Non-dimensional analysis

We work with non-dimensional versions of (2.1) in the remainder of this article. Different scalings of the lengths and time are used in accordance with the dominant physical effects for the work presented in the subsequent sections. In each case, (2.1) is left with one non-dimensional parameter, $Q$ , that is directly proportional to $\langle |\delta \hat {h}|\rangle$ , the average amplitude of thermal fluctuations in the film (Dhaliwal et al. Reference Dhaliwal, Pedersen, Kadri, Miquelard-Garnier, Sollogoub, Peixinho, Salez and Carlson2024). Here, we outline the non-dimensionalisation used for each case; further details can be found in Appendix B.

In § 3, we study (2.1) on a one-dimensional (1-D) domain (which represents a 2-D film since $h$ is the height in the $z$ -direction), with $\gamma =0$ in (2.3) and no adhesive molecules/proteins. In this case, the horizontal length $\hat {x}$ is scaled by the domain size $L$ , as it is the only relevant horizontal length scale, the film height $\hat {h}(\hat {x},\hat {t})$ is scaled by the initial film height $\hat {h}_0$ and time $\hat {t}$ is scaled by $\tau _\mu = {12\mu L^6}/{B\hat {h}_0^3}$ , which is the time scale of viscous relaxation of an elastohydrodynamic thin film and can be obtained by scaling the left-hand side of (2.1) with the bending term on the right (Carlson Reference Carlson2018). With these scalings, the dimensionless version of (2.1) is

(2.6) \begin{equation} \frac {\partial h}{\partial t} = \frac {\partial }{\partial x} \left (h^3\frac {\partial }{\partial x} \left (\frac {\partial ^4}{\partial x^4} h\right)\right) + Q_{1\text{-}\textrm{D}}\frac {\partial }{\partial x} \left (h^{3/2}\eta \right)\!, \end{equation}

where $Q_{\textrm{1-D}}= ({L}/{\hat {h}_0})\sqrt { {2k_B T L}/{Bw}}$ . Here, $w$ is the width of the quasi-one-dimensional (quasi-1-D) film in the y-direction.

In § 4, we study (2.1) on a 2-D domain (i.e. a three-dimensional (3-D) film), with either  $\gamma$ or $B$ set to zero in (2.3) and the protein pressure contribution described in § 2.3. In this case, a physical horizontal length scale can be obtained by balancing the relevant membrane pressure with the protein spring pressure (Carlson & Mahadevan Reference Carlson and Mahadevan2015a ). When bending dominates, this gives $L_B= (B/c_0\kappa)^{1/4}$ and while tension dominates, $L_\gamma =(\gamma /c_0\kappa)^{1/2}$ . We scale $\hat {x}$ and $\hat {y}$ by this length scale, $\hat {h}$ by the protein equilibrium length $l$ , the protein concentration $\hat {c}$ by $c_0$ , and $\hat {t}$ by a time scale constructed as for the 1-D case, but with the gradients scaled by the new horizontal length scale $L_\gamma$ or $L_B$ . For the bending-driven case, inserting this scaling, i.e. $\hat {x}=L_Bx$ , $\hat {y}=L_By$ and $\hat {t}=t( {12\mu B^{1/2}}/{l^3(c_0\kappa)^{3/2}})$ , gives the following dimensionless equation:

(2.7) \begin{equation} \text{Bending} \,(\gamma =0): \quad \frac {\partial h}{\partial t} = \boldsymbol{\nabla }\boldsymbol{\cdot }\left (h^3\boldsymbol{\nabla }\left ({\nabla} ^4h+ (h-1)c\right) \right) + Q_{B}\boldsymbol{\nabla }\boldsymbol{\cdot }\left (h^{3/2}\boldsymbol{\eta }\right)\!, \end{equation}

where $Q_{B}= ({1}/{lL_B})\sqrt { {2k_B T }/{c_0\kappa }}$ . For tension-driven films, the equivalent scalings, i.e. $\hat {x}=L_\gamma x$ , $\hat {y}=L_\gamma y$ and $\hat {t}=t({12\mu \gamma }/{l^3(c_0\kappa)^{2}})$ , give us the appropriate form of the non-dimensional thin film equation:

(2.8) \begin{equation} \text{Tension}\, (B=0): \quad \frac {\partial h}{\partial t} = \boldsymbol{\nabla }\boldsymbol{\cdot }\left (h^3\boldsymbol{\nabla }\left (-{\nabla} ^2h+ (h-1)c\right) \right) + Q_{\gamma }\boldsymbol{\nabla }\boldsymbol{\cdot }\left (h^{3/2}\boldsymbol{\eta }\right)\!, \end{equation}

where $Q_{\gamma }=({1}/{lL_{\gamma }})\sqrt { {2k_B T }/{c_0\kappa }}$ .

2.5. Finite element solver

In this article, (2.6), (2.7) and (2.8) are solved using the finite element method on a rectangular domain. In all simulations, periodic boundary conditions are imposed in the horizontal directions at the four boundaries. Since the film height $h$ is a dependent variable, simulations are performed in both 1-D and 2-D domains, which represent 2-D and 3-D films, respectively. The source code for the simulations in this article is available on GitHub (Dhaliwal Reference Dhaliwal2025).

To solve (2.1), we set up a separate equation for the pressure $p$ , with contributions coming from (2.3) and (2.4). Equation (2.1) is divided into a system of two (plus an additional one for the film curvature ${\nabla} ^2 h(x,y,t)$ when there is a bending pressure) coupled second-order partial differential equations. The system is reformulated into weak form where boundary terms can be neglected due to the periodic boundary conditions. The scalar fields ${\nabla} ^2 h(x,y,t)$ , $p(x,y,t)$ and $h(x,y,t)$ are then discretised with linear elements, and the system of equations is solved using a Newton’s method solver from the FEniCS finite element library (Logg et al. Reference Logg2012). Time integration is performed using an implicit first-order finite difference scheme. The domain and its discretisation in space and time are chosen according to the relevant non-dimensional form presented in § 2.4. For the 1-D simulations presented in § 3.2, a domain of length $L=1$ is used with $\Delta x=0.01$ and $\Delta t =1\times 10^{-8}$ . For the 2-D simulations presented in § 4.2, a square $100\times 100$ -cell grid is used with both $\Delta x$ and $\Delta t$ varying in the range $1{-}3$ (which means that $L$ is in the range $100{-}300$ ) for the various results presented.

The random vector $\boldsymbol{\eta }(x,y,t)$ is implemented by choosing random numbers using the ‘normal’ function in the ‘random’ class of NUMPY (Harris et al. Reference Harris2020). For the numerical solution at each time step, the two components of $\boldsymbol{\eta }(x,y,t)$ are each assigned a new value at every point in the mesh. Values are drawn from a Gaussian distribution with zero mean and a variance of $1/(\Delta x \Delta t)$ in 1-D and $1/(\Delta x^2 \Delta t)$ in 2-D. Due to the stochastic nature of the problem, each individual run is not to be considered as deterministic. For each set of parameters, we run multiple independent realisations and report the ensemble averages of the extracted/predicted parameters.

3.1. Rare-event theory

It is well known that thermal fluctuations can generate waves across the membrane, where the average wave amplitudes $\delta h_q$ of frequency $q$ , or ‘roughness’ of the membrane, depends on bending moduli and tension coefficient (Helfrich & Servuss Reference Helfrich and Servuss1984; Rautu et al. Reference Rautu, Orsi, Michele, Rowlands, Cicuta and Turner2017). If the $h_0-h^*$ is of the same order as $\delta h_q$ , then thermal fluctuations can easily initiate binding. If $h_0-h^*$ is sufficiently larger than $\delta h_q$ , attachment of the molecules then requires the membrane profile $h(x,t)$ to attain a rather unlikely shape, which may take a long time. The process can thus be thought of as $h(x,y,t)$ fluctuating in an energy landscape that does not favour large deviations from $h_0$ , until it eventually gets ‘lucky’ and reaches $h^*$ at some point in the domain. Although the waiting time for protein binding is a random variable, the rare event theory can provide a prediction for the ensemble-averaged waiting time, known as the Eyring–Kramers law, which states that $\langle t_B \rangle \sim C \exp ({2(F[h_B(x)]-F[h_0])/Q_{\text{1-D}}^2})$ , where $\langle t_B \rangle$ is the ensemble-averaged binding time, $C$ is a prefactor, $F[h(x,t)]$ is the energy of a profile $h(x,t)$ and $h_B(x)$ is the final binding profile (Eyring Reference Eyring1935; Kramers Reference Kramers1940; Liu et al. Reference Liu, Sprittles and Grafke2024).

If the membrane primarily exhibits tension rather than bending, i.e. $L\gt l_{BC}$ , the energy is determined by the second term in (2.3), and the problem is analogous to the rupture of a thin viscous liquid film with a free surface due to van der Waals forces. This problem was recently studied by Sprittles et al. (Reference Sprittles, Liu, Lockerby and Grafke2023), where the rare-event description based on the Eyring–Kramers formula was validated by both numerical simulations of the stochastic thin film equations and in molecular dynamics simulations. In this paper, these methods are adapted to determining the binding time between two membranes when rather the bending dominates the pressure term in (2.3).

When bending dominates over tension, i.e. $L\lt l_{BC}$ , the non-dimensional energy functional for the 1-D membrane becomes

(3.1) \begin{equation} F[h(x)] = \int ^1_0 \frac {1}{2} \left (\frac {\partial ^2 h}{\partial x ^2}\right)^2 \,{\rm d}x. \end{equation}

The pressure term in (2.3) can be recovered by taking a functional derivative of $F$ with respect to $h(x)$ . Finding the average waiting time for adhesion requires finding the profile $h_B(x)$ that minimises $F[h]$ while also maintaining conservation of mass, which enters as the constraint $\int _0^1(h-h_0)\,{\rm d}x = 0$ , and satisfying the periodic boundary condition. This constrained optimisation problem can be solved using the Euler–Lagrange equation, which gives the fourth-order polynomial profile shown by the dashed blue line in figure 2 $(a)$ . The energy corresponding to this profile is $F_{B}= 360 (h_0-h^*)^2$ . With a sharp asymptotics analysis (see Appendix C for details), the predicted average attachment time is given by

(3.2) \begin{equation} \langle t_B \rangle = \frac {1}{\beta (h_B)}\sqrt {\frac {Q_{\text{1-D}}^2}{(2\pi)^7}} \exp \left (720\left (\frac {h_0-h^*}{Q_{\textrm{1-D}}}\right)^2\right)\!, \end{equation}

where the parameter $\beta$ is given by

(3.3) \begin{equation} \beta (h_B) = (h_0-h^*)(2\pi)^6\left (\frac {1}{2}h_0^3 +\frac {3}{8}h_0(h_0-h^*)^2\right)\!. \end{equation}

We note that the expressions in (3.2) and (3.3) predict the attachment time across the channel without any fitting parameters. They are valid in the limit of small noise, i.e. either large $h_0-h^*$ or small $Q_{\text{1-D}}$ , such that the trajectory of initial attachment must cross $h_B(x)$ in the energy landscape.

Figure 2.

$(a)$ Film profiles at time of attachment obtained from 15 independent solutions (centred around the point of ‘contact’) with parameters $Q_{\text{1-D}}=5$ and $h^*=0.3$ . The dotted black line represents the average of the individual simulations and the dashed blue line represents the theoretical prediction from the Euler–Lagrange equation. $(b)$ Average waiting time for adhesion $\langle t^* \rangle$ as a function of $(h_0-h^*)^2$ for different values of the noise amplitude $Q_{\text{1-D}}$ . The lines represent the predicted value of $\langle t_B \rangle$ from (3.2) with no free parameter. Error bars represent the standard deviation for a set of $N=15$ simulations for each data point. The shaded blue colour is intended as a guide to the eye to emphasise the region where rare-event theory is valid, i.e. attachment events are sufficiently unlikely.

3.2. Numerical predictions

Numerical simulations of (2.6) are used to test the prediction in (3.2) considering only the effect of membrane bending. The simulations are initiated with a flat membrane of height $h_0=1$ . Since there are no adhesive molecules, i.e. no force pulling/pushing on the membrane, the minimum film height fluctuates randomly with time, occasionally reaching a lower value than it has before. The simulations are run until the lowest point in the film reaches the cutoff height $h^*$ , indicating that the molecules are now within binding range for initial attachment, and we denote this time as $t^*$ . The point in the periodic domain at which attachment occurs is random. If the deformation required for film rupture, $h_0-h^*$ , is small, attachment occurs rapidly and the profile shape at $t^*$ varies significantly from simulation to simulation. When $h_0-h^*$ is larger, however, the minimum film height fluctuates for a long time before reaching attachment (the number of time steps being many orders of magnitude), often coming close many times before finally reaching $h^*$ . When this is the case, the individual profiles at $t^*$ are consistent, as shown in figure 2 $(a)$ , which suggests that attachment does indeed tend to occur at a specific minimum energy profile.

When the value of $h_0-h^*$ is large, as is the case for the profiles shown in figure 2 $(a)$ , one would expect that the final profile is similar to the fourth-order polynomial predicted theoretically in § 3.1. The dashed blue line in figure 2 $(a)$ shows that this is indeed the case if the profiles are centred about their minimum and then averaged. This indicates that it is truly rare for fluctuations to make the film reach $h^*$ and that this almost always occurs very close to the minimum energy profile that allows attachment.

The average attachment time, $\langle t^* \rangle$ , is plotted for various values of $Q_{\text{1-D}}$ and $h^*$ in figure 2 $(b)$ . It is clear that the time increases rapidly when $h^*$ is reduced and eventually grows exponentially with $(h_0-h^*)^2$ for constant $Q_{\text{1-D}}$ . The theoretical prediction for the attachment time across a membrane channel from (3.2) provides an excellent quantitative prediction for the average rupture time when $h_0-h^*$ is large. For smaller values of $h_0-h^*$ , the rare-event theory is not valid, as the membrane profile does not necessarily cross the saddle point on its way to attachment. We note that the performance of the theoretical prediction also depends on $Q_{\text{1-D}}$ , as fluctuations are more likely to cause large deformations to the film profile, making the event less rare as well as introducing larger errors in the Laplace method used to approximate the mean first passage time (see Appendix C). Generally, the rare-event prediction is valid when the attachment event is sufficiently rare, meaning that the average attachment time is sufficiently large. This can be achieved either by having a small value of $Q_{\text{1-D}}$ or a large height difference $h_0-h^*$ , as in the blue shaded region of figure 2 $(b)$ .

4.1. Domain coarsening of adhered patches

Equation (2.1) can be seen as the governing equation of a dynamical phase field, in which the film height $h(x,y,t)$ is the sole order parameter (McLeish et al. Reference McLeish, Cates, Higgins, Olmsted and Bray2003; Cates Reference Cates2022). When coupled with the pressure term given by (2.3), it is in fact quite similar to the widely studied dynamical model B, which is used to describe the coarsening of a conserved order parameter (Hohenberg & Halperin Reference Hohenberg and Halperin1977). Nevertheless, a number of features distinguish our system. First, a nonlinear ${\sim} h^3$ mobility (as described in § 2.1) arises from the lubrication flow in the narrow channel. Second, rather than the typical double-well potential, our system imposes an energy landscape arising from the adhesive molecules/proteins described in § 2.3 and Appendix A. Finally, and perhaps most significantly, the classical Cahn–Hilliard Laplacian free energy term (which, in this case, represents isotropic membrane tension) is supplemented by the fourth-order bending term as described in (2.3).

A quantitative understanding of phase separating systems is achieved through the dynamical scaling hypothesis, which states that in the later stage of a coarsening process, the domain structure of a system (as quantified by the structure function $S(\boldsymbol {k},t)=\langle h_{\boldsymbol {k}}(t)h_{\boldsymbol {-k}}(t)\rangle$ ) is constant in time if one rescales the lengths by a single characteristic length scale $L_c(t)$ (Bray Reference Bray1994; Camley & Brown Reference Camley and Brown2011). This concept has been shown to be valid both numerically and experimentally in many systems, with power-law growth displayed when $L_c(t)$ is computed as a moment of $S(\boldsymbol {k},t)$ or the radius of individual particles (Furukawa Reference Furukawa1985; Shinozaki & Oono Reference Shinozaki and Oono1993; Sung et al. Reference Sung, Karim, Douglas and Han1996; Tateno & Tanaka Reference Tateno and Tanaka2021; Saiseau et al. Reference Saiseau, Truong, Guérin, Delabre and Delville2024; Su et al. Reference Su, Ho, Gettel, Rowland, Keating and Parikh2024). For coarsening of the aforementioned model B system, $L_c(t)$ is known to grow according to a ${\sim} t^{1/3}$ power law, as has been validated by a number of numerical studies (Vladimirova, Malagoli & Mauri Reference Vladimirova, Malagoli and Mauri1998; Tiribocchi et al. Reference Tiribocchi, Wittkowski, Marenduzzo and Cates2015), and can also be predicted theoretically using scaling arguments, renormalisation group theory or Lifshitz–Slyozov–Wagner theory (in the limit where the minority phase occupies a small volume fraction) (Lifshitz & Slyozov Reference Lifshitz and Slyozov1961; Wagner Reference Wagner1961; Bray Reference Bray1994; Camley & Brown Reference Camley and Brown2011).

For the case where the bending energy drives coarsening rather than surface tension, it is unclear if the coarsening rate will follow the ${\sim} t^{1/3}$ power law. This essentially corresponds to replacing the surface tension contribution to the classical free energy functional, $F_{\gamma }(h(x,y)) = \int ({1}/{2}) \lvert \boldsymbol{\nabla }h\rvert ^2 \,{\rm d}A$ , by the corresponding bending term $F_{\textit{bend}}(h(x,y)) = \int ({1}/{2}) ({\nabla} ^2 h)^2 \,{\rm d}A$ . We now follow the derivation of Bray (Reference Bray1994) to make a scaling prediction for how bending changes the coarsening dynamics (details are provided in Appendix D.2). We ignore the effects of the nonlinear mobility as well as the specifics of the potential function. In analogy with the coarsening literature, we note that the pressure in our system takes a form that can be interpreted as a chemical potential $\varPhi$ (Hohenberg & Halperin Reference Hohenberg and Halperin1977; Bray Reference Bray1994; Cates Reference Cates2022). We thus adopt the notation $ p\equiv \varPhi =\delta F/\delta h$ , where $F$ consists of the aforementioned bending term and a symmetric double-well potential function $V(h)$ . The chemical potential is thus

(4.1) \begin{equation} \varPhi = \frac {\partial V}{\partial h} +{\nabla} ^4h \end{equation}

and the film height then evolves according to a continuity equation with flux $\boldsymbol{j}=-\boldsymbol{\nabla }\varPhi$ .

In late-stage coarsening when motion of the interface between two domains is slow, diffusion of the order parameter $h$ and chemical potential $\varPhi$ in the bulk is fast, so both should satisfy Laplace’s equation ${\nabla} ^2h={\nabla} ^2\varPhi =0$ in the bulk regions. The flux through the interface, and thus by continuity the velocity of the interface, can then be determined by finding $\varPhi$ at the interface.

At an interface between the two phases with radius of curvature $R$ , $\varPhi$ can be re-expressed in terms of the coordinate $g$ representing distance along the unit vector $\boldsymbol{\hat {g}}$ in the direction perpendicular to the interface ( $g=\pm \infty$ in the bulk and $g=0$ at the centre of the interface). Noting that $\boldsymbol{\nabla }h = (\partial h/\partial g)\boldsymbol{\hat {g}}$ and $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\hat {g}}=1/R$ near the interface, we find that

(4.2) \begin{equation} \varPhi = \frac {\partial V}{\partial h} +\frac {\partial ^4 h}{\partial g^4}+\frac {2}{R}\frac {\partial ^3 h}{\partial g^3}-\frac {1}{R^2}\frac {\partial ^2 h}{\partial g^2}+\frac {1}{R^3}\frac {\partial h}{\partial g}. \end{equation}

The value of $\varPhi$ at the interface is then determined by multiplying (4.2) by $\partial h/\partial g$ and integrating across the interface from $g=-\infty$ to $g=\infty$ . By imposing the boundary condition that $h$ is constant in the bulk and assuming that the chemical potential is a double well with equal potential on both sides, integration gives

(4.3) \begin{equation} \varPhi \Delta h = -\frac {2}{R}\int ^\infty _{-\infty }\left (\frac {\partial ^2 h}{\partial g^2}\right)^2\,{\rm d}g+\frac {1}{R^3}\int ^\infty _{-\infty }\left (\frac {\partial h}{\partial g}\right)^2\,{\rm d}g, \end{equation}

where $\Delta h$ is the height difference between the domains on the inside and outside of the interface, and the integrals $-2\int ^\infty _{-\infty} ({\partial ^2 h}/{\partial g^2})^2\,{\rm d}g$ and $\int ^\infty _{-\infty} ({\partial h}/{\partial g})^2\,{\rm d}g$ can effectively be interpreted as 2-D line tension-like coefficients $\varGamma _1$ and $\varGamma _2$ that represent an energy per unit length associated with the existence of the interface between the two phases in 2-D space. Equation (4.3) thus represents a Gibbs–Thomson-like boundary condition that determines the value of $\varPhi$ on interfaces with radius of curvature $R$ . Since the flux of fluid is proportional to the gradient of the chemical potential (assuming a constant mobility), the velocity of the interface scales as $-\partial \varPhi /\partial g$ . Assuming that the domains are circular with radius $R$ corresponding to the macroscopic length scale $L_c$ , this then gives the following growth law for bending-driven coarsening:

(4.4) \begin{equation} \frac {{\rm d}L_c}{{\rm d}t}\sim \frac {\varGamma _1}{L_c^2}+\frac {\varGamma _2}{L_c^4}. \end{equation}

Equation (4.4) stands in contrast to the model B case which has only a single contribution to the line tension, giving rise to an ${\sim} L_c^{-2}$ growth rate. This suggests that bending-driven coarsening should have a coarsening rate with an exponent somewhere between $1/3$ and $1/5$ with time, with the value being determined by the relative sizes of $\varGamma _1$ and $\varGamma _2$ , which in turn depend on the shape of the potential function. We note that this simple scaling prediction only describes the difference between surface tension-driven and bending-driven coarsening. It does not take into account the other differences between our model and the classical model B system described by Bray (Reference Bray1994) such as the ${\sim} h^3$ mobility and the protein binding potential energy. Since the potential energy function $V(h)$ is not specified, we expect these results to be valid for other adhesive potentials such as van der Waals interactions. Equation (4.4) suggests a fundamental distinction between surface tension-driven and bending-driven coarsening, which will be investigated numerically for our specific system in § 4.2.

4.2. Numerical investigation of bending-driven coarsening dynamics

We now present the results of numerical simulations of the phase separation of adhesive molecules/proteins during membrane adhesion when dominated by bending of the membrane, as described in (2.7). The membrane is always initialised as a flat profile with non-dimensional initial height $h(x,y,0)= \hat {h}_0/l\gt 1$ and the non-dimensional width of the protein binding rate distribution, $\sigma _{{on}}/l=0.2$ . Figure 3( $a$ ) shows the contour plots (see supplementary movie 1 available at https://doi.org/10.1017/jfm.2025.10872) of the film profile at three different times when the initial film height is just at the edge of the range in which the proteins can bind. Fluctuations are clearly present for the first time steps, where the film is pulled to the substrate by the adhesive molecules. Most of the film profile moves towards the equilibrium length of the adhesive molecules ( $h=1$ in non-dimensional terms) to reduce the energy of the system. Due to conservation of liquid mass, excess liquid must flow somewhere else, but periodic boundary conditions prevent fluid from leaving the domain. Thus, in some regions, the film is pushed up, forming circular ‘pockets’ or ‘lumens’ in which proteins are unbound and the film height is significantly larger than the initial height, as can be seen in the first panel to the left in figure 3( $a$ ).

Figure 3.

Contour plots illustrating the height $h(x,y,t)$ of thin films in a 2-D domain under the influence of protein binding and unbinding for different times $t$ , with $h_0=1.4$ . In panel $(a)$ , the fluctuation parameter is $Q_{B}=0.005$ , whereas in panel $(b)$ , it is set to $Q_{B}=0.5$ . Although similar coarsening dynamics occur at late times regardless of $Q_{B}$ , the size of the domains at $t=350$ is somewhat larger for $Q_{B}=0.5$ due to early-time coalescence.

With time, the detached domains coarsen in an Ostwald ripening-like process, resulting in fewer, larger domains as can be seen in the second and third panels of figure 3( $a$ ). Once the circular pockets form, their centres remain essentially fixed as fluid is slowly transported from smaller to larger domains in a diffusion-like manner. When the size of a single pocket drops below a critical threshold, it suddenly ‘implodes’, which causes a minor horizontal rearrangement of nearby domains (see supplementary movies). Figure 4( $a$ ) shows the time evolution of the characteristic size $L_c$ for these domains, as computed using the method of Shinozaki & Oono (Reference Shinozaki and Oono1993), averaging the structure factor $S(\boldsymbol {k},t)$ across 15 independent realisations (details of how we calculte $L_c$ are provided in Appendix D.1). Although the distinct ‘implosion’ events lead to discrete jumps in the growth for an individual simulation, the ensemble-averaged $L_c$ grows smoothly with time following a power law. The growth is significantly slower than the ${\sim} t^{1/3}$ power law observed for a model B system (Bray Reference Bray1994; Tiribocchi et al. Reference Tiribocchi, Wittkowski, Marenduzzo and Cates2015), instead, the exponent is closer to $1/5$ , which can be expected based on our analysis in § 4.1.

Figure 4.

Length scale $L_c$ computed using (D5) from the average of $N=15$ individual simulations under the same conditions as the data in figure 3. In panel $(a)$ , the fluctuation parameter is $Q_{B}=0.005$ , whereas in panel $(b)$ , it is set to $Q_{B}=0.5$ . Power law coarsening with an exponent well below $1/3$ is observed in both cases, but starts off with a larger domain size when $Q_{B}=0.5$ .

The contour plots in figure 3( $b$ ) show the evolution of height profiles of a coarsening film when the amplitude of the thermal fluctuations, represented by the parameter $Q_{B}$ , is increased by two orders of magnitude. With higher noise amplitude, the lumens are irregular in shape and constantly deforming (see supplementary movie 2). The fluctuations also cause significant lateral motion of the domains, which leads to some instances of coalescence at very early times. Despite this, at the later time, large droplets seem to repel each other and domain growth is primarily caused by fluid flow through adhered patches, as was observed in the low noise amplitude case. Interestingly, the power law for domain growth is relatively unaffected by the strength of the fluctuations, as shown in figure 4( $b$ ). The size of the domains before the late-stage power law growth is somewhat higher when the fluctuations are larger, perhaps due to coalescence events at early times, as can be seen in the first panels of figure 3(a,b). The primary effect of fluctuations is thus to provide the perturbations that trigger the instabilities giving rise to the initial configuration, rather than driving coarsening.

To better understand how and why the coarsening rate of our system differs from the more familiar model B system, we perform additional simulations in which we vary the initial membrane height $h_0$ , the fluctuation strength $Q$ , the mobility prefactor in the thin film equation and which term in (2.3) we use. In general, we find that coarsening behaviour is a preserved feature of the system even when these parameters are changed. To ensure an initial configuration with many small pockets in the domain, $h_0$ should be close to the edge of the binding range of the adhesive molecules. If $h_0$ is too small, the entire domain is already in a bound state and few pockets form. If, however, $h_0$ is too large, it takes a very long time for contact to occur, which is typically at only one point in the domain, leading to a small number of initial pockets when coarsening starts.

To begin with, we simulate coarsening using the tension term of (2.3) (corresponding to the well-studied Cahn–Hilliard free energy) instead of the bending term, as in (2.8). As expected, coarsening behaviour is observed, as shown in figure 5. Comparing the first snapshots of figures 5 $(a)$ and 5 $(b)$ demonstrates that larger initial heights lead to fewer domains at early times. At later times, we observe that the coarsening actually appears to be slower when $h_0$ is increased.

Figure 5.

Contour plots illustrating the height $h(x,y,t)$ of thin films in a 2-D domain under the influence of protein binding and unbinding for different times $t$ , with $Q_{\gamma }=0.01$ . In panel $(a)$ , the initial height is $h_0=1.25$ , whereas in panel $(b)$ , it is set to $h_0=1.45$ . Increasing $h_0$ leads to a larger initial domain size, but slower coarsening at late times.

Figure 6 $(a)$ shows how the characteristic length scale $L_c$ of a viscous film profile grows with time when the membrane exhibits only interfacial tension. When $h_0$ is small, the growth exponent is close to ${\sim} t^{1/3}$ , in accordance with what one expects for the well-known model B system. For higher values of $h_0$ , the growth rate for late-stage coarsening decreases. This seems to be caused by the nonlinear ${\sim} h^3$ term arising from the viscous resistance to flow in the governing equation. To confirm this, we perform additional simulations in which we replace the nonlinear mobility by a constant mobility with value  $1$ . The inset in figure 6 $(a)$ shows that $L_c\sim t^{1/3}$ growth is observed regardless of $h_0$ when the mobility is constant, as expected.

Figure 6.

$(a)$ Scaling length $L_c$ as a function of time for tension-driven coarsening in a viscous film with varying $h_0$ when $Q_{\gamma }=0.01$ . $L_c$ is calculated from the average $S(\boldsymbol {k},t)$ for $N=22$ independent simulations. The inset shows the results when a constant mobility is used instead. $(b)$ Scaling length $L_c$ as a function of time for bending-driven coarsening in a viscous film with varying $h_0$ when $Q_B=0.01$ . $L_c$ is calculated from the average $S(\boldsymbol {k},t)$ for $N=20$ independent simulations. The inset shows the results when a constant mobility is used instead.

In figure 6 $(b)$ , similar results are presented, but this time, the tension term in (2.3) is replaced by the bending term, i.e. we set $\gamma \gt 0$ and $B=0$ , and solve 2.8. As suggested by the theory in § 4.1, the power law for bending-driven coarsening is reduced from $L_c\sim t^{1/3}$ . The inset in figure 6 $(b)$ shows that this is indeed the case when the mobility is kept constant, with a growth exponent of approximately $1/4$ consistently observed. For the viscous case, there is a decrease in the growth exponent when the value of $h_0$ is increased, but this effect is not as pronounced as it is for the tension-driven films.

The decreased coarsening rate as $h_0$ is increased may seem counterintuitive at first glance, as a thicker film should lead to less viscous resistance to fluid flow and thus a higher mobility. In the coarsening process, however, the fluid flow between pockets must pass through the adhered region, where the film height is fixed at the equilibrium height of the adhesive molecules ( $h=1$ ), thus restricting the mobility of the flow through this region. Increasing $h_0$ thus only increases the depth of the pockets, which decreases the flux through the interface between the pocket and the attached region.

Figure 7.

$(a)$ Scaling length $L_c$ as a function of time for tension-driven coarsening in a viscous film with varying $Q_{\gamma }$ . These results are for $h_0=1.25$ and the $L_c$ is calculated from the average $S(\boldsymbol {k},t)$ for $N=22$ independent simulations. $(b)$ Scaling length $L_c$ as a function of time for bending-driven coarsening in a viscous film with varying $Q_B$ . These results are for $h_0=1.4$ and the $L_c$ is calculated from the average $S(\boldsymbol {k},t)$ for $N=15$ independent simulations.

Figure 7 shows the coarsening behaviour for both tension-driven and bending-driven viscous films as $Q$ is varied by two orders of magnitude. In both cases, the value of $L_c$ at the beginning of coarsening is larger when the fluctuations are increased to $Q=0.5$ , due to increased coalescence at very early times. In the bending case, large fluctuations also lead to an earlier onset of the power-law coarsening regime, as can be seen in figure 7 $(b)$ . During the late coarsening stage, however, fluctuations do not play a significant role. In figure 7 $(a)$ , we observe that increasing the magnitude of fluctuations leads to a slight decrease in the power law for tension-dominated films, while figure 7 $(b)$ shows that changing $Q_B$ does not seem to significantly affect the power law of bending-dominated films. This seems to suggest that the main role of the fluctuations is simply to provide perturbations at early times which initiate the adhesion process.