# Fluctuation-driven dynamics in nanoscale thin-film flows: physical insights
from numerical investigations

Chengxi Zhao<EMAIL_ADDRESS>Department of Modern Mechanics,
University of Science and Technology of China, Hefei 230026, China  School of
Engineering, University of Warwick, Coventry CV4 7AL, United Kingdom  Jingbang
Liu<EMAIL_ADDRESS>Mathematics Institute, University of Warwick,
Coventry CV4 7AL, United Kingdom  Duncan A. Lockerby<EMAIL_ADDRESS>School of Engineering, University of Warwick, Coventry CV4 7AL, United Kingdom
James E. Sprittles<EMAIL_ADDRESS>Mathematics Institute,
University of Warwick, Coventry CV4 7AL, United Kingdom

###### Abstract

The effects of thermal fluctuations on nanoscale flows are captured by a
numerical scheme that is underpinned by fluctuating hydrodynamics. A
stochastic lubrication equation (SLE) is solved on non-uniform adaptive grids
to study a series of nanoscale thin-film flows. The Fornberg scheme is used
for high-resolution spatial discretisation and a fully-implicit time-marching
scheme is designed for numerical stability. The accuracy of the numerical
method is verified against theoretical results for thermal capillary waves
during the linear stage of their development. The framework is then used to
study the nonlinear behaviour of three bounded thin-film flows: (i) droplet
spreading, where new power laws are derived; (ii) droplet coalescence, where
molecular dynamics results are reproduced by the SLE at a fraction of the
computational cost and it is discovered that thermal fluctuations decelerate
the process, in contrast to previously investigated phenomena; and (iii) thin-
film rupture, where, in the regime considered, disjoining pressure dominates
the final stages of rupture.

††preprint: APS/123-QED

##  I INTRODUCTION

Bounded planar thin-film flows are common in both nature and technology. A
variety of interesting flow behaviours fall into this class, including rupture
[1], dewetting [2], droplet spreading [3] and sessile droplet coalescence [4].
Dynamics of these processes can often be well described by lubrication
equations (LE), derived using a long-wave approximation to the Navier-Stokes
equations [5], which reduces the modelling problem to solving a single partial
differential equation. Motivated by emerging technologies in micro and
nanoscale fluid dynamics [6], thin-film nanoflows have attracted considerable
interest recently and challenged conventional theories, due to the emergence
of new dominant physics at these scales. An important physical factor at the
nanoscale is thermal noise, which has been shown experimentally to influence
interfacial dynamics through the observation of thermal capillary waves on
interfaces of (ultra-low surface-tension) colloid-polymer mixtures [7].
Recently, by the use of molecular dynamics (MD) simulations, it has been
discovered that fluctuation-driven nanowaves dominate a range of nanoscale
interfacial-flow phenomena, such as nanothread breakup [8], nanojet instablity
[9], nanodroplet coalescence [10] and development of rough interfaces on
nanoscale thin films [11]. Notably, the observations of nanothread breakup in
MD were further confirmed by an analytical model [12] and experiments with
colloid-polymer mixtures [13, 14].

To model thin-film flows with thermal noise mathematically, Grün et al. [15]
derived a stochastic lubrication equation (SLE) by applying a long-wave
approximation to the Landau-Lifshitz-Navier-Stokes equations (the fluctuating
hydrodynamics equations) [16]. Linear instability analysis was then conducted
on the SLE by Mecke & Rauscher [17] to obtain an evolving spectrum for thermal
capillary waves of the film interface. It was shown that thermal noise changes
the spectrum of thermal capillary waves from an exponential decay to a power
law for large wavenumbers, which was then confirmed by the experimental
observations on the dewetting of polymer films [18]. The same behaviour of
spectrum was also confirmed in MD simulations, first done by Willis & Freund
[19] in 2010 and recently, followed by Zhang et al. [20], who showed that
interfacial roughening falls into a universality class [11]. However, theories
based on small deviations from equilibrium cannot, inevitably, describe
inherently nonlinear events such as film breakup.

Besides the theoretical analysis, numerical studies of the fully nonlinear SLE
have also been carried out with uncorrelated noise [21] and a spatially
correlated noise model [22], whose solutions at the linear stage are confirmed
by the theory for thermal capillary waves and experimental data, respectively.
The nonlinear dynamics was also investigated in these numerical studies with
particular attention paid to thin-film rupture time [22] and droplet size
distribution after dewetting [21]. However, analytic results for these
nonlinear behaviours are non-trivial. The only successful attempt is on nano-
droplet spreading [23, 24], where a similarity solution was derived for
‘fluctuation-dominated spreading’ and then verified numerically. Therefore,
numerical investigations have become an important approach to better
understand the nonlinear behaviours of thin-film flows. Interestingly, the
spatial discretisation in previous works [21, 22, 25] was developed only on
uniform grids, meaning that the locally nonlinear behaviours (e.g. multiscale
rupture dynamics) could not be easily resolved with high accuracy, while, in
the deterministic cases (LE), these local behaviours have been captured
accurately on non-uniform grids in a range of different numerical schemes [26,
27, 28, 29].

In this work, we develop an accurate and efficient numerical scheme (on non-
uniform grids) to overcome the aforementioned drawbacks and employ it to
investigate different kinds of bounded thin-film flows; in particular for
flows where we expect locally nonlinear dynamics (such as in rupture and
coalescence).

The article is laid out as follows. In Sec .II, the SLE is introduced, the
numerical scheme is proposed and the correlated-noise model is presented. In
Sec. III, numerical solutions for the SLE are verified against known
analytical results (Sec. III). In Sec. IV, we use the SLE solver to study
three different kinds of thin-film flows with nonlinear dynamics: (i) droplet
spreading on a precursor film verified by new theories developed (Sec. IV.1);
(ii) bounded droplet coalescence, validated by MD (Sec. IV.2) and (iii) thin-
film rupture (Sec. IV.3).

##  II Mathematical Modelling and Numerical Approach

In this section, we introduce the mathematical model which describes nanoscale
thin-film flows and then propose a numerical framework for solving the
associated system of equations. We first present the non-dimensionalised SLE.

###  II.1 Stochastic lubrication equations (SLE)

The SLE for a two-dimensional bounded film was first derived by Grün et al.
[15] and Davidovitch et al. [23], who applied a long-wave approximation to the
fluctuating hydrodynamics equations to describe the dynamics of the interface
by the film height $h(x,t)$ (see Figure 1). Recently, Zhang et al. [30]
proposed a more general SLE using a slip boundary condition, which is crucial
at the nanoscale, so this approach is employed here.

Figure 1: Schematic of a thin film with a perturbed interface

To identify the governing dimensionless parameters, we non-dimensionalise the
SLE with the characteristic scales shown below:

$h=\tilde{h}/h_{0},\quad t=\tilde{t}/(3\mu
h_{0}/\gamma),\quad\Pi=\tilde{\Pi}/(\gamma/h_{0}),\quad\mathcal{N}=\mathcal{\tilde{N}}\sqrt{3\mu
h_{0}^{2}/\gamma},$

where $\tilde{h}$, $\tilde{t}$, $\tilde{\Pi}$ and $\tilde{\mathcal{N}}$
represent the dimensional interface height, time, disjoining pressure and
normally distributed random variable (the model for thermal fluctuations),
respectively. Note that the dimensional material parameters are not given
tildes. $h_{0}$ is the characteristic film height, $\mu$ is the liquid’s
dynamic viscosity and $\gamma$ is surface tension. Following Zhang et al.[30],
the dimensionless SLE is given by:

$\frac{\partial h}{\partial t}=-\frac{\partial}{\partial
x}\left[M(h)\left(\frac{\partial^{3}h}{\partial
x^{3}}-\frac{\partial\Pi}{\partial x}\right)-\sqrt{2\varphi
M(h)}\,\mathcal{N}\right]\,.$ (1)

Here, the mobility $M(h)=h^{3}+3\,l_{\mathrm{s}}h^{2}$, so that the slip
length $l_{\mathrm{s}}=0$ the ‘orginal SLE’ from Ref. [15] is recovered. The
disjoining pressure $\Pi=A/(6\pi h^{3})$, where $A$ is the dimensionless
Hamaker constant, reflecting the strength of the van der Waals forces between
liquid and a substrate [20]. The noise term $\mathcal{N}$ has zero mean and
covariance
$\langle\mathcal{N}(x,t)\mathcal{N}(x^{\prime},t^{\prime})\rangle=\delta(x-x^{\prime})\delta(t-t^{\prime})$.
The dimensionless parameter $\varphi=l^{2}_{\mathrm{T}}/(Wh_{0})$ relates to
the intensity of interface fluctuations, where
$l_{\mathrm{T}}=\sqrt{k_{\mathrm{B}}T/\gamma}$ is the characteristic thermal
fluctuation length and $W$ is the initial thickness of the film
($z$-direction). When $\varphi=0$, the deterministic LE [5] is recovered. In
this work, periodic boundary conditions are considered in all cases and the
initial conditions depend on the type of thin-film flow we simulate. The
remaining parts of this section are concerned with the numerical scheme for
solving Equation (1).

###  II.2 Spatial discretisation on non-uniform grids

To resolve locally large gradients in flow variables, without slowing down the
computation dramatically, we use a non-uniform grid for spatial
discretisation. We choose the well-known scheme proposed by Fornberg [31] to
approximate spatial derivatives using finite differences,

$\frac{d^{k}}{dx^{k}}h(x)\approx\sum_{j=0}^{n}\mathcal{L}_{j,n}^{k}(x)h(x_{j}),$
(2)

where $\\{x_{j}\\}^{s}_{j=0}$ are the grid points chosen for approximation and
$\mathcal{L}_{j,s}^{k}(x)$ are the Fornberg coefficients calculated at $x$.
Details of the calculation and error estimate of the Fornberg coefficients can
be found in Appendix A.

To apply the Fornberg scheme for spatial discretisation we first rearrange the
SLE into a conservative form:

$\frac{\partial h}{\partial t}=-\frac{\partial F(h)}{\partial
x},\quad\mathrm{where}\,\,F(h)=M(h)\left(\frac{\partial^{3}h}{\partial
x^{3}}+\frac{A}{2\pi h^{4}}\frac{\partial h}{\partial x}\right)-\sqrt{2\varphi
M(h)}\,\mathcal{N}.$ (3)

The film height $h(x,t)$ is to be solved at the $m$-th step, $t^{m}$,
according to a certain spatial discretisation $\\{x_{i}\\}_{i=0}^{n}$ to give
an array $\\{h_{i}^{m}\\}_{i=0}^{n}$. We use 3 ($n=2$) and 5 ($n=4$) points to
approximate the first-order and the third-order derivatives, respectively,
with the target point located at the centre to provide second-order accuracy.
The expression for the spatial discretisation of equation (3) at the $i$-th
node is then written as:

$\partial_{t}h_{i}=-\sum^{2}_{j=0}\mathcal{L}^{1}_{j,2}F_{i-1+j},$ (4)

where

$\displaystyle
F_{i}=M(h_{i})\left(\sum^{4}_{j=0}\mathcal{L}^{3}_{j,4}h_{i-2+j}+\frac{A}{2\pi
h_{i}^{4}}\sum^{2}_{j=0}\mathcal{L}^{1}_{j,2}h_{i-1+j}\right)-\sqrt{2\varphi
M(h_{i})}N_{i}\,.$ (5)

Here, $\mathcal{L}^{1}_{j,2}(x_{i})$ using $\\{x_{i-1},x_{i},x_{i+1}\\}$ and
$\mathcal{L}^{3}_{j,4}(x_{i})$ using
$\\{x_{i-2},x_{i-1},x_{i},x_{i+1},x_{i+2}\\}$ are the Fornberg coefficients
for the first and third derivative, respectively.

###  II.3 Implicit time-marching method

Equation (4) is of the form

$\frac{\partial\mathbf{h}}{\partial t}=\mathbf{D}\,\mathbf{h},$ (6)

where $\mathbf{h}$ is the solution vector and $\mathbf{D}(\mathbf{h})$ is a
matrix representing the nonlinear differential operator and has size
$(n+1)\times(n+1)$. Employing the implicit Euler time marching method [32] to
Equation (6) gives

$\frac{\mathbf{h}^{m+1}-\mathbf{h}^{m}}{\triangle
t}=\mathbf{D}(\mathbf{h}^{m+1})\,\mathbf{h}^{m+1},$ (7)

where the superscript $m$ denotes the time-step level. This is equivalent to

$\mathbf{G}=\left[\mathbf{I}-\triangle
t\,\mathbf{D}(\mathbf{h}^{m+1})\right]\mathbf{h}^{m+1}-\mathbf{h}^{m}=\mathbf{0}\,,$
(8)

and this root-finding problem, $\mathbf{G}(\mathbf{h}^{m+1})=\mathbf{0}$, can
be solved by the Newton-Kantorovich method [33], using an initial guess of the
solution $\mathbf{h}^{g}$ (the superscript $g$ denotes ‘guess’) obtained by an
explicit Euler time marching method:

$\mathbf{h}^{g}=\mathbf{h}^{m}+\triangle
t\,\mathbf{D}(\mathbf{h}^{m})\,\mathbf{h}^{m}.$ (9)

If $\mathbf{h}^{g}$ is a distance $\mathbf{q}$ away from the solution, i.e.,
$\mathbf{h}^{m+1}=\mathbf{h}^{g}+\mathbf{q}$, so that
$\mathbf{G}(\mathbf{h}^{g}+\mathbf{q})=\mathbf{0}$, then in the linear
approximation one has

$\mathbf{G}(\mathbf{h}^{g}+\mathbf{q})\approx\mathbf{G}(\mathbf{h}^{g})+\frac{\partial\mathbf{G}(\mathbf{h}^{g})}{\partial\mathbf{h}^{g}}\mathbf{q}=\mathbf{0}\,.$
(10)

From Equation (8) we have

$\displaystyle\frac{\partial\mathbf{G}(\mathbf{h}^{g})}{\partial\mathbf{h}^{g}}=\mathbf{I}-\triangle
t\mathbf{J}\,,$ (11)

where the $\mathbf{J}$ is the Jacobian of $\mathbf{D}(\mathbf{h}^{g})$.
Combining Equation (10) and (11) yields a linear approximation of
$\mathbf{q}$,

$\mathbf{q}\approx-(\mathbf{I}-\triangle
t\mathbf{J})^{-1}\mathbf{G}(\mathbf{h}^{g}),$ (12)

so that we can update $\mathbf{h}^{g}$ as $\mathbf{h}^{g}+\mathbf{q}$ and
repeat the process. If the iteration converges, both $\mathbf{q}$ and
$\mathbf{G}(\mathbf{h}^{g})$ will decrease quadratically, eventually to zero.
The final solution, $\mathbf{h}^{m+1}=\mathbf{h}^{g}$.

To control numerical errors caused by the temporal discretisation, we
implement the following criteria:

  * •

The new solution is not negative at any point, i.e., $h_{i}>0$;

  * •

The maximum value of the time derivative error has to be smaller than a
prescribed upper limit, where the error can be computed using the second-order
derivative in the Taylor expansion; i.e., $\max\limits_{1\leq i\leq
n}\left(\left|\frac{(\triangle
t^{m})^{2}}{h^{m}_{i}}\frac{d^{2}h^{m}_{i}}{dt^{2}}\right|\right)<10^{-3}$;

  * •

The iteration process converges, i.e., specifically, $|q_{i}|/h_{i}$ decreases
monotonically to below $10^{-4}$;

  * •

The number of iteration steps is smaller than 100.

If any condition is not satisfied, we go back to the initial value
($h^{g}_{i}=h^{m}_{i}$) and restart the iteration process with a reduced time
step. If the new $\triangle t$ is smaller than $10^{-16}$, the iteration halts
and the code ends in a failed state. Unavoidably, the implicit Euler Scheme
applied to stochastic differential equations converges slower with reducing
time steps than the equivalent deterministic system [32], forcing a smaller
time step for the SLE (compared to the LE) that significantly increases the
computational costs.

###  II.4 Correlated-noise model

The covariance of uncorrelated fluctuations are described by a Dirac delta
function. The delta function could be approximated by a 2D rectangular
(boxcar) function (in $t$ and $x$) that is non-zero over a time step
($\triangle t$) and grid spacing ($\triangle x$), i.e.,
$\mathcal{N}(x,t)\approx N_{i}^{t}/\sqrt{\triangle t\triangle x}$. Here,
$N_{i}^{t}$ represents computer-generated random numbers, that follows a
normal distribution with zero mean and unit variance. However, this model has
been shown prone to numerical instability for the SLE [34]; problems that are
exacerbated as $\triangle x$ and $\triangle t$ become smaller and the
amplitude of noise becomes larger.

To achieve a robust numerical scheme with non-uniform spatial and temporal
discretisation, we combine the spatially correlated noise model from Ref. [15]
and the temporal one from Ref. [34], so that the noise becomes correlated
beneath the spatial correlation length $L_{\mathrm{c}}$ and the time
correlation length $T_{\mathrm{c}}$. Uncorrelated behaviour can then be
approximated by taking the limit of these lengths to zero, ensuring they are
numerically well resolved throughout the limiting process.

Following [15], the stochastic term, $\mathcal{N}(x,t)$ is expanded using
separation of variables in the Q-Wiener ($\mathcal{W}(x,t)$) process as,

$\mathcal{N}(x,t)=\frac{\partial\mathcal{W}(x,t)}{\partial
t}=\sum^{q\rightarrow+\infty}_{q\rightarrow-\infty}\chi_{q}\,\dot{c}_{q}(t)\,g_{q}(x)\,.$
(13)

Here, the constant $\chi_{q}$ are the eigenvalues of the correlation function
$F_{\mathrm{cor}}$,

$\chi_{q}=\int_{-L/2}^{L/2}F_{\mathrm{cor}}(x)e^{-i\,2\pi qz/L}dx\,$ (14)

where $q$ represents an integer sequence. The expressions for
$F_{\mathrm{cor}}$ and details about this model can be found in Appendix B.
The coefficient $\dot{c}_{q}(t)$ represents a temporally correlated noise
process in our setup, in contrast to [15] where uncorrelated noise is
considered, modelled by a simple linear interpolation between uncorrelated
random noise at the endpoints of the temporal correlation interval [34] (see
Figure 17). As the full implicit time marching method usually provides stable
numerical performance with large time steps ($\triangle t>T_{\mathrm{c}}$),
this temporally correlated model will be only activated, and is only
practically necessary, when capturing local dynamics near singularities for
very small $\triangle t$, e.g., for the final stage of film rupture (see more
details in Appendix B).

###  II.5 Grid-size convergence

To test the grid-size convergence of the numerical approach above, we consider
the simulation of a relatively short film $L=15$ with an initial perturbation
$h(x,0)=-0.7\,\mathrm{cos}(2\,\pi\,x/L)$, which relaxes to a relatively flat
film in the absence of disjoining pressure ($A=0$). We set
$L_{\mathrm{c}}=1.5$ and run the simulations with an increasingly fine grid
spacing for different values of the noise strength $\varphi$, whose influence
is shown in Figure 2.

Figure 2: Interface profiles at two time instants ($t_{1}=20$ and $t_{2}=70$)
from numerical solutions of the SLE with different values of the noise
strength $\varphi$. Note that the same random numbers are employed in all the
cases. Figure 3: (a) Interface profiles obtained with different grids at three
time instances: $t_{1}=20$, $t_{2}=40$ and $t_{3}=70$. Here,
$\varphi=10^{-2}$. (b) Convergence characteristics for decreasing grid size,
where Dev. is average deviation of interface profiles to the finest resolution
profile in (a). Figure 4:  (a) Ensemble-averaged interface profiles obtained
from 200 independent realisations at three time instances: $t_{1}=20$,
$t_{2}=50$ and $t_{3}=80$. Here, $\varphi=10^{-2}$. (b) Convergence
characteristics for decreasing grid size, where Dev. is average deviation of
ensemble-averaged interface profiles to the finest resolution profile in (a).

Notably, besides numerical errors caused by discretisations in the simulation,
statistical errors are also introduced by the stochastic term in the SLE. To
better demonstrate the convergence, we measure two different kinds of errors.

The first one comes from the numerical discretisation scheme. To exclude the
influence of statistical variability, we perform simulations with the same
$\mathcal{N}(x,t)$ for all cases, which is generated before the simulations by
fixing $N_{q}^{t}$ (see Appendix B) at each time instant. Figure 3 (a) shows
simulation results with $\varphi=10^{-2}$ at three time instances for varying
grid size. These interface profiles in Figure 3 (a) agree well with each
other, indicating that the numerical errors indeed converge. We check the
average deviation of each interface profile to that with the finest resolution
calculated (i.e., $\triangle x=0.0094$) shown in Figure 3 (b), where the
second-order convergence, as expected, is confirmed, with the prefactor
depending, unsurprisingly, on noise strength.

Since convergence studies for noisy systems are naturally stochastic, the
second target quantity to test convergence is ensemble-averaged; here the
ensemble consists of 200 independent simulations. The ensemble-averaged
profiles at three-time instances are plotted in Figure 4 (a) for varying grid
size. Here, all profiles are presented relative to the minimum point (i.e., we
plot $h$ against $x-x_{\mathrm{min}}$). Then we check the average deviation of
each ensemble-averaged profile to that with the finest resolution (i.e.,
$\triangle x=0.0094$); see Figure 4 (b), where the convergence is confirmed.
Note the deterministic result (black lines) follow second-order convergence,
while, as also expected, the stochastic cases converge with approximately
first-order convergence.

##  III Numerical verification: thermal capillary waves

In this section, the numerical scheme is verified by comparison to analytical
results known for thermal capillary waves, with particular attention paid to
the effect of the non-uniform grid, which are known to create spurious effects
in some stochastic partial differential equation systems [35].

The fluctuations of the interface in thermal equilibrium [17, 11] can be
described by the classical theory for thermal capillary waves [7, 18] with the
static spectra of each surface mode (wavenumber) derived from the
equipartition theorem. To capture the development of thermal capillary waves
towards the equilibrium state, more advanced theories that include the time
evolution of wave spectra were derived by applying linear stability analysis
to the SLE [22, 20]; we will use the model from Ref. [20] to validate our
numerical solutions. Starting with a smooth initial surface and neglecting
disjoining pressure ($A=0$), in dimensionless form the time evolution of the
spectrum is given by (see Appendix C for derivation):

$|H|_{\mathrm{rms}}=\sqrt{-\frac{\varphi
L}{k^{2}}\left(e^{-2k^{4}t}-1\right)}\,,$ (15)

where $|H|_{\mathrm{rms}}$ represents the root-mean-square of wave
perturbations in the frequency domain, $k$ is the dimensionless wave number
and $L$ is the dimensionless length of the film.

Figure 5: Interface profiles of a bounded film at three time instants, i.e.,
$t_{1}=0.0$ (purple dashed line), $t_{2}=0.001$ (black solid line),
$t_{3}=0.1$ (blue dotted line) and $t_{4}=1.0$ (red dash-dotted line). The
parameters chosen are $\varphi=10^{-3}$ and $L=100$. The inset shows how the
waves develop on the non-uniform grid (with exponential distribution and
markers representing nodal positions), with a higher density as one goes from
left to right. Figure 6: The root-mean-square of (wave) disturbance amplitude
versus the wavenumber at three time instants, i.e., $t_{1}=0.2$, $t_{2}=0.4$
and $t_{3}=0.6$; a comparison of ensemble-averaged SLE simulations (dashed
lines with markers) and the analytical result (solid lines).

Figure 5 shows a typical simulation result of the development of thermal
capillary waves, where perturbations, driven by thermal fluctuations, grow
against time to generate significant capillary waves at the later stage (see
the dash-dotted red line in Figure 5). Note that, to test our non-uniform grid
implementation, the grid nodes are non-uniformly distributed with the largest
grid size $\triangle x_{\mathrm{max}}=0.1$ at $x=0$ and the smallest grid size
$\triangle x_{\mathrm{min}}=0.001$ at $x=L$. We use the spatial correlated
noise with $L_{\mathrm{c}}=0.1\geq\triangle x$. The noise is uncorrelated in
time with $\triangle t=10^{-3}$.

To gather statistics, 50 independent simulations (or ‘realisations’) are
performed. The solution is divided into two regions: $x<L/2$ with coarse grids
and $x>L/2$ with dense ones, to check the influence of the different grid
sizes. On each region, a discrete Fourier transform of $h(x,t)$ is applied to
get the power spectral density. We then ensemble average the power spectral
density at each time instant over the realisations and take the square root to
produce the numerical results in Figure 6 (dash lines with markers). Good
agreement with the theory for thermal capillary waves (Equation (15)) can be
found for both regions (see red and blue line), validating the accuracy of our
numerical scheme on extreme non-uniform grids.

##  IV Numerical Simulations for Physical Insight

Having verified our numerical scheme, in this section we use the SLE solver
established to explore three bounded thin-film flows with strong nonlinear
dynamics: (i) droplet spreading on a precursor film; (ii) sessile droplet
coalescence and (iii) thin-film rupture, all of which are ubiquitous phenomena
observed both in nature and many industrial processes [36, 37].

###  IV.1 Droplet spreading on a precursor film

Figure 7: Spreading profiles at three time instants, i.e., $t_{1}=10$ (black
lines), $t_{2}=10^{2}$ (blue lines) and $t_{3}=10^{3}$ (red lines) with the
no-slip boundary. where the solid lines represent one selected realisation.
The dash-dotted lines are the average from 50 realisations. The dotted lines
are the (LE) solution of the deterministic equation. The inset shows the local
behaviours of the selected realisation (solid lines in the full plot). The
lines symbols represent the adaptive grid nodes.

We now consider the nonlinear dynamics of a droplet spreading on a precursor
film, for which the deterministic scaling is given by Tanner’s law [3], and
the more recent work of Davidovitch [23] resulted in a law that accounted for
thermal fluctuations. In particular, these are given by

$\displaystyle\left\\{\begin{array}[]{ll}\ell\sim
t\,^{1/7},\quad\textrm{Hydrodynamic (deterministic) spreading (Tanner's law)
\,\,\cite[cite]{[\@@bibref{Number}{tanner1979spreading}{}{}]}},\\\ \ell\sim
t\,^{1/4},\quad\textrm{Spreading driven by thermal fluctuations (Fluctuation-
enhanced Tanner's law)
\cite[cite]{[\@@bibref{Number}{davidovitch2005spreading}{}{}]}}\,\,.\end{array}\right.$
(18)

Here, $\ell$ represents a characteristic lateral scale, which is estimated
using the average second moment of $h$ (see Equation (53) in Appendix D )
[23]. The larger power law ($1/4$) in Davidovitch’s theory demonstrates that
the nanoscale spreading is enhanced by thermal fluctuations. However, these
two theories are only valid for spreading with the no-slip boundary condition
($l_{\mathrm{s}}=0$), which is often inaccurate at the nanoscale. Therefore,
we use the SLE with the slip boundary condition ($l_{\mathrm{s}}\neq 0$) to
derive ‘slip-modified’ power laws for the spreading (see Appendix D for the
detailed derivation):

$\displaystyle\left\\{\begin{array}[]{ll}\ell\sim t\,^{1/6},\quad\textrm{Slip-
modified Tanner's law,}\\\ \ell\sim t\,^{1/3},\quad\textrm{Slip-modified
fluctuation-enhanced Tanner's law}.\end{array}\right.$ (21)

Note that the new power laws are larger than those based on the no-slip
boundary condition, predicting that the spreading is accelerated by slip, as
we might expect.

The initial droplet profile for the numerical solutions is given by a section
of a sinusoidal function with a precursor film set over the whole simulation
domain. This approach has been widely used [22, 38, 23, 24] and designed not
only for numerical convenience, but also to circumvent the contact-line
dynamics (which are fascinating, but not the focus of this article). To avoid
rupture of the thin precursor film (the thickness of precursor film is less
than one percent of the initial height of the droplet), we follow Ref. [23]
and switch off the thermal fluctuations on it. Notably, it is unclear whether
this ‘artificial limiter’ modifies the numerical path in the stochastic
process. More advanced approaches, like a Brownian bridge technique [39] or
the inclusion of disjoining pressures, are certainly worthy of future
investigation in this regard, but given the slightly artificial nature of the
precursor film here they are considered beyond the scope of this article. When
the droplet spreads to a precursor-film node and ‘pulls it up’, namely,
$h_{i}>h^{*}$, the fluctuations on this node are activated. Moreover, since we
are only interested in the spreading driven by surface tension and
fluctuations in this work, the disjoining pressure is neglected ($A=0$).

To accurately capture the spreading power-laws, one needs to simulate the drop
spreading for several decades of length. Thus it is computationally expensive
to perform simulations on grids uniformly distributed in both the ‘drop’ and
‘precursor film’, in which the height of many nodes will vary very little over
each time step. To reduce computational costs, we utilise our non-uniform
adaptive grids with $\triangle x_{\mathrm{min}}=0.01$ on the ‘drop’ and
$\triangle x_{\mathrm{max}}=0.1$ on the ‘film’, with the nodes in front of the
advancing contact angles refined automatically ($\triangle x_{i}=\triangle
x_{\mathrm{min}}$) to capture the spreading dynamics (see the inset in Figure
7). Here, each SLE realisation on the uniform grids need about 4 core hours,
while one simulation on the non-uniform adaptive grids costs less than 0.2
core hours; i.e., the speed-up is about 20.

Figure 8: Influence of the slip on the characteristic lateral scale, where
black lines and red lines represent no-slip and slip results respectively. The
slip length $l_{\mathrm{s}}=1.0$. The dash-dotted lines are the solutions of
the deterministic cases. The solid lines are the solutions of the stochastic
cases with $\varphi=10^{-3}$. Dashed lines represent similarity solutions in
Equation (18) and (21).

Figure 7 shows the droplet profiles at different time steps, where the
stochastic profiles are the average of the 50 independent realisations with
$\varphi=10^{-3}$. Note that the initial spreading is much ‘faster’ than that
at the later stage ( $t_{3}-t_{2}$ is much larger than $t_{2}-t_{1}$) due to
the initially stronger capillary forces (from the larger curvatures near where
the drop meets the precursor), while at the later stage the thermal
fluctuations play a significant role and accelerate the process (see the
comparison between dotted lines and dash-dotted lines in Figure 7).

The characteristic lateral scales (i.e., the drop’s width) are plotted in
Figure 8, where the numerical solutions match not only previous analytical
solutions with the no-slip boundary condition (black lines) [23], but also our
new ‘slip-modified’ power laws (red lines) very well, giving us further
confidence that our numerical scheme captures nanoscale flow physics both
accurately and efficiently. Notably, in both boundary conditions, there exists
a transfer from hydrodynamic spreading (Tanner’s law) to the fluctuation-
driven spreading (Davidovitch’s law), showing that the noise dominates over
deterministic relaxation only at the later stages of the spreading ($t\gg 1$).

###  IV.2 Bounded droplet coalescence

Previous studies for droplet coalescence on a substrate have been carried out
in the thin-film regime [4, 40, 41] without any influence of thermal
fluctuations taken into account (i.e., using the LE), where surface tension is
considered as the main driving force. However, Perumanath et al [10] have
shown that the fluctuations are crucial to the dynamics of coalescence of two
‘free’ nano-droplets in a vacuum with MD. A similar influence of noise can be
expected in the coalescence of two ‘bounded’ nano-droplets, which will be
explored by the SLE solver in this section. To verify our numerical
predictions, we perform independent MD realisations for both symmetric-drop
and asymmetric-drop coalescence (see details of MD settings in Appendix E).
Note that we focus on 2D cases here.

Figure 9: MD results for the coalescence. (a) Symmetric-drop coalescence with
the same initial droplet radius (at $t=t_{0}-16.7$), $R=10\,$nm (b)
Asymmetric-drop coalescence with different initial droplet radii (at
$t=t_{0}-12.2$), i.e., $R_{1}=15\,$nm and $R_{2}=10\,$nm.

The MD results of two coalescence cases are presented in Figure 9 with two
separate droplets on the substrate set as the initial conditions. The initial
distance between the two droplets is $50\,$nm. Because of the fully wettable
substrate, both droplets spread first until their contact lines meet. The
moment when the two droplets first connect and form a ‘molecular bridge’ is
defined as $t_{0}$, shown in Figure 9. After that, the liquid bridge grows and
the two droplets eventually merge into one.

Figure 10: Asymmetric-drop coalescence profiles predicted in a representative
simulation by the LE (dash-dotted lines) and SLE (solid lines) at three time
instants: $t_{1}=0.8$, $t_{2}=6.4$ and $t_{3}=13.8$. The inset shows the grid
nodes distributions.

Similar to Sec. IV.1, here the disjoining pressure is also neglected ($A=0$).
The dimensionless fluctuation intensity ($\varphi=5.1\times 10^{-3}$) is
calculated from the liquid transport properties of MD (see the details in
Appendix E). To set the initial configuration for the SLE, $h(x,0)$, the MD
interface profile of each realisation at $t_{0}$ is extracted from uniform
bins along the $x$-axis based on a threshold density. Then we shift all the
coalescence points to the same position and use the averaged interface
profiles as $h(x,0)$ for the SLE. Simulations are carried out on fixed non-
uniform grids with more nodes interpolated near the coalescence point of
$h(x,0)$ (see the inset in Figure 10). The grids are exponentially distributed
with $\triangle x_{\mathrm{min}}=0.01$ at $h_{\mathrm{min}}(t_{0})$ and
$\triangle x_{\mathrm{max}}=0.1$ on the precursor film on both sides to
capture large gradients near the coalescence point. Additionally, we set the
spatial correlation length $L_{\mathrm{c}}=0.1\geq\triangle x$ and use a
constant time step $\triangle t=10^{-3}$. The numerical process is stable with
this time step for the configuration of coalescence, so the temporally
correlated noise model is not activated here.

Figure 11: Time evolution of the bridge height during drop coalescence:
comparison between MD and the numerical solutions for the LE and SLE, for
symmetric-drop coalescence in (a) and asymmetric-drop coalescence in (b). The
SLE results are an average from 50 realisations and the MD ones come from 20
realisations. The error bars and shadows in (a) and (b) represent the standard
deviations of the MD and the SLE, respectively. Figure 12: Coalescence profile
predicted by the SLE/LE and the MD at three time instants, $t_{1}$ (black
lines), $t_{2}$ (blue lines) and $t_{3}$ (red lines). (a) symmetric-drop
coalescence: $t_{1}=t_{0}+0.83$, $t_{2}=t_{0}+12.52$ and $t_{3}=t_{0}+22.15$.
(b) Asymmetric-drop coalescence: $t_{1}=t_{0}+0.83$, $t_{2}=t_{0}+6.37$ and
$t_{3}=t_{0}+13.85$.

The time evolution of the the minimum bridge height, $h_{\mathrm{min}}(t)$, is
shown in Figure 11, where the MD coalescence time, $t_{0}$ is set as zero for
the comparison. From a theoretical aspect, $h_{\mathrm{min}}(t_{0})$ is
expected to be zero. However in MD, with finite-sized particles, the thickness
of the initial ‘molecular bridge’, $h_{\mathrm{min}}(t_{0})$, is approximately
equal to the molecular scale. In both symmetric and asymmetric cases, good
agreement between SLE and MD predictions is found at all times for the mean
values, but also, importantly, for the standard deviation. Notably, the
deterministic model (LE) is not able to accurately capture the dynamics at
these scales, highlighting the significant role that thermal fluctuations must
play. Moreover, the stochastic $h_{\mathrm{min}}(t)$ always appears smaller
than that predicted by the (deterministic) LE for the same time, demonstrating
that noise decelerates the coalescence. This finding is contrary to previous
results for other nanoscale interfacial flows, where typically thermal
fluctuations accelerate processes; such as in jet instability [9], thread
rupture [34] and droplet spreading [23]. Therefore, it is concluded that the
thermal noise cannot always assumed to be a driving force in the interface
dynamics at the nanoscale, and its role is determined by the particular fluid
configuration. Hernández-Sánchez et al. proposed a power law from the LE to
describe coalescence dynamics and demonstrated their model with experiments in
[4]. However, this power law is not found in Figure 11, even in the LE
solutions (black dashed lines in Figure 11). At present, the reason is unclear
and should be the subject of future investigation.

The ensemble-averaged profiles plotted in Figure 12 show good overall
agreement between the MD results (solid lines) and the SLE solutions (dashed
lines). For the asymmetric case, the averaged bottom point
($h_{\mathrm{min}}$) moves towards the smaller droplet with time, driven by
the surface tension. As noted above, the deterministic predictions (dash-
dotted lines) do not capture the details of the profiles of the MD results,
while the SLE results show remarkably good agreement. Moreover, the SLE
solution can reproduce the MD result at a fraction of the computational cost
of MD. For the asymmetric-drop coalescence simulation in this section, one MD
realization takes more than 200 core hours, while one SLE simulation takes
less than 0.2 core hour; i.e., the speed-up is about $10^{3}$.

###  IV.3 Thin-film rupture

The bounded film can become unstable (and rupture) due to the van der Waals
forces, as found in experiments with polymeric liquids [42]. This phenomenon
has been successfully described by the LE [42] with a similarity solution,
$h_{\mathrm{min}}(t)\sim(t_{\mathrm{b}}-t)^{1/5}$ [43], where the van der
Waals forces are modelled by the disjoining pressure term
($\partial\Pi/\partial x$) in Equation (1). As the scale at the final stage of
the rupture reaches several nanometers, thermal fluctuations are expected to
play a significant role in rupture with a potentially modified similarity
solution due to the fluctuations, which is, remarkably, not included in
previous studies [42, 21, 22, 38]. Therefore, we use our SLE solver in this
section to investigate the influence of thermal fluctuations on the rupture
dynamics.

Figure 13: Thin-film rupture profiles at three time instants. The dashed lines
are the (deterministic) LE result with $t_{1}=1.0$ (black), $t_{2}=51.0$ (red)
and $t_{3}=52.6$ (blue). The solid lines are from one realisaton for the SLE
($\varphi=10^{-3}$) with $t_{1}=1.0$ (black), $t_{2}=13.0$ (red) and
$t_{3}=14.0$ (blue). The inset shows the adaptive grids at the rupture point.

We set the correlated length $L_{\mathrm{c}}=0.15$ and correlated time scale
$T_{\mathrm{c}}=10^{-5}$. The dimensionless Hamaker constant $A=0.2$. As
$\Pi\sim 1/h^{3}$, the disjoining pressure term increases rapidly at the final
stage of the rupture, resulting in two challenges for capturing the dynamics
numerically: (i) the spatial derivatives are very large at the rupture point
(see the sharp ‘spikes’ in Figure 13) and (ii) the breakup happens extremely
fast at the final stage (the time for the final breakup, $t_{3}-t_{2}$, is
much smaller than the time for the early perturbation development,
$t_{2}-t_{1}$, in Figure 13). To overcome these challenges, we employ both
spatial and temporal refinement with the simple criteria $\triangle x\sim
h(x,t)$ and $\triangle t\sim h_{\mathrm{min}}(t)$. Initially, $\triangle
x=0.15$ and $\triangle t=10^{-5}$, and they are decreased automatically to
$1.5\times 10^{-4}$ and $10^{-8}$ respectively, to capture the dynamics at the
length scale of $10^{-2}h_{0}$ (see the inset in Figure 14).

Figure 14: The temporal evolution of the minimum film height for different
values of $\varphi$. Here, green lines represent the simulation result of LE.
Black, red and blue lines are the average of 50 SLE realisations, where the
error bars are the standard deviation. Inset shows the comparison with the
similarity solution [43].

Figure 14 shows the time evolution of the minimum film height with different
fluctuation intensities. Since we focus on the dynamics near rupture,
$h_{\mathrm{min}}$ is plotted against time to rupture, $t_{\mathrm{r}}-t$,
where $t_{\mathrm{r}}$ is the time of rupture. The error bars represent the
standard deviation from 50 independent SLE realisations. It is not surprising
to find stronger fluctuations (larger $\varphi$) lead to larger standard
deviations, but, it transpires that, these do not affect the ‘dynamic path’
significantly in this regime (see the mean values in Figure 14). A further
interesting finding is that all numerical predictions (with different values
of $\varphi$) match the similarity solution proposed by Zhang & Lister [43]
fairly well to the instant of rupture (see inset in Figure 14), indicating
that van der Waals forces, rather than thermal fluctuations, dominate the
final stages of thin-film rupture at the nanoscale. This finding is perhaps
not surprising given the singular nature of the disjoining pressure term. The
conclusion is also supported by the profiles in Figure 13, where the
stochastic spike is very similar to the deterministic one, despite the obvious
oscillations located on two sides of the film. Therefore, in this particular
regime, thermal fluctuations can play a role in accelerating instability
generation, as described in Ref. [17, 22, 20], but once disjoining pressure
become significant, it overwhelmes the influence of fluctuations and dominates
the dynamics.

##  V Conclusions and future work

In this work, a simple yet effective numerical scheme for the SLE has been
developed to predict the interface dynamics of different classes of nanoscale
bounded thin-film flows. The Fornberg scheme and correlated-noise model is
employed in the solver for the non-uniform adaptive grids, offering the
capability of capturing local dynamics accurately and efficiently. Based on
verification with theoretical models and comparisons to MD results, this
solver is demonstrated to be a powerful tool for studying both linear and
nonlinear thin-film flows.

Potential directions for future research are related to the physics predicted
by the SLE in spreading, coalescence and rupture. In this article, thermal
fluctuations are found to (i) accelerate the droplet spreading, (ii)
surprisingly decelerate the bounded droplet coalescence and (iii) not affect
the final stages of film rupture. In the future, it will be interesting to see
whether these findings can be verified experimentally by using colloid-polymer
mixture with ultra-low surface tension, which has been applied to generate
significant thermal fluctuations in interfacial flows at several micrometers
[7, 13, 18, 14]. Additionally, the coalescence dynamics has been shown to be
described by the analytical (similarity) solution in the deterministic (LE)
cases [4]. So a similar analytical solution is expected to predict the
‘fluctuation-dominated’ coalescence, which would be worthy of further
investigations.

Moreover, the results presented in this paper are obtained for two-dimensional
flows. Interesting dynamics in three-dimensional thin-film flows, such as
different film-rupture patterns [42] and fingering instability in
wetting/dewetting [44, 2], has been thoroughly studied by the three-
dimensional LE. Therefore, it would be interesting to extend our numerical
scheme to the three-dimensional SLE and explore the influence of thermal
fluctuations on these three-dimensional fluctuating hydrodynamics.

###### Acknowledgements.

This work was supported by the EPSRC (Grants No. EP/N016602/1, No.
EP/P020887/1, No. EP/S022848/1 and No. EP/P031684/1).

##  Appendix A The Fornberg scheme

Here we present a quick derivation of the Fornberg scheme [31] to show that
the errors in the derivative approximations can be well controlled. Suppose we
are given the values of target function $h(x)$ at $x_{0}<x_{1}<\ldots<x_{n}$
and we would like to approximate the $k$th derivative of $h(x)$ at
$x\in\Omega$ where $\Omega=[x_{0},x_{n}]$. The Lagrange interpolation
polynomial of $h(x)$ based on ${h_{i}=h(x_{i})}_{i=0}^{n}$, $i=0,1,\ldots,j$
is given by

$p_{n}(x)=\sum_{i=0}^{n}\mathcal{L}_{i,n}h_{i},$ (22)

where $\\{\mathcal{L}_{i,n}\\}_{i=1}^{n}$ are the Lagrange polynomials defined
by

$\mathcal{L}_{i,n}(x)=\prod_{j=0,j\neq i}^{n}\frac{x-x_{j}}{x_{i}-x_{j}}.$
(23)

It is then immediately known that if $h(x)$ is $n+1$ times continuously
differentiable on $\Omega$, then for each $x\in\Omega$,

$\|h^{(k)}(x)-p_{n}^{(k)}(x)\|\leq\|\pi^{(k)}(x)\|\frac{\|h^{(n+1)}(x)\|}{k!(n+1-k)!},$
(24)

where $k\leq n$, $h^{(k)}$ is the $k$th derivative of $h(x)$,
$\pi(x)=(x-x_{0})(x-x_{1})\cdots(x-x_{n})$ and $\|$ $\|$ denotes the supreme
norm on $\Omega$ [45]. Thus $h^{(k)}(x)$ can be well approximated by
$p_{n}^{(k)}(x)$ for $k\leq n$ if $|x_{n}-x_{0}|$ is small,

$h^{(k)}(x)\approx
p_{n}^{(k)}(x)=\sum_{i=0}^{n}\frac{d^{k}\mathcal{L}_{i,n}(x)}{dx^{k}}h_{i}=\sum_{i=0}^{n}\mathcal{L}_{i,n}^{k}(x)h_{i},$
(25)

where $\mathcal{L}_{i,n}^{k}(x)$ is the expression for the Fornberg
coefficients and can be calculated recursively. We know

$\mathcal{L}_{i,n}(x)=\frac{x-x_{n}}{x_{i}-x_{n}}\mathcal{L}_{i,n-1}(x),\quad\text{for}\quad
i\neq n,$ (26)

and

$\mathcal{L}_{n,n}(x)=\frac{\prod_{j=0}^{n-2}(x_{n-1}-x_{j})}{\prod_{j=0}^{n-1}(x_{n}-x_{j})}(x-x_{n-1})\mathcal{L}_{n-1,n-1}(x).$
(27)

Then for $i\neq n$

$\displaystyle\mathcal{L}_{i,n}^{k}(x)$
$\displaystyle=\frac{d^{k}}{dx^{k}}\left(\frac{x-x_{n}}{x_{i}-x_{n}}\mathcal{L}_{i,n-1}(x)\right)$
(28)
$\displaystyle=\frac{x-x_{n}}{x_{i}-x_{n}}\mathcal{L}_{i,n-1}^{k}(x)+\frac{k}{x_{i}-x_{n}}\mathcal{L}^{k-1}_{i,n-1}(x),$
(29)

and

$\displaystyle\mathcal{L}_{n,n}^{k}(x)$
$\displaystyle=\frac{d^{k}}{dx^{k}}\left(\frac{\prod_{j=0}^{n-2}(x_{n-1}-x_{j})}{\prod_{j=0}^{n-1}(x_{n}-x_{j})}(x-x_{n-1})\mathcal{L}_{n-1,n-1}(x)\right)$
(30)
$\displaystyle=\frac{\prod_{j=0}^{n-2}(x_{n-1}-x_{j})}{\prod_{j=0}^{n-1}(x_{n}-x_{j})}\left[(x-x_{n-1})\mathcal{L}_{n-1,n-1}^{k}(x)+k\mathcal{L}_{n-1,n-1}^{k-1}(x)\right].$
(31)

Since $\mathcal{L}_{0,0}$ represents a single point, it is natural to set
$\mathcal{L}_{0,0}^{0}=0$. Further more $\mathcal{L}_{j,m}^{s}=0$ for $s>m$
because the $s$th derivative of a $m$th order polynomial is always zero.
Without loss of generality we set $x=0$ to get

$\displaystyle\mathcal{L}_{i,n}^{k}$
$\displaystyle=\frac{1}{x_{n}-x_{i}}\left(x_{n}\mathcal{L}_{i,n-1}^{k}-k\mathcal{L}^{k-1}_{i,n-1}\right),$
(32) $\displaystyle\mathcal{L}_{n,n}^{k}$
$\displaystyle=\frac{\prod_{j=0}^{n-2}(x_{j}-x_{n-1})}{\prod_{j=0}^{n-1}(x_{j}-x_{n})}\left[x_{n-1}\mathcal{L}_{n-1,n-1}^{k}(x)-k\mathcal{L}_{n-1,n-1}^{k-1}(x)\right],$
(33) $\displaystyle\|h^{(k)}$ $\displaystyle-
p_{n}^{(k)}\|\leq\|\pi^{(k)}\|\frac{\|h^{(n+1)}\|}{k!(n+1-k)!}\simeq\mathcal{O}({\delta
x}^{n+1-k}),$ (34)

where $\delta x=\max\\{|x_{0}|,|x_{n}|\\}$. Generally speaking the Fornberg
scheme has at least accuracy of order $n+1-k$.

##  Appendix B Correlated noise model

In this appendix, we introduce the spatially correlated noise model, proposed
by Grün in [15], where an exponential correlation function is employed

$\displaystyle F_{\mathrm{cor}}(x,L_{\mathrm{c}})=$
$\displaystyle\left\\{\begin{array}[]{ll}X^{-1}\mathrm{exp}\left[-\frac{1}{2}\left(\frac{L}{L_{\mathrm{c}}}\mathrm{sin}\left(\pi
x/L\right)\right)^{2}\right],&\mathrm{for}\,\,L_{\mathrm{c}}>0\,,\\\
\delta(x),&\mathrm{for}\,\,L_{\mathrm{c}}=0\,.\end{array}\right.$ (37)

Figure 15: Linear spectrum of eigenvalues for several values of
$L_{\mathrm{c}}$ from Equation (38). Here, the wavenumber, $k=2\pi q/L$

Here, $L_{\mathrm{c}}$ is the spatial correlation length, $L$ is the domain
length, $X$ is such that $\int_{0}^{L}F_{\mathrm{cor}}(x,L_{\mathrm{c}})dx=1$.

Diez et al.[22] calculated the integral and found that $\chi_{q}$ could be
expressed by the Bessel function,

$\chi_{q}=I_{k}(\alpha)/I_{0}(\alpha)\,,$ (38)

where

$\alpha=\left(\frac{L}{2L_{\mathrm{c}}}\right)^{2}\,\quad\mathrm{and}\quad
k=2\pi q/L\,.$

Figure 15 shows the eigenvalue spectrum for several values of
$L_{\mathrm{c}}$. Note that for $L_{\mathrm{c}}\rightarrow 0$ (i.e.,
$\alpha\rightarrow\infty$), we have $\chi_{q}\rightarrow 1$ for all $q$,
leading to the limiting case of the white (uncorrelated) noise.

Figure 16: Spatially correlated noise with different $L_{\mathrm{c}}$. Figure
17: The temporally correlated stochastic term $N^{t}$ using linear
interpolation.

The term $g_{q}$ corresponds to the set of orthonormal eigenfunctions
according to

$\displaystyle g_{q}(x)=$
$\displaystyle\left\\{\begin{array}[]{ll}\sqrt{\frac{2}{L}}\mathrm{cos}(\frac{2\pi
qx}{L}),&\mathrm{for}\,\,q>0\\\ \sqrt{\frac{1}{L}},&\mathrm{for}\,\,q=0\\\
\sqrt{\frac{2}{L}}\mathrm{sin}(\frac{2\pi
qx}{L}),&\mathrm{for}\,\,q<0\end{array}\right.$ (42)

Therefore, the discretised expression of the noise term is

$N^{t}_{i}=\frac{1}{\sqrt{\triangle
t}}\sum^{q=\frac{M+1}{2}}_{q=-\frac{M+1}{2}}\chi_{q}\,N^{t}_{q}\,g_{q}(x)\,,$
(43)

where $M$ is the number of nodes. Samples of $N^{t}_{i}$ are illustrated in
Figure 16 with different spatial correlation lengths. Note that a larger
$L_{\mathrm{c}}$ leads to smooth large-wavelength and small-amplitude noise.

The temporally correlated noise model proposed by Zhao et al. [34] is employed
in this work, shown in Figure 17, where $T_{\mathrm{c}}$ is the correlation
timescale. When the adaptive time steps are reduced to less than
$T_{\mathrm{c}}$ to capture local dynamics (e.g., for the final stage of the
film rupture), this temporally correlated model is activated with time step
added using linear interpolation. The discretised expression of the noise term
becomes

$N^{t}_{i}=\frac{1}{\sqrt{T_{\mathrm{c}}}}\sum^{q=\frac{M+1}{2}}_{q=-\frac{M+1}{2}}\chi_{q}\,N^{t}_{q}\,g_{q}(x)\,.$
(44)

##  Appendix C Derivation of theory for thermal capillary waves

Without the disjoining pressure term, Equation (1) becomes,

$\partial_{t}h=-\partial_{x}\left(h^{3}\partial_{x}^{3}h-\sqrt{2\varphi
h^{3}}\,\mathcal{N}\right)\,.$ (45)

For the linear instability, we set $h=1+\hat{h}$ with $\hat{h}\ll 1$ to
linerize Equation (45):

$\partial_{t}\hat{h}+\partial_{x}^{4}\hat{h}^{4}=\sqrt{2\varphi}\,\partial_{x}\mathcal{N}.$
(46)

Then a finite Fourier transform is applied to Equation (46) to get

$\partial_{t}H+k^{4}H=ik\sqrt{2\varphi}N\,,$ (47)

where the transformed variables are defined as follows:

$H(k,t)=\int^{L}_{0}\hat{h}(x,t)e^{-ikx}dx\quad\textrm{and}\quad
N(k,t)=\int^{L}_{0}\mathcal{N}(x,t)e^{-ikx}dx.$

The solution of Equation (47) is linearly decomposed into two parts:

$H={H}_{\mathrm{LE}}+{H}_{\mathrm{fluc}}\,.$ (48)

The first part is the solution to the homogenous form of Equation (47) (i.e.,
with $\varphi=0$) with some stationary initial disturbances (i.e.,
$H$=$H_{\mathrm{i}}$ at $t=0$). Since we start from a smooth initial surface,
$H_{\mathrm{LE}}(k,t)=0$. The second component of the solution arises from
solving the full form of Equation (48) without any initial disturbances; this
part of the solution is solely due to fluctuations, and is thus denoted
$H_{\mathrm{fluc}}$. This is obtained by determining the homogeneous
equation’s impulse response,

$\displaystyle H_{\mathrm{res}}(k,t)=e^{-k^{4}t}\,,$ (49)

which due to the linear, time-invariant nature of the system, allows to write

$\displaystyle\overline{|H_{\mathrm{fluc}}|^{2}}$
$\displaystyle=\left(ik\sqrt{2\varphi}\right)^{2}\overline{\left|\int^{t}_{0}N(k,t-\tau)H_{\mathrm{res}}(k,\tau)d\tau\right|^{2}}$
(50)

$H$ is both a random and complex variable with zero mean. So the root mean
square (rms) of $H$ is sought, which from Equation (48) is given by

$\displaystyle|H|_{\mathrm{rms}}$
$\displaystyle=\sqrt{\overline{|H_{\mathrm{LE}}+H_{\mathrm{fluc}}|^{2}}}=\sqrt{\overline{|H_{\mathrm{fluc}}|^{2}}}\,.$
(51)

Because $N$ is uncorrelated Gaussian white noise, and the variance of the norm
of the white noise $\overline{|N|^{2}}=L$, Equation (49) and (50) combine to
give

$\displaystyle\overline{|H_{\mathrm{fluc}}|^{2}}$ $\displaystyle=2\varphi
k^{2}\int^{t}_{0}\overline{|N(k,t-\tau)|^{2}}H(k,\tau)^{2}d\tau\,,$
$\displaystyle=2\varphi k^{2}L\int^{t}_{0}H^{2}d\tau\,,$
$\displaystyle=-\frac{\varphi L}{k^{2}}\left(e^{-2k^{4}t}-1\right)\,.$ (52)

Equation (51) and (52) constitute the theory for thermal capillary waves used
in Sec. III.

##  Appendix D Derivation for the spreading power-laws

In this appendix, we present the derivations for both no-slip and slip
theories, i.e., $\ell\sim t^{n}$. The average second moment of $h$ is applied
to estimate the drop’s width $\ell$ [23] whose explicit expression is

$\ell(t)=\left\langle\left[\frac{1}{V}\int(x-X)^{2}h(x,t)dx\right]\right\rangle\,,$
(53)

where $V=\int hdx$ is the constant volume of the droplet, $X=\left(\int
xhdx\right)/V$ is the instantaneous position of the droplet center, and
$\left\langle\ldots\right\rangle$ represents the ensemble average of all the
realisations.

In order to decide the power law ($n$), a similarity transform is made with
the change of variables:

$x=b\breve{x},\qquad h=b^{\alpha}\breve{h},\qquad t=b^{\eta}\breve{t},$ (54)

where $b$ is an arbitrary factor, and $\alpha$ and $\eta$ are constants that
remain to be fixed. The symbol ‘$\,\breve{\textbullet}\,$’ means ‘transformed’
variables. According to the scaling relation above
($\breve{x}\sim\breve{h}\sim\breve{t}\sim\mathrm{O}(1)$), we can easily obtain
$\ell\sim x\sim t^{1/\eta}$, namely, the power law is equal to $1/\eta$. For
the value of $\eta$, we substitute the transform relations above into the SLE
(Equation (1)) and obtain

$\frac{\partial\breve{h}}{\partial\breve{t}}=-\left(b^{3\alpha+\eta-4}\right)\partial_{\breve{x}}\left(\breve{h}^{3}\partial^{3}_{\breve{x}}\breve{h}\right)+\left[b^{(\alpha+\eta-3)/2}\right]\sqrt{2\varphi}\partial_{\breve{x}}\left[\breve{h}^{3/2}\breve{\mathcal{N}}(\breve{x},\breve{t})\right]\,.$
(55)

Here, we have two independent force terms on the right-hand side: (i) the
deterministic term due to the surface tension and (ii) the stochastic term due
to the thermal fluctuations. To hold the ‘similarity’ of the transform, the
scaling powers of the arbitrary scaling factor $b$ should always be zero,
namely,

$\displaystyle\left\\{\begin{array}[]{ll}3\alpha+\eta-4=0\,,\quad\textrm{for
the deterministic term,}\\\ \alpha+\eta-3=0\,,\quad\textrm{for the stochastic
term.}\end{array}\right.$ (58)

In addition, no matter which force drives the spreading, the droplet volume,

$V=\int hdx=b^{\alpha+1}\int\breve{h}d\breve{x}\,,$

should always be conserved, requiring $\alpha=-1$. So, we can obtain the value
of the left coefficient $\eta$ from Equation (58), i.e., $\eta=7$ in the
surface tension term; and $\eta=4$ in the stochastic term, implying two power-
law spreading regimes:

$\displaystyle\left\\{\begin{array}[]{ll}\ell\sim
t\,^{1/7},\quad\textrm{Tanner's law,}\\\ \ell\sim
t\,^{1/4},\quad\textrm{fluctuation enhanced Tanner's law},\end{array}\right.$
(61)

which have been proposed by Tanner [3] and Davidovitch et al. [23],
respectively.

To take the slip effect into account, Equation (58) is modified:

$\frac{\partial\breve{h}}{\partial\breve{t}}=-\left(b^{2\alpha+\eta-4}\right)\partial_{\breve{x}}\left(3\tilde{\ell_{\mathrm{s}}}\breve{h}^{2}\partial^{3}_{\breve{x}}\breve{h}\right)+b^{(\eta-3)/2}\sqrt{2\varphi}\partial_{\breve{x}}\left[\sqrt{3\tilde{\ell_{\mathrm{s}}}\breve{h}^{2}}\breve{\mathcal{N}}(\breve{x},\breve{t})\right].$
(62)

By the same approach, we get

$\displaystyle\left\\{\begin{array}[]{ll}2\alpha+\eta-4=0\,,\quad\textrm{for
the deterministic term,}\\\ \eta-3=0\,,\quad\textrm{for the stochastic
term,}\end{array}\right.$ (65)

where $\alpha$ is still equal to $-1$. Therefore, we can obtain a ‘slip-
modified’ power law for the spreading:

$\displaystyle\left\\{\begin{array}[]{ll}\ell\sim t\,^{1/6},\quad\textrm{Slip-
modified Tanner's law,}\\\ \ell\sim t\,^{1/3},\quad\textrm{Slip-modified
fluctuation enhanced Tanner's law}.\end{array}\right.$ (68)

##  Appendix E MD settings for the droplet coalescence

In this work, we choose the mW model [46] to simulate liquid water. The model
mimics the hydrogen-bonded structure of water through the introduction of a
non-bond angular dependent term that encourages tetrahedral configurations.
The model contains two terms: (i) $\phi_{ij}$ depending on the distances
between pairs of atoms (represented by $r_{ij}$ and $s_{ik}$) and (ii)
$\phi_{ijk}$ depending on the angles formed by triplets of atoms (represented
by $\theta_{ijk}$). The full expression is given by

$\displaystyle\mathcal{U}$
$\displaystyle=\sum\limits_{i}\sum\limits_{j>i}\phi_{ij}(r_{ij})+\sum\limits_{i}\sum\limits_{j\neq
i}\sum\limits_{k>j}\phi_{ijk}(r_{ij},s_{ik},\theta_{ijk})\,,$ (69)
$\displaystyle\phi_{ij}(r_{ij})=A\epsilon\left[B\left(\frac{\sigma}{r_{ij}}\right)^{p}-\left(\frac{\sigma}{r_{ij}}\right)^{q}\right]\mathrm{exp}\left(\frac{\sigma}{r_{ij}-a\sigma}\right)\,,$
$\displaystyle\phi_{ijk}(r_{ij},s_{ik},\theta_{ijk})=\kappa\epsilon(\cos\theta_{ijk}-\cos\theta_{0})^{2}\,\mathrm{exp}\left(\frac{\chi\sigma}{r_{ij}-a\sigma}\right)\,\mathrm{exp}\left(\frac{\chi\sigma}{s_{ik}-a\sigma}\right)\,,$

where $A$, $B$, $p$, $q$, $\chi$ and $\kappa$ respectively give the form and
scale to the potential, and $\theta_{0}$ represents the tetrahedral angles.
All the parameters are presented in Table 1.

Table 1:  Parameters of the mW model $\epsilon$ (kJ mol-1) |  $\sigma$ (nm) | $A$ | $B$ | $p$ | $q$ | $\chi$ | $\kappa$ | $a$ |  $\theta_{0}$ (degree)
---|---|---|---|---|---|---|---|---|---
25.87 | 0.2390 | 7.050 | 0.6022 | 4 | 0 | 1.2 | 23.15 | 1.8 | $109.47$

The platinum substrate is assumed to be rigid with an atomic mass of
$3.24\times 10^{-25}\,\mathrm{kg}$ [47]. The liquid-solid interaction is
modelled by the 12-6 LJ potential with
$\epsilon_{\mathrm{ls}}/k_{\mathrm{B}}=444\,$K and
$\sigma_{\mathrm{ls}}=0.28\,$nm to create a fully wettable substrate (zero
contact angle). The initial configurations of droplets are cut from a liquid
bulk, created from equilibrium NVT simulations with a Nosé-Hoover thermostat
at $T=400\,$K. The same ensemble and thermostat is used for the main
simulations with the time step, 2.5 femtoseconds. The entire domain width
along the $z$-axis, $W=2\,$nm and the length scale for nondimensionalisation,
$h_{0}=10\,$nm.

To compare MD results with the predictions of the SLE, the liquid transport
properties are calculated are required. Here, dynamic viscosity is found to be
$\mu=1.64\times 10^{-4}\,\mathrm{kg}\,\mathrm{m}^{-1}\,\mathrm{s}^{-1}$ by the
Green Kubo method [48, 49], which integrating the time-autocorrelation
function of the off-diagonal elements of the pressure tensor $P_{ij}$ so that

$\mu=\frac{V_{\mathrm{bulk}}}{k_{\mathrm{B}}T}\int_{0}^{\infty}\left\langle
P_{ij}(t)\cdot P_{ij}(0)\right\rangle dt\quad(i\neq j),$ (70)

where $V_{\mathrm{bulk}}$ is the volume of the bulk fluid, $k_{\mathrm{B}}$ is
the Boltzmann constant, and $T$ is temperature. The pressure tensor components
are obtained using the definition of [50] and the angular brackets indicate
the expectation. The surface tension is calculated from the profiles of the
components of the pressure tensor in a simple liquid-vapour system, using the
mechanical definition [51]:

$\gamma=\frac{1}{2}\int_{0}^{L_{z}}\left[P_{\mathrm{n}}(z)-P_{\mathrm{t}}(z)\right]dz,$
(71)

where $L_{z}$ is the length of the MD domain, and subscripts ‘n’ and ‘t’
denote normal and tangential components, respectively. Finally, we have
$\gamma=5.45\times 10^{-2}\,\mathrm{N}\,\mathrm{m}^{-1}$.

## References

  * Jacobs  _et al._ [1998] K. Jacobs, S. Herminghaus, and K. R. Mecke, Thin liquid polymer films rupture via defects, Langmuir 14, 965 (1998).
  * Bäumchen  _et al._ [2014] O. Bäumchen, L. Marquant, R. Blossey, A. Münch, B. Wagner, and K. Jacobs, Influence of slip on the rayleigh-plateau rim instability in dewetting viscous films, Phys. Rev. Lett. 113, 014501 (2014).
  * Tanner [1979] L. H. Tanner, The spreading of silicone oil drops on horizontal surfaces, J. Phys. D: Appl. Phys. 12, 1473 (1979).
  * Hernández-Sánchez  _et al._ [2012] J. F. Hernández-Sánchez, L. A. Lubbers, A. Eddi, and J. H. Snoeijer, Symmetric and asymmetric coalescence of drops on a substrate, Phys. Rev. Lett. 109, 184502 (2012).
  * Oron  _et al._ [1997] A. Oron, S. H. Davis, and S. G. Bankoff, Long-scale evolution of thin liquid films, Rev. Mod. Phys. 69, 931 (1997).
  * Stone  _et al._ [2004] H. A. Stone, A. D. Stroock, and A. Ajdari, Engineering flows in small devices: microfluidics toward a lab-on-a-chip, Annu. Rev. Fluid Mech. 36, 381 (2004).
  * Aarts  _et al._ [2004] D. G. Aarts, M. Schmidt, and H. N. Lekkerkerker, Direct visual observation of thermal capillary waves, Science 304, 847 (2004).
  * Moseler and Landman [2000] M. Moseler and U. Landman, Formation, stability, and breakup of nanojets, Science 289, 1165 (2000).
  * Zhao  _et al._ [2019] C. Zhao, J. E. Sprittles, and D. A. Lockerby, Revisiting the rayleigh-plateau instability for the nanoscale, J. Fluid Mech. 861 (2019).
  * Perumanath  _et al._ [2019] S. Perumanath, M. K. Borg, M. V. Chubynsky, J. E. Sprittles, and J. M. Reese, Droplet coalescence is initiated by thermal motion, Phys. Rev. Lett. 122, 104501 (2019).
  * Zhang  _et al._ [2021] Y. Zhang, J. E. Sprittles, and D. A. Lockerby, Thermal capillary wave growth and surface roughening of nanoscale liquid films, J. Fluid Mech. 915 (2021).
  * Eggers [2002] J. Eggers, Dynamics of liquid nanojets, Phys. Rev. Lett. 89, 084502 (2002).
  * Hennequin  _et al._ [2006] Y. Hennequin, D. Aarts, J. H. van der Wiel, G. Wegdam, J. Eggers, H. N. W. Lekkerkerker, and D. Bonn, Drop formation by thermal fluctuations at an ultralow surface tension, Phys. Rev. Lett. 97, 244502 (2006).
  * Petit  _et al._ [2012] J. Petit, D. Rivière, H. Kellay, and J. Delville, Break-up dynamics of fluctuating liquid threads, Proc. Natl. Acad. Sci. U.S.A. 109, 18327 (2012).
  * Grün  _et al._ [2006] G. Grün, K. Mecke, and M. Rauscher, Thin-film flow influenced by thermal noise, J. Stat. Phys. 122, 1261 (2006).
  * Landau  _et al._ [1980] L. D. Landau, E. M. Lifshits, and L. P. Pitaevskii, in  _Statistical Physics, Part 2: Theory of the Condensed State_ , Course of theoretical physics, Vol. 9 (Pergamon Press, 1980) Chap. 9, pp. 86–91, 2nd ed.
  * Mecke and Rauscher [2005] K. Mecke and M. Rauscher, On thermal fluctuations in thin film flow, J. Phys.: Condens Matter 17, S3515 (2005).
  * Fetzer  _et al._ [2007] R. Fetzer, M. Rauscher, R. Seemann, K. Jacobs, and K. Mecke, Thermal noise influences fluid flow in thin films during spinodal dewetting, Phys. Rev. Lett. 99, 114503 (2007).
  * Willis and Freund [2010] A. M. Willis and J. B. Freund, Thermal capillary waves relaxing on atomically thin liquid films, Phys. Fluids 22, 022002 (2010).
  * Zhang  _et al._ [2019] Y. Zhang, J. E. Sprittles, and D. A. Lockerby, Molecular simulation of thin liquid films: Thermal fluctuations and instability, Phys. Rev. E 100, 023108 (2019).
  * Nesic  _et al._ [2015a] S. Nesic, R. Cuerno, E. Moro, and L. Kondic, Fully nonlinear dynamics of stochastic thin-film dewetting, Phys. Rev. E 92, 061002 (2015a).
  * Diez  _et al._ [2016] J. A. Diez, A. G. González, and R. Fernández, Metallic-thin-film instability with spatially correlated thermal noise, Phys. Rev. E 93, 013120 (2016).
  * Davidovitch  _et al._ [2005] B. Davidovitch, E. Moro, and H. A. Stone, Spreading of viscous fluid drops on a solid substrate assisted by thermal fluctuations, Phys. Rev. Lett. 95, 244505 (2005).
  * Nesic  _et al._ [2015b] S. Nesic, R. Cuerno, E. Moro, and L. Kondic, Dynamics of thin fluid films controlled by thermal fluctuations, Eur. Phys. J. Special Top. 224, 379 (2015b).
  * Shah  _et al._ [2019] M. S. Shah, V. van Steijn, C. R. Kleijn, and M. T. Kreutzer, Thermal fluctuations in capillary thinning of thin liquid films, J. Fluid Mech. 876, 1090 (2019).
  * Diez and Kondic [2002] J. A. Diez and L. Kondic, Computing three-dimensional thin film flows including contact lines, J. Comput. Phys. 183, 274 (2002).
  * Pahlavan  _et al._ [2018] A. A. Pahlavan, L. Cueto-Felgueroso, A. E. Hosoi, G. H. McKinley, and R. Juanes, Thin films in partial wetting: Stability, dewetting and coarsening, J. Fluid Mech. 845, 642 (2018).
  * Zitz  _et al._ [2019] S. Zitz, A. Scagliarini, S. Maddu, A. A. Darhuber, and J. Harting, Lattice boltzmann method for thin-liquid-film hydrodynamics, Phys. Rev. E 100, 033313 (2019).
  * Lam  _et al._ [2019] M. Y. Lam, L. J. Cummings, and L. Kondic, Computing dynamics of thin films via large scale gpu-based simulations, J. Comput. Phys. 2, 100001 (2019).
  * Zhang  _et al._ [2020] Y. Zhang, J. E. Sprittles, and D. A. Lockerby, Nanoscale thin-film flows with thermal fluctuations and slip, Phys. Rev. E 102, 053105 (2020).
  * Fornberg [1998] B. Fornberg, Classroom note: Calculation of weights in finite difference formulas, SIAM Rev. 40, 685 (1998).
  * Kloeden and Platen [1992] P. E. Kloeden and E. Platen, Stochastic differential equations, in  _Numerical Solution of Stochastic Differential Equations_ (Springer, 1992) pp. 103–160.
  * Polyak [2006] B. T. Polyak, Newton-kantorovich method and its global convergence, J Math. Sci 133, 1513 (2006).
  * Zhao  _et al._ [2020] C. Zhao, D. A. Lockerby, and J. E. Sprittles, Dynamics of liquid nanothreads: Fluctuation-driven instability and rupture, Phys. Rev. Fluids 5, 044201 (2020).
  * Hairer and Voss [2011] M. Hairer and J. Voss, Approximations to the stochastic burgers equation, J. Nonlinear Sci. 21, 897 (2011).
  * Herminghaus  _et al._ [2008] S. Herminghaus, M. Brinkmann, and R. Seemann, Wetting and dewetting of complex surface geometries, Annu. Rev. Mater. Res. 38, 101 (2008).
  * Bonn  _et al._ [2009] D. Bonn, J. Eggers, J. Indekeu, J. Meunier, and E. Rolley, Wetting and spreading, Rev. Mod. Phys. 81, 739 (2009).
  * Durán-Olivencia  _et al._ [2019] M. A. Durán-Olivencia, R. S. Gvalani, S. Kalliadasis, and G. A. Pavliotis, Instability, rupture and fluctuations in thin liquid films: Theory and computations, J. Stat. Phys. 174, 579 (2019).
  * Russo  _et al._ [2021] A. Russo, S. P. Perez, M. A. Durán-Olivencia, P. Yatsyshin, J. A. C., and S. Kalliadasis, A finite-volume method for fluctuating dynamical density functional theory, J. Comput. Phys. 428, 109796 (2021).
  * Ristenpart  _et al._ [2006] W. D. Ristenpart, P. M. McCalla, R. V. Roy, and H. A. Stone, Coalescence of spreading droplets on a wettable substrate, Phys. Rev. Lett. 97, 064501 (2006).
  * Karpitschka and Riegler [2014] S. Karpitschka and H. Riegler, Sharp transition between coalescence and non-coalescence of sessile drops, J. Fluid Mech. 743 (2014).
  * Becker  _et al._ [2003] J. Becker, G. Grün, R. Seemann, H. Mantz, K. Jacobs, K. R. Mecke, and R. Blossey, Complex dewetting scenarios captured by thin-film models, Nat. Mater. 2, 59 (2003).
  * Zhang and Lister [1999] W. W. Zhang and J. R. Lister, Similarity solutions for van der waals rupture of a thin film on a solid substrate, Phys. Fluids 11, 2454 (1999).
  * Troian  _et al._ [1989] S. M. Troian, X. L. Wu, and S. A. Safran, Fingering instability in thin wetting films, Phys. Rev. Lett. 62, 1496 (1989).
  * Howell [1991] G. W. Howell, Derivative error bounds for lagrange interpolation: An extension of cauchy’s bound for the error of lagrange interpolation, J. Approx. Theory 67, 164 (1991).
  * Molinero and Moore [2008] V. Molinero and E. B. Moore, Water modeled as an intermediate element between carbon and silicon, J. Chem. Phys 113, 4008 (2008).
  * Shi and Dhir [2009] B. Shi and V. K. Dhir, Molecular dynamics simulation of the contact angle of liquids on solid surfaces, J. Chem. Phys. 130, 034705 (2009).
  * Green [1954] M. S. Green, Markoff random processes and the statistical mechanics of time-dependent phenomena. ii. irreversible processes in fluids, J. Chem. Phys 22, 398 (1954).
  * Kubo [1957] R. Kubo, Statistical-mechanical theory of irreversible processes. i. general theory and simple applications to magnetic and conduction problems, J. Phys. Soc. Japan 12, 570 (1957).
  * Kirkwood and Buff [1949] J. G. Kirkwood and F. P. Buff, The statistical mechanical theory of surface tension, J. Chem. Phys 17, 338 (1949).
  * Trokhymchuk and Alejandre [1999] A. Trokhymchuk and J. Alejandre, Computer simulations of liquid/vapor interface in lennard-jones fluids: Some questions and answers, J. Chem. Phys 111, 8510 (1999).


# LIMSI_UPV at SemEval-2020 Task 9: Recurrent Convolutional Neural Network for
Code-mixed Sentiment Analysis

Somnath Banerjee1, Sahar Ghannay1, Sophie Rosset1, Anne Vilnat1, Paolo Rosso2
1 LIMSI, CNRS, Université Paris-Saclay, Orsay, France
2 Universitat Politècnica de València, Spain
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
###### Abstract

This paper describes the participation of LIMSI_UPV team in SemEval-2020 Task
9: Sentiment Analysis for Code-Mixed Social Media Text. The proposed approach
competed in SentiMix Hindi-English subtask, that addresses the problem of
predicting the sentiment of a given Hindi-English code-mixed tweet. We propose
Recurrent Convolutional Neural Network that combines both the recurrent neural
network and the convolutional network to better capture the semantics of the
text, for code-mixed sentiment analysis. The proposed system obtained 0.69
(best run) in terms of F1 score on the given test data and achieved the 9th
place (Codalab username: somban) in the SentiMix Hindi-English subtask.

##  1 Introduction

00footnotetext:  This work is licensed under a Creative Commons Attribution
4.0 International Licence. Licence details:
http://creativecommons.org/licenses/by/4.0/. 00footnotetext: The code of this
work is available at https://github.com/somnath-banerjee/Code-
Mixed_SentimentAnalysis.

In this digital era, users express their personal thoughts and opinions
regarding a wide range of topics on social media platforms such as blogs,
micro-blogs (e.g., Twitter), and chats (e.g., WhatsApp and Facebook messages).
Multilingual societies like India with a decent amount of internet penetration
widely adopted such social media platforms. However, the regional language
influences the proliferation of the Hindi-English Code-Mixed (CM) data.
Sentiment analysis of these end-user data from social media is a crucial
resource for commerce and governance. However, in contrast to the classical
sentiment analysis methods, which were originally designed for dealing with
well-written product reviews, CM texts from social media often contain
misspellings (often intentional), badly cased words, letter substitutions,
ambiguities, non standard abbreviations, improper use of grammar, etc.

CM poses several unseen difficulties to natural language processing (NLP)
tasks such as word-level language identification, part-of-speech tagging,
dependency parsing, machine translation and semantic processing. In the last
few years, a number of workshops such as Linguistic Code-Switching
Workshops111https://code-switching.github.io/2020/ and shared tasks such as
Mixed Script Information Retrieval [Banerjee et al., 2020] have been organized
due to the emerging popularity of code-mixing. To promote research in this
area, Task 9 of SemEval-2020 was devoted to CM sentiment analysis in Twitter.
The goal of the task was to automatically classify the polarity of a given CM
Twitter post into one of the three predefined categories: positive, negative
and neutral. The CM languages are English-Hindi and English-Spanish; for a
more detailed description of the task see [Patwa et al., 2020].

In this paper, we present a deep learning approach, using a Recurrent
Convolutional Neural Network for the task of automatic CM sentiment
classification of tweets.

The rest of the paper is structured as follows. Section 2 provides background
in brief. Section 3 provides the system overview and Section 4 describes our
approach in detail. In Section 5, we discuss the analysis and evaluation
results for our system. We conclude our work in Section 6.

##  2 Background

Sentiment classification is the task of detecting whether a textual item
(e.g., a product review, a blog post, an editorial, etc.) expresses a POSITIVE
or a NEGATIVE opinion in general or about a given entity, e.g., a product, a
person, a political party, or a policy [Nakov et al., 2016]. Classifying
tweets according to sentiment has many applications in political science,
social sciences, market research, and many others [Martínez-Cámara et al.,
2014, Mejova et al., 2015]. Although initially sentiment identification was
focused on newswire text [Baccianella et al., 2010], later research turned
towards social media [Rosenthal et al., 2015]. Since 2013, a sentiment
classification task on Twitter data have been organized in SemEval campaigns.

Most of the earlier approaches to this problem were based on hand crafted
features and sentiment lexicons [Pak and Paroubek, 2010, Mohammad et al.,
2013]. These features were then used as input to classifiers (e.g., Support
Vector Machines). However, such approaches required extensive domain
knowledge, were laborious to define, and can lead to incomplete or over-
specific features.

Recently, researchers pay their attention to sentiment polarity detection on
CM data. However, a few research work have been carried out in particular
Hindi-English CM data with different approaches: lexicon lookup [Sharma et
al., 2015], sub-word with CNN-LSTM [Joshi et al., 2016], Siamese networks
[Choudhary et al., 2018], dual Encoder Network with features [Lal et al.,
2019]. Lai et al. [Lai et al., 2015] proposed Recurrent Convolutional Network
for text classification which is a foundational task in many NLP applications.
We followed this model in our task.

##  3 System overview

We are inspired by the model proposed in [Lai et al., 2015] particularly
proposed for the text classification task. The proposed model takes sequence
of CM words as input and provides sentiment polarity class as output. The
recurrent structure of the proposed model captures the contextual information
during the learning of the word representation, and the max-pooling layer
identifies the key CM words. If $T$, $s$ and $\theta$ denote a CM tweet made
up of sequence of CM words (T = $w_{1}$, $w_{2}$, $w_{3}$, …, $w_{n}$), the
sentiment polarity class and parameters of the neural network respectively,
$p(s|T,\theta))$ denotes the probability of the tweet $T$ having sentiment
polarity $s$, where $s$ could be any one of sentiment polarity classes, i.e.,
$s\in${positive, negative, neutral}.

In this model, a CM word is represented by combining the word and its both
side contexts. If $c_{l}(w_{i})$ and $c_{r}(w_{i})$ denote the left-side and
right-side contexts of a CM word $w_{i}$, a word is represented as a
concatenation of the left-side context vector $c_{l}(w_{i})$, the word
embedding $e(w_{i})$ and the right-side context vector $c_{r}(w_{i})$.

$x_{i}=[c_{l}(w_{i});e(w_{i});c_{r}(w_{i})]$

where, $x_{i}$ is the representation of the $i^{th}$ word (i.e., $w_{i}$) in
$T$.

Both the contexts (i.e., $c_{l}(w_{i})$ and $c_{r}(w_{i})$) are calculated as
follows:

$c_{l}(w_{i})=f(W^{(l)}c_{l}(w_{i-1})+W^{(s_{l})}e(w_{i-1}))$
$c_{r}(w_{i})=f(W^{(r)}c_{r}(w_{i+1})+W^{(s_{r})}e(w_{i+1}))$

Where,

$f$: is a non-linear activation function;

$e(w_{i-1})$: the word embedding of word $w_{i-1}$;

$W^{(l)}$: a matrix that transforms the hidden layer (context) into the next
hidden layer;

$W^{(sl)}$: a matrix that is used to combine the semantic of the current word
with the next word’s left context.

In the forward scan, the recurrent structure of the model obtains all the
$c_{l}$ of the CM tweet, whereas, it obtains all the $c_{r}$ in a backward
scan of the CM tweet. After obtaining the $x_{i}$ for the word $w_{i}$, a
linear transformation together with the tanh activation function is applied to
$x_{i}$ and the result is sent to the next layer:

$y^{(2)}_{i}=\textrm{tanh}(W^{(2)}x_{i}+b^{(2)})$

After calculating all of the representations of words, a max-pooling layer is
applied. Hence, the pooling layer utilizes the output of the recurrent
structure as the input. This layer attempts to find the most important latent
semantic factors in the CM tweet:

$y^{(3)}=\max^{n}_{i=1}y^{(2)}_{i}$

Finally, in the output layer, a softmax function is applied that provides the
sentiment polarity probabilities:

$y^{(4)}=W^{(4)}y^{(3)}+b^{(4)}$
$p_{i}=\frac{\textrm{exp}(y_{k}^{(4)})}{\displaystyle\sum_{k=1}^{n}\textrm{exp}(y_{k}^{(4)})}$

Figure 1: Recurrent Convolutional Neural Network [Lai et al., 2015] model

The Figure 1 depicts the architecture of the proposed model. From Figure 1, we
could see that $c_{l}(w_{n})$ captures the left-side context, whereas,
$c_{l}(w_{1})$ captures the right-side context of a word $w_{i}$.

##  4 Experimental setup

Data: For sentiment analysis on Hindi-English CM tweets, we used the dataset
provided by the organizers of Task 9 at SemEval-2020. The training dataset
consists of 14 thousand tweets. Whereas, the validation dataset as well as the
test dataset contain 3 thousand tweets each. The details of the dataset are
given in [Patwa et al., 2020]. For this task, we did not use any external
dataset.

Preprocessing: The CM Twitter data is not formal in nature. Also, often casing
is not followed properly. As a result, it misleads the language processing
task. Therefore, we converted the entire tweets into lower case. After
observing the CM dataset, we found that the urls are not contributing any
information. While writing, the writer of a tweet addresses a person by
mentioning the person’s name followed by @ or by twitter id which starts with
@. Similarly, to specify a particular topic, the writer uses the topic name
that starts with #. To capture the semantics and syntax of a tweet, we
replaced the person and a topic being addressed with a MENTION and a TOPIC
token respectively. The following steps were applied on each CM tweet:

  * •

Tokens were converted to lowercase.

  * •

If a token is http or https, the following tokens are merged to identify as
links and were deleted.

  * •

Garbage tokens (e.g., â€¦, $\grave{a}$ ¥¥, etc.) were deleted.

  * •

If a token is @, the following token is merged with the first token (i.e., @)
and the new token is replaced with a MENTION token. For example, the tokens @
and bomanirani are merged and replaced with a MENTION token.

  * •

If a token is #, the following token is merged with the first token (i.e., #)
and the new token is replaced with a TOPIC token. For example, the tokens #
and LoveIsLove are merged and replaced with a TOPIC token.

  * •

Emoji’s with text were divided into two tokens. For example, he☺ becomes he
and ☺.

  * •

If a token contains more than one emoji, each emoji was considered as a token.
For example, ☺☹ becomes ☺ and ☹.

Embeddings: Following Collobert et al. [Collobert et al., 2011], a lot of
authors argued that word embedding plays a vital role to improve natural
language task performance. Hence, we experimented the use of word embeddings
to improve the performance of our proposed models. Using the fastText
[Bojanowski et al., 2017], we prepared two embedding models: Skip-gram and
Cbow. After empirically evaluating the performance on validation set, the
embeddings‘ dimensionality was set to 300 for all the embeddings. The
embeddings are trained on training data using the parameters: lr=0.05, context
window=5, epochs=5, minimal number of word occurences=5, dimensionality=300.

Experiment: We carried out two experiments with similar settings except
different word embedding approaches: Skip-gram for SkipGRun, and Cbow for
CbowRun.

Hyper-parameters: After evaluating the model performance on the validation
data, the optimal values of the hyper-parameters were set. We used the
following list of hyper-parameters: learning rate = 0.6, word embedding vector
size = 300, hidden layers= 2, hidden layer size = 64, context vector size = 5,
dropout rate = 0.1, optimizer = stochastic gradient descent, loss function=
negative log likelihood, and batch size = 64.

##  5 Results and Analysis

As mentioned in the previous section, we carried out experiments with our
proposed model using two settings. All the presented experiments are evaluated
on the test data for the given task. The performance of the systems were
evaluated using F1 averaged across the positive, negative and the neutral.

| Precision | Recall | F1 score
---|---|---|---
Baseline | - | - | 0.654
SkipGRun | 0.6952 | 0.6893 | 0.6913
CbowRun | 0.6566 | 0.6556 | 0.6560
Table 1: Evaluation results on test data

The organizer baseline F1 scores for the validation and test data are 0.58 and
0.654 respectively. The details of the baseline are given in [Patwa et al.,
2020]. The obtained results with our submitted runs are given in Table 1. For
SkipGRun, we achieved 0.6913 of F1 score with 0.6952 and 0.6893 of precision
and recall respectively. The SkipGRun outperformed the CbowRun by around 0.40
in terms of F1 score. CbowRun outperformed the organizers’ baseline by a
slight margin, however, SkipGRun outperformed the baseline by around 0.4 in
terms of averaged F score.

For SkipGRun, the confusion matrix and the performance of class-wise sentiment
polarities are presented in Table 3 and Table 3, respectively. The SkipGRun is
the best among the two submitted runs. From the confusion matrix (cf. Table
3), we can observe that the system is more successful in identifying the
positive polarity class (F1 score 0.76) than the polarity classes: negative
(F1 score 0.71) and neutral (F1 score 0.62). This run performed well on
distinguishing a positive tweet from negative tweet and vice versa. Therefore,
only 35 positive tweets were misclassified as negative and 24 negative tweets
were misclassified as positive. However, it has a problem to separate positive
and negative tweets from neutral. Hence, a lot of positive and negative tweets
were misclassified as neutral. From Table 3, we could say that around 24% of
positive tweets and 28% of negative tweets were identified as neutral.
Similarly, a number of neutral tweets were misclassified as positive and
negative. From Table 3, it is evident that the positive tweets were the
easiest class to predict. The F1 scores were 0.76 (positive class), 0.71
(negative class) and 0.62 (neutral class).

| Positive | Negative | Neutral
---|---|---|---
Positive | 729 | 35 | 236
Negative | 24 | 624 | 252
Neutral | 175 | 210 | 715
Table 2: SkipGRun: Confusion matrix

| Precision | Recall | F1 score
---|---|---|---
Positive | 0.79 | 0.73 | 0.76
Negative | 0.72 | 0.69 | 0.71
Neutral | 0.59 | 0.65 | 0.62
Table 3: SkipGRun: Polarity class-wise performance

Table 5 and Table 5 present the confusion matrix and the performance of class-
wise sentiment polarities for CbowRun. Like SkipGRun, the CbowRun is also able
to successfully identify the positive polarity class (F1 score 0.72) in
comparison to others (cf. Table 5). For CbowRun, we observed similar
characteristics of results like SkipGRun. Although the results obtained for
negative tweets are almost similar with SkipGRun (F1 score 0.71) and CbowRun
(F1 score 0.70), the difference in performance is notable for positive and
neutral tweets. From Table 5, we can see that for CbowRun the F1 scores were
0.73 (positive class), 0.70 (negative class) and 0.57 (neutral class).

| Positive | Negative | Neutral
---|---|---|---
Positive | 709 | 41 | 250
Negative | 33 | 630 | 237
Neutral | 232 | 240 | 628
Table 4: CbowRun: Confusion matrix

| Precision | Recall | F1 score
---|---|---|---
Positive | 0.73 | 0.71 | 0.72
Negative | 0.69 | 0.70 | 0.70
Neutral | 0.56 | 0.57 | 0.57
Table 5: CbowRun: Polarity class-wise performance

To get a deeper analysis of the results, we also performed the error analysis.
We observed that the system could not identify the sentiment when Hindi lyrics
are used in tweets. For example: RT MENTION Aankho ki hai ye khawaise ki
chehre se teri na hate … (id:36925, gold:positive, predicted:negative); RT
MENTION Saare jahan se achha # Hindustan hamara Ham bulbulain hai iski yeh
gulsitan hamara… (id:32707, gold:positive,predicted:neutral).

The system often failed to identify the sentiment, when a long tweet consists
of a number of complete or incomplete sentence with mixed sentiments. For
example: RT MENTION I miss childhood days … No problems … No hates … No shames
… No stress … No heartbreaks … Go school … Life was easy (id:25011,
gold:neutral predicted:positive); RT MENTION PM Modi won 356 seats & thanked
all 130 crore people for faith in democracy . Sonia won 52 seats & thanked
only the 12 cro (id:30674, gold:positive, predicted:neutral)

The proposed system often failed to identify the tweets that have any
punctuation or delimiter to separate the sentence clauses. For example: Jb Koi
aapke liye kam krta h aapke liye 5sal din rat ek krta h aap pe vishvas krta h
aapka vishvas jit ta h jis (id:31020, gold:neutral, predicted:positive)

We observed that mostly the system misclassified the positive and negative
tweets as neutral and vice versa. However, we observed that there are some
tweets that may arise some arguments, such as MENTION Beautiful words and its
one of my fav song (id:31662), MENTION Sir me from Bihar Love you You have
Great job for Bjp continue We are with you (id:39755), etc. tagged as neutral
and our system identified as positive. Identifying the polarity of an entity
or entity’s aspect in the tweets along with more training data could help to
resolve these issues.

##  6 Conclusion

This paper describes the approach we proposed for SemEval-2020 Task 9:
Sentiment Analysis for CM Social Media Text (SentiMix Hindi-English). In our
approach, we pre-processed the CM tweets and proposed a Recurrent
Convolutional Neural Network for the sentiment analysis of CM tweets. We
submitted two runs and obtaining promising results: our best run obtained
0.691 of F1 averaged across the positives, negatives and the neutral. We
observed that the proposed architecture occasionally strives to separate the
positive and negative polarities from the neutral and vice versa.

For future work, we will explore the performance of our model with larger
corpora against the testing set. Also, we would like to investigate other
embedding choices such as BERT [Devlin et al., 2019]. Moreover, due to the
impact that irony and sarcasm have on sentiment analysis [Hernández Farıas and
Rosso, 2016] it would be interesting to apply deep learning techniques to
detect irony [Zhang et al., 2019] but in a code-mixed scenario.

## Acknowledgements

The research work of the first four authors was supported by ERA-Net CHIST-ERA
LIHLITH Project funded by ANR (France) project ANR-17-CHR2-0001-03. The
research work of the last author was partially funded by the Spanish MICINN
under the project MISMIS-FAKEnHATE on Misinformation and Miscommunication in
social media: FAKE news and HATE speech (PGC2018-096212-B-C31).

## References

  * [Baccianella et al., 2010]  Stefano Baccianella, Andrea Esuli, and Fabrizio Sebastiani.  2010\.  Sentiwordnet 3.0: an enhanced lexical resource for sentiment analysis and opinion mining.  In LREC, volume 10, pages 2200–2204.
  * [Banerjee et al., 2020]  Somnath Banerjee, Monojit Choudhury, Kunal Chakma, Sudip Kumar Naskar, Amitava Das, Sivaji Bandyopadhyay, and Paolo Rosso.  2020\.  MSIR@FIRE: a comprehensive report from 2013 to 2016.  SN Computer Science, 1(1):1–15.
  * [Bojanowski et al., 2017]  Piotr Bojanowski, Edouard Grave, Armand Joulin, and Tomas Mikolov.  2017\.  Enriching word vectors with subword information.  Transactions of the Association for Computational Linguistics, 5:135–146.
  * [Choudhary et al., 2018]  Nurendra Choudhary, Rajat Singh, Ishita Bindlish, and Manish Shrivastava.  2018\.  Sentiment analysis of code-mixed languages leveraging resource rich languages.  arXiv preprint arXiv:1804.00806.
  * [Collobert et al., 2011]  Ronan Collobert, Jason Weston, Léon Bottou, Michael Karlen, Koray Kavukcuoglu, and Pavel Kuksa.  2011\.  Natural language processing (almost) from scratch.  Journal of Machine Learning Research, 12(Aug):2493–2537.
  * [Devlin et al., 2019]  Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova.  2019\.  BERT: Pre-training of deep bidirectional transformers for language understanding.  In Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers), pages 4171–4186, Minneapolis, Minnesota, June. Association for Computational Linguistics.
  * [Hernández Farıas and Rosso, 2016]  DI Hernández Farıas and Paolo Rosso.  2016\.  Irony, sarcasm, and sentiment analysis. chapter 7.  Sentiment Analysis in Social Networks. Morgan Kaufmann, pages 113–127.
  * [Joshi et al., 2016]  Aditya Joshi, Ameya Prabhu, Manish Shrivastava, and Vasudeva Varma.  2016\.  Towards sub-word level compositions for sentiment analysis of hindi-english code mixed text.  In Proceedings of COLING 2016, the 26th International Conference on Computational Linguistics: Technical Papers, pages 2482–2491.
  * [Lai et al., 2015]  Siwei Lai, Liheng Xu, Kang Liu, and Jun Zhao.  2015\.  Recurrent convolutional neural networks for text classification.  In Twenty-ninth AAAI Conference on Artificial Intelligence.
  * [Lal et al., 2019]  Yash Kumar Lal, Vaibhav Kumar, Mrinal Dhar, Manish Shrivastava, and Philipp Koehn.  2019\.  De-mixing sentiment from code-mixed text.  In Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics: Student Research Workshop, pages 371–377.
  * [Martínez-Cámara et al., 2014]  Eugenio Martínez-Cámara, M Teresa Martín-Valdivia, L Alfonso Urena-López, and A Rturo Montejo-Ráez.  2014\.  Sentiment analysis in twitter.  Natural Language Engineering, 20(1):1–28.
  * [Mejova et al., 2015]  Yelena Mejova, Ingmar Weber, and Michael W Macy.  2015\.  Twitter: a digital socioscope.  Cambridge University Press.
  * [Mohammad et al., 2013]  Saif Mohammad, Svetlana Kiritchenko, and Xiaodan Zhu.  2013\.  NRC-Canada: building the state-of-the-art in sentiment analysis of tweets.  In Second Joint Conference on Lexical and Computational Semantics (* SEM), Volume 2: Proceedings of the Seventh International Workshop on Semantic Evaluation (SemEval 2013), pages 321–327.
  * [Nakov et al., 2016]  Preslav Nakov, Alan Ritter, Sara Rosenthal, Fabrizio Sebastiani, and Veselin Stoyanov.  2016\.  SemEval-2016 task 4: sentiment analysis in twitter.  In Proceedings of the 10th International Workshop on Semantic Evaluation (SemEval-2016), pages 1–18, San Diego, California, June. Association for Computational Linguistics.
  * [Pak and Paroubek, 2010]  Alexander Pak and Patrick Paroubek.  2010\.  Twitter as a corpus for sentiment analysis and opinion mining.  In LREc, volume 10, pages 1320–1326.
  * [Patwa et al., 2020]  Parth Patwa, Gustavo Aguilar, Sudipta Kar, Suraj Pandey, Srinivas PYKL, Björn Gambäck, Tanmoy Chakraborty, Thamar Solorio, and Amitava Das.  2020\.  Semeval-2020 task 9: Overview of sentiment analysis of code-mixed tweets.  In Proceedings of the 14th International Workshop on Semantic Evaluation (SemEval-2020), Barcelona, Spain, December. Association for Computational Linguistics.
  * [Rosenthal et al., 2015]  Sara Rosenthal, Saif M Mohammad, Preslav Nakov, Alan Ritter, Svetlana Kiritchenko, and Veselin Stoyanov.  2015\.  SemEval-2015 task 10: Sentiment analysis in twitter.  In Proceedings of the 9th International Workshop on Semantic Evaluation (SemEval-2015), page 451–463. Association for Computational Linguistics.
  * [Sharma et al., 2015]  Shashank Sharma, PYKL Srinivas, and Rakesh Chandra Balabantaray.  2015\.  Text normalization of code mix and sentiment analysis.  In 2015 International Conference on Advances in Computing, Communications and Informatics (ICACCI), pages 1468–1473. IEEE.
  * [Zhang et al., 2019]  Shiwei Zhang, Xiuzhen Zhang, Jeffrey Chan, and Paolo Rosso.  2019\.  Irony detection via sentiment-based transfer learning.  Information Processing & Management, 56(5):1633–1644.