Instability, Rupture and Fluctuations in Thin Liquid Films: Theory and Computations

Thin liquid films are ubiquitous in natural phenomena and technological applications. They have been extensively studied via deterministic hydrodynamic equations, but thermal fluctuations often play a crucial role that needs to be understood. An example of this is dewetting, which involves the rupture of a thin liquid film and the formation of droplets. Such a process is thermally activated and requires fluctuations to be taken into account self-consistently. In this work we present an analytical and numerical study of a stochastic thin-film equation derived from first principles. Following a brief review of the derivation, we scrutinise the behaviour of the equation in the limit of perfectly correlated noise along the wall-normal direction, as opposed to the perfectly uncorrelated limit studied by Grün et al. (J Stat Phys 122(6):1261–1291, 2006). We also present a numerical scheme based on a spectral collocation method, which is then utilised to simulate the stochastic thin-film equation. This scheme seems to be very convenient for numerical studies of the stochastic thin-film equation, since it makes it easier to select the frequency modes of the noise (following the spirit of the long-wave approximation). With our numerical scheme we explore the fluctuating dynamics of the thin film and the behaviour of its free energy in the vicinity of rupture. Finally, we study the effect of the noise intensity on the rupture time, using a large number of sample paths as compared to previous studies.


Introduction
A thin liquid film can be understood as a layer of liquid with thickness typically ranging from fractions of a nanometer to several micrometers, typically resting or flowing on a substrate.Thin liquid films are quite often found in nature and nowadays technological applications: from gravity currents under water and lava flows, falling films of water on the surface of windows and sloppy streets on a rainy day, to films used in industrial coating processes for decorative, insulating or protective purposes and the cooling of microelectronic devices, to name but a few examples [8,29,34,47,51].Thin films have attracted considerable attention since the pioneering work of Reynolds on lubrication [58], who realised their significance in both applications and fluid dynamics fundamentals.Over the last few decades, understanding the behaviour of thin films and being able to make reliable predictions of their stability and dynamics has been crucial in the rapidly growing field of microfluidics, i.e. in the art of miniaturizing chemical and lab-on-chip devices (e.g.Ref. [55]).Such small-scale setups have shown a tremendous relevance to model and replicate biological systems, e.g.blood circulation systems [6,56,62], or biological processes, such as in-vivo protein crystallisation and bone formation [2,24,37].
It is therefore no surprise that the problem of modelling and predicting the behaviour of thin liquid films has attracted the attention of numerous researchers from engineering, physics, and applied mathematics, over the last few decades.In most practical cases, the thin films are subject to additional effects and complexities, such as body forces, often gravity as in the case of falling liquid films [14, 15, E S ij (x, y; t)S lm (x 0 , y 0 ; t 0 ) = 2 k B T µ q x (x x 0 )q y (y y 0 ) (t t 0 )( il jm ⇢(r; t) ⇢(r; t)u(r; t) e(r; t) ⇢(r; t)u(r; t) ⇢(r; t)(u(r; t) ⌦ u(r; t)) + p(r; t)1 (e(r; t) + p(r; t))u(r; t) + Q(r; t)u(r; t) ⌃(r; t)u(r; t) + ⌅(r; t) Fluctuating Hydrodynamics 0 @ ⇢(r; t) u(r; t) e(r; t) t) ê(r; t) Micro-and Mesoscopic Hydrodynamic Fields ṙi = @H @p i , ṗi = @H @q i Hamiltonian System Equivalent Stochastic Thin-Film Equation SPDE yielding the equation in Ref. [10], the most widely-accepted version of the stochastic thin-film equation (e.g.[20,43,55].This equation opens the door to numerical scrutiny of the effect of noise on dewetting. The present work introduces an alternative and more systematic derivation of the stochastic thinfilm equation from FH.Our derivation differs from the previous works fundamentally in two points.First, we do not require the noise to be delta-correlated in both the streamwise and the wall-normal direction, unlike the work of Grün et al.Indeed, we obtain the same SPDE by only imposing that the noise is perfectly correlated along the crossstream direction.We believe that this condition is more physically meaningful than imposing an uncorrelated noise along the crossstream direction, as is the case in previous works.Indeed, as we already emphasised, it is precisely the coherence of the film in this direction due to the action of viscosity that forms the basis of the long-wave approximation.The correlation coefficients along the streamwise direction are ultimately obtained by imposing the detailed-balance condition.Second, we present an efficient numerical scheme based on spectral collocation methods, which we believe to be more convenient as it gives us the opportunity of selecting in a straightforward way the frequency modes of the noise (following the spirit of the long-wave approximation).Along the way, we also perform the necessary theoretical analysis that is missing in the literature and provide a formal definition of the rupture time.This includes a detailed linear stability analysis of the initially uniform thin film, which is necessary to understand the origin and nature of the thin-film rupture.We show that the film is unconditionally stable to sufficiently small perturbations in the case of a negligible interface potential.In the case of a general interface potential, φ(h) = α h −3 − β h −2 , which is the sum of a non-negative convex and a concave term, we find that the condition required for the film to become linearly unstable is β > 2α.We also prove that the stochastic thin-film equation with the simpler multiplicative noise fulfils the detailed-balance condition, which is required from thermodynamical arguments [10,30], for a very particular noise structure.Moreover, we present simulations of the dynamics of rupture by using our numerical scheme and study how the noise intensity affects the rupture time.Our results and definition of rupture time seem to perfectly agree with the numerical results obtained by Grün et al.Our work thus brings closure, from the theoretical point of of view but also in the practical sense, to the stochastic thin-film problem and provides a general theoretical and numerical apparatus that will be relevant in a wide variety of natural processes which involve a similar stochastic gradient-flow structure.
In the same spirit, we believe that our work is of relevance not only for the thin-film community, for obvious reasons, but also to the current state-of-the-art in stochastic modelling and statistical physics.This is due to the fact that the theoretical analysis and the numerical method introduced here can be useful to research fields involving a similar SPDE, which can be found in a wide variety of problems.Examples of these are the study of nanojets via phenomenological equations [18,45], the dynamics of colloidal systems and/or systems with long-range interactions via a fluctuating dynamic density-functional theory (DDFT) [1,4,5,13,17,35,36], phase transitions in complex systems with coarse-grained models from fluctuating DDFT [40,41], crystallisation via stochastic phase-field-crystal (PFC) models [19,26,27], and many more.Indeed, in recent years there has been an explosion of interest in the formal study of PDEs and SPDEs of the same structure, ever since the pioneering work of Otto [48], where it was shown that the solutions of such kind of equations for a particular mobility tensor represent a gradient flow of the Dirichilet free energy with respect to the 2-Wasserstein metric (which metrizes the weak topology of probability measures [22,68]) when the problem is formulated into a discrete time variational scheme.Similar results were obtained by Lisini, Matthes, and Savaré [39] for concave mobility operators, but do not exist for the kind of mobility tensor involved in this work.More recently, Reina and Zimmer [57] found a general fluctuation-dissipation relation between the drift and the fluctuating term of a general stochastic gradient flow which also includes the stochastic thin-film equation considered here, so that both the maximum-entropy production (MEP) and the large-deviation (LD) principle are satisfied.Such a relation is indeed satisfied by the structure of the noise we obtain from imposing the detailed-balance condition to the equivalent multiplicative noise.It is therefore evident that theoretical and numerical methods presented in this work should be of interest to practitioners of many other fields, who might be interested in analysing the sensitivity to fluctuations or in simulating the dynamics of the system at hand, as already noted earlier.
In Sect. 2 we present the Landau-Lifshitz FH equations for a fluid film flowing on a horizontal substrate.We continue with the long-wave (or lubrication) approximation, for which we appropriately nondimensionalise the FH equations, and subsequently take a regular perturbation expansion for ε 1.
This yields the stochastic thin-film equation with a conservative noise involving a convolution integral.
A flow diagram to summarise our derivation and the relationship with previous approaches is given in Fig. 1.Before substituting the convoluted noise with a simpler multiplicative version, the gradient flow structure of the deterministic part of the equation is shown in Sect.3. In this section we also give the effective energy functional related to the thin-film equation and study the linear stability of an initially uniform thin film.In Sect. 4 a simpler multiplicative noise with free parameters, which is equivalent to the one derived from first principles, is proposed.We then study the implications of imposing the detailed-balance condition to the resultant SPDE and show the equivalence with the original noise in the limit of perfect correlation along the crossstream direction.In Sect. 5 we introduce a numerical method to integrate the stochastic thin-film equation based on a spectral collocation scheme.Finally, concluding remarks and discussion are offered in Sect.6.
2 Theoretical framework

Stochastic Navier-Stokes equation
The fluctuating dynamics of a two-dimensional thin film of Newtonian fluid flowing on a horizontal substrate can be described by the incompressible Navier-Stokes equations [30] (see also Fig. 1): where u = (u, v) is a two-dimensional velocity vector field, with u and v being the streamwise and crossstream components, respectively, µ is the dynamic viscosity, ρ is the fluid density, and p is the pressure of the fluid.The operator is the convective derivative, and S is the fluctuating stress tensor, which represents the effect of random thermal fluctuations on the film dynamics.The stress tensor S is symmetric and has zero mean, besides having the following correlation structure: where k B is the Boltzmann constant, T is a positive constant fixing the temperature/noise intensity of the system and the functions q x and q y are left undefined for the time being.At the wall, we apply the standard no-slip boundary condition: whereas at the fluid-air interface y = h(x; t), with h being the film height at the position x and time t, we apply the stress-balance boundary condition: where σ is the viscous stress tensor, κ is the mean curvature of the surface, γ is the surface tension coefficient, n is the normal vector to the interface, and Π = −φ (h) = − ∂φ(h) ∂h is the disjoining pressure, with φ(h) the interface potential, which models molecular interactions between liquid molecules and air [30].Additionally, we apply the following kinematic boundary condition at the interface: which states that a fluid particle on the interface will remain there for all times, thus preventing matter from leaving the interface via e.g.evaporation or any other mechanism.
Table 1 Nondimensionalisation of the FH equation, as proposed by Grün et al. [30].

Long-wave approximation: The stochastic thin film equation
In the following we simplify the FH equations by nondimensionalising them via the parameters shown in Tab. 1 following [30].This way we rewrite the equations in terms of two fundamental parameters, the characteristic film height d along the crossstream direction and the characteristic length scale λ along the streamwise direction.We then take the limit ε = d/λ 1 and retain terms up to and including O(ε).This step, and as already mentioned in the Introduction, is what is widely known as the long-wave approximation.In Tab. 1, U is the characteristic velocity scale of the flow and quantities with tildes are non-dimensional and taken to be of O( 1) with respect to ε.
The scaling of the deterministic terms follows e.g. from the analysis introduced by Dussan V. and Davis [67].The noise tensor is scaled in such a way that it is of the same order as the corresponding leading order terms in the viscous stress tensor.This ensures that the terms are retained to leading order in ε.For the Reynolds number, Re = ρU d/µ, we assume Re ∼ O(1).Thus, the equations for the streamwise and crossstream velocities are given by: where the tildes were removed for the sake of simplicity.To leading order in ε, this reduces to: On the other hand, the continuity equation (1a), no-slip condition (3) and the kinematic boundary condition (4) remain unchanged.It can also be shown that the curvature κ becomes ∂ 2 x h + O(ε 2 ), while the normal component of the normal stress tensor, n • (σ + S) • n, is simply −p + O(ε 2 ).Therefore, to leading order in ε, the interface stress balance condition (5) reduces to: The tangential component of the normal stress, given by t • (σ + S) • n, then becomes: which yields the following boundary condition: Integrating now Eq.(7a) with respect to y and applying the boundary conditions in Eqs. ( 3) and ( 10), we obtain the following expression for the streamwise velocity: Substituting Eq. ( 11) into the kinematic boundary condition Eq. ( 5), yields the time-evolution equation for the thin-film height: To simplify the fluctuating term, we can integrate the nested integral by parts and define S(r; t) = S xy (r; t), so that we finally obtain the central equation of this work: where the correlation structure of the noise is given by: Obviously, Eq. ( 13) shows exactly the same structure as the SPDE derived by Grün et al. [30], as we have closely followed their derivation.But at the same time there is a remarkable difference in that we are not imposing so far any specific structure for the correlators q x and q y , which in the study of Grün et al. [30] are imposed to be Dirac's delta.Indeed we will study later a different limit, where the noise is perfectly correlated along the crossstream direction, which as already noted in the Introduction it is more physically meaningful.
As also already emphasised in the Introduction, the main advantage of this approach to the thin-film dynamics is that the fluctuations are derived ab initio.That sets aside any controversy regarding the best noise to represent real fluctuations, as the noise is naturally derived from first principles within the self-consistent framework of FH.However, the fluctuating term in Eq. ( 13) is too convoluted for practical purposes and, also, against the spirit of the long-wave approximation to eliminate the dependence in the crossstream direction since it still shows explicitly a dependency in this direction.In Sec. 4 we will discuss how to derive an equivalent dynamics by proposing a much simpler version of the noise that is more convenient for analytical and numerical purposes.

Gradient flow structure
In this section we will demonstrate the gradient-flow structure of the deterministic part of the SPDE (13).To do that, we need to show that there exists an effective energy functional H(h) such that the deterministic part of the equation can be rewritten as the gradient of that functional with a certain metric.If we consider the functional: where L is the length of the domain in the x-direction, Eq. ( 13) can be rewritten as: with δ δh representing the functional derivative.Neglecting the noise for the time being we have that, which can be given a geometric interpretation as gradient flow on the energy surface in h-space with a metric of g −1 h = ∂ x M (h)∂ x , where M = h 3 /3 is the effective mobility, and DH := δH δh .Indeed, setting φ(h) = 0 for the time being, Eq. ( 17) belongs to a more general family of fourth-order degenerate PDEs of the form, which models the dynamics of different physical systems depending on the choice of the "mobility" M (h).For M (h) = h 3 + αh, the equation models the spread of a small viscous drop [28], and in the case of M (h) = |h| it models the behaviour of a thin neck of fluid in a Hele-Shaw cell [7] (see also the discussion in Ref. [9]).Similar equations are also obtained as the mean-field limits of interacting particle systems, with the resulting PDEs governing the evolution of the density of the system [1,4,5,13,17,35,36].Focusing only on the deterministic part of Eq. ( 16), i.e. without the noise term and setting φ(h) = 0, the resulting deterministic equation ( 17) is also conservative.This means that if the system has initially a given mass given by c = L 0 h(x; t = 0) dx, then the mass will remain constant over time.Moreover, for non-negative initial conditions, the positivity of the solution is preserved [3].
Let us now consider the time derivative of the energy functional which will give us information about the behaviour of the energy function over time.For this purpose, let us assume periodic boundary conditions in x, hence: with c being the mass of the initial condition, as mentioned above.The sign of the temporal derivative of the energy functional informs us about the dissipative nature of the functional for all acceptable film heights (i.e. with h ≥ 0 and mass with mass c) except for h * = c/L.Also, the functional H is bounded from below by 0 and is continuous.Thus, H satisfies the conditions required for it to be a Lyapunov functional [60].This implies that h * = c L is the only steady state for φ(h) = 0 and that it is globally attractive.
In the presence of a non-zero interaction potential φ, we can explicitly compute the first and second variations of H: where η belongs to a suitably defined space of admissible variations.It can be easily shown that any critical point of H, determined by δH(h, η) = 0, ∀η, is a stationary solution of the deterministic PDE (17).What is more, given a critical point of H we can solve a variational problem in η for the functional δ 2 H(h, η).Assuming η is sufficiently smooth, we obtain the following Euler-Lagrange equation, Multiplying by η both sides of Eq. ( 23), and after some manipulations, we obtain: Substituting this result into Eq.(22b) and assuming periodic boundary conditions in x, we get that δ 2 H = 0 for all critical points η (and, thus, by extension for all minimisers) of this functional.This results in δ 2 H(h, η) ≥ 0 for all critical points h of H and all admissible η.Hence, the critical points satisfy at least the necessary conditions for being minimisers of the functional H.

Linear stability analysis of the uniform state
In what follows we study the linear stability of an initial uniform state with constant height.To do that, we can linearise Eq. ( 17) about the state h * = 1 with mass c = L, by inserting perturbations of the form: Let us first consider a system with zero interaction potential, φ(h) = 0. Retaining only the terms of ( ) in the equation, the resulting linearised operator N can be written as: Considering again periodic boundary conditions, the eigenvalues λ n and eigenfunctions g n of the linearised operator are the following: for all n ∈ Z\{0}.As can be immediately seen, the eigenvalues λ n are strictly negative.As a result, the system is unconditionally stable to sufficiently small perturbations.
Let us consider now the interesting case of a non-zero interaction potential that can be written as the sum of a non-negative convex and a concave function of the height.Consider, e.g., the case of φ(h) = αh −3 − βh −2 with α, β ∈ R + .[The physical basis of this form goes back to the work by Derjaguin and Framkin, while using elements from the statistical mechanics of classical fluids, namely DFT, it can be shown that the (local) Derjaguin-Framkin disjoining pressure is an asymptote to that obtained from DFT as the distance of the chemical potential from its saturation value vanishes [70].]The linearised operator can then be written as follows, For x ∈ [0, ∞), we assume the factorisation g(x; t) = e ikx G(t), with k = (2π/L) ∈ R the wavenumber and L the wavelength.This yields Therefore, one branch of the neutral curve is given by k = 0 and the other is given by which defines the minimum wavelength (and, hence, domain length) for instability:  31b) and (32).Solid lines represent the values associated with the maximum growth rate, Eq. ( 33).Now, we can readily compute the wavenumber associated with the maximum growth rate, k max , and the corresponding domain length, L max , From (31b) it is straightforward that the condition for the system to be unstable is β > 2α.These results are summarised in Fig. 2.

Equivalent stochastic dynamics
Having characterised the linear stability conditions for the deterministic part of the stochastic thin-film equation derived in Sect.2, we turn to the effect of fluctuations.As repeatedly mentioned so far, the structure of the noise in Eq. (13) is not yet convenient for practical purposes.That is why there is a need for a statistically equivalent equation involving a simpler noise term.Particularly, Grün et al. [30] have shown that there exists an equivalence in law between the above following two SPDEs: where S and N are zero-mean spatiotemporal noise functions, defined by their correlation structure: Specifically, the equivalence is obtained by using a finite-difference discretisation of (34a) and, then, computing the related Kramers-Moyal coefficients to show that they both have associated the same Fokker-Planck equation (for delta-correlated noise in the x and y directions).In what follows, we use a finite-difference discretisation only to show that the simplified SPDE satisfies the detailed-balance condition (which is a requirement from thermodynamics [16]).This is done in Sect.4.1.After that, in Sect.4.2, we demonstrate that there exists an equivalency between the two SPDEs using a Q-Wiener representation of the noise function, S.

Detailed-balance condition
In [30] it was stated without proof that the simplified SPDE in Eq. (34b) satisfies the detailed-balance condition.In this section we present a proof of this result.Our proof outlined below is based on a spatial discretisation of a more general version of Eq. (34b) with the multiplicative noise given by νh n .
In the end, we find a general condition that the noise must fulfil based on the Fokker-Planck equation associated with the resultant set of SDEs.
Let us consider the SPDE: where ν, n ∈ R are free parameters for the moment, H being the energy functional given in Eq. ( 15).Let us also consider a finite difference discretisation of the x-axis such that x = a, with a ∈ R + the grid-spacing length and = 0, . . ., d − 1 the grid index.In that case, we have that: Using now the differentiation matrix for a central-difference scheme A ∈ R d×d , we can rewrite Eq. ( 36) as follows: where we have also introduced the discretised version of the drift, j = (j ) ∈ R d , and of the noise intensity, given by the diagonal matrix K = (K ) ∈ R d×d : We also define the diagonal matrix, K = (K ) ∈ R d×d as follows, In the above expressions, the discretised version H of the functional H is defined by the Riemman sum H = a l φ (h l ) − γ 2 (Ah) 2 l and ξ = (ξ 1 , . . ., ξ d ) is a Gaussian white noise with the following correlation structure: The adjoint of the infinitesimal generator, L † , of the Markov semigroup associated with this SDE is then given by: with Σ = AK(AK) T .Thus, the Fokker-Planck equation governing the temporal evolution of the probability density function (PDF) to observe a certain state at a given time, ρ(h; t), can be written as: with the initial condition ρ 0 = ρ(0).Here we have used the fact that A is antisymmetric and that ∇ h (AKK T A) = 0. Inserting the definition of the discretised drift given in Eq. (39a) into Eq.( 42), we obtain, If we now set ρ = ρ s to be the Gibbs measure associated with H (h), i.e., ρ s = 1 Z e − 1 T H (h) for some normalisation constant Z, we obtain: Now, we need to show that if K K T − KK T = 0, ∀h ∈ R d there exists and admissible, h ∈ R d for which Γ (h) = 0.It is sufficient to show this for only one component of the vector, [Γ h] l .To do this, we define the function, f (h ), as follows, Then, after some straightforward but lengthy calculations, we obtain an explicit expression for We can now pick h 0 = h 2 = h d−2 = h d−4 = κ/(2da), h 4 = 2h 2 for some 0 < κ < 1 and pick h 1 < 1/(ad) such that f (h 1 ) = 0.The other components of h can be chosen such that h a = 1.This gives us [Γ (h 0 )] = 0. Thus the flux J(ρ s ) can be zero, ∀h ∈ R d ,if and only if ν = ± 1 √ 3 and n = 3 2 .That is, the Fokker-Planck equation satisfies the detailed-balance condition and, thus, the generator L is symmetric [49], if and only if the general noise in Eq. ( 36) has exactly the same noise coefficient, ν and dependency on h, as proposed in Refs.[10,30].
This result is in agreement with the general fluctuation-dissipation relation for general stochastic gradient flows dz = g −1 (z)DH(z)dt + σ(z)dW(t), with K the metric, σ an operator acting on dW(t) and W(t) a Brownian sheet, recently discussed by Reina and Zimmer [57].Following these authors, for the MEP and the LD principles to be fulfilled concomitantly, the relationship g −1 ∝ σσ * must be satisfied.In our particular case, that means: whence, ν = ± 1 √ 3 and n = 3 2 if the proportionality constant is assumed to be unity.As can be seen, the detailed-balance condition is more restrictive as it imposes the proportionality constant without doubt.
In that case, the stationary process associated with the SDE in Eq. ( 38) is reversible, i.e. for every T ∈ [0, ∞), h(t) and the time-reversed process h(T − t) have the same finite-dimensional distribution [49].This means that, given any finite sequence of times t 0 < t 1 < t 2 < • • • < t k = T , and corresponding measurable subsets, A 0 , A 1 , A 2 , . . ., A k , the following identity is true: which means that, statistically, the stationary process is insensible to the time-arrow direction.

Representation of the noise in two dimensions
We turn now our attention back to Eq. (34a) with the aim of representing the spatiotemporal fluctuations, S, in a more convenient way, as we already mentioned before.When the noise term is finally represented as an infinite expansion in terms of independent real-valued Brownian motions, we will be able to impose the long correlation-length limit along the crossstream direction, i.e. l y → ∞, and show that the noise in Eq. (34a) (hence, the SPDE itself) converges to the much simpler one-dimensional noise structure of Eq. (34b).Consider a finite domain of length L along the streamwise direction x, and height h(x; t) along the wall-normal direction y, with periodic boundary conditions along x.Let H be the Hilbert space of all square integrable functions on [0, L] × [0, h].Assume that for every time t ∈ [0, ∞), S takes values in H.It is known [12,65] that a spatiotemporal noise process can be written as the formal time derivative of a Q-Wiener process, W, such that: where W is an infinite-dimensional zero-mean Gaussian process which takes values in H and is defined entirely by the covariance operator Q, which is symmetric, positive and of finite trace.Consider the following result [54]: where β k (t) are independent real-valued Brownian motions and where the convergence is in L 2 (Ω, F, P; H), the space of Bochner-square integrable functions, f : Ω → H.
Applying Prop. 1, we can obtain a formal representation for S as follows: where βk : Ω × [0, ∞) → R are independent white-noise processes, i.e. zero-mean Gaussian processes with correlation determined by E( βk (t) βl (s)) = δ kl δ(t − s), and λ k and g k are the eigenvalues and eigenvectors, respectively, of the operator Q, defined by its action on a field: We propose the functions q x and q y to be: where l x , l y are the correlation lengths, and Z x , Z y the normalisation constants, along x and y, respectively.We set Z x such that L/2 −L/2 q x dx = √ 2T , while Z y is left undefined for the time being.Q is symmetric and nonnegative, as expected.The trace is given by y dx dy and is, thus, finite.To close the alternative representation of the noise we need to solve the eigenvalue problem: The structure of the operator Q motivates us to write the eigenfunctions g k as the product of two functions g k (x, y) = X(x)Y (y), so that the eigenvalue problem is rewritten as: The solution to this eigenvalue problem has the following eigenfunctions: with the eigenvalues given by, At this stage, we require that b 0 must be constant and independent of l y .Thus, the normalisation constant Z y must be defined as, Finally, we have all the ingredients for the easier representation of W, and hence for S according to Eq. ( 52), together with the definitions given before in Eqs. ( 58)-( 60) and with β mn (t) as independent 1D Brownian motions.

Long correlation-length limit: Perfectly correlated noise along the wall-normal direction
As we mentioned at the beginning of Sect.4.2, we are going to show that to obtain Eq. ( 34b) what is left is to take the limit of W when the correlation length goes to infinity.Thus, Consider now the behaviour of the eigenvalues, b n , ∀n = 0, in the same limit: By dominated convergence we have, Now consider the Q-Wiener process W having the following representation, Fix t > 0, and consider the sequence of random variables, W ly (t) : Ω → H.Then, we have: for some which can be made arbitrarily small.Using this result and Eq. ( 64), we have W → W as l y → ∞, in L 2 (Ω; H).Thus, we can formally write that S ly → D = dW dt as l y → ∞.Finally, inserting this representation into (34a) one gets: Integrating the last term of the latter equation we eventually obtain the expression we were after: As can be seen, the infinite sum in this equation represents the spatiotemporal noise N of Eq. (34b).Thus, we can conclude that in the long correlation-length limit the SPDE (34a) converges to a b 0parametrised family of SPDEs with a simpler noise structure.We can now select a particular value of b 0 , which was left intentionally undefined before.The choice of that constant must be such that the resulting Fokker-Planck equation satisfies the detailed-balance condition with the invariant measure, ρ s , as specified in the previous section.Under such circumstances, b 0 should be equal to 4ν 2 = 4 3 and, therefore, the term multiplying the noise N will become ± h  [10] and Grün et al. [30].It is also interesting to note that numerical computation of the eigenvalues b n reveals that they are distributed according to a discrete Gaussian distribution of the form (see Fig. 3):

Numerical experiments
Having shown the equivalency between the two equations (34a) and (34b), we perform detailed numerical experiments using a spectral method as an alternative to the finite-element discretisation used in the previous study by Grün et al. [30].Spectral methods give us the opportunity of selecting, in a straightforward way, the frequency modes of the noise (following the spirit of the long-wave approximation), besides the convenience of smaller but denser matrices.SPDEs need to be discretised carefully, since different discretisations can converge to different limiting processes as the mesh size goes to zero.
A good discussion about this issue is offered in Ref. [31], where discrepancies are shown for different finite-difference approximations of the stochastic Burgers equation.
In what follows, we provide a brief description of how our numerical method works and, then, proceed to validate it by computing statistics for the rupture times.Considering the SPDE in (34b), we know that there exists a trace-class, nonnegative covariance operator Q such that its eigenfunctions can be used to provide a spectral decomposition of the noise term N .We argue that this set of eigenfunctions should form the natural setting to discretise the SPDE.

Description of the numerical method
Consider the uniformly spaced partition {x n } N −N ⊂ [0, L] with spacing ∆x.Again, we define the film height vector h ∈ R 2N +1 , h i = h(x i−N −1 ).Let X m and a m be the eigenfunctions and eigenvalues of Q respectively.Thus we can now define the spectral matrix, C ∈ R (2N +1)×(2N +1) such that C ij = X j−N −1 (x i−N −1 ) and the vector of eigenvalues, a ∈ R 2N +1 such that a i = a i−N −1 .This way, Cx for any x ∈ R 2N +1 is a Galerkin projection onto the space of eigenfunctions of Q .If we choose Q as follows: then X m and a m are exactly those in (58) and (59b).This considerably simplifies the differentiation operation.We can use this to define the spectral differentiation matrices D n ∈ R (2N +1)×(2N +1) that will be required for the numerical experiment.Let B, A ∈ R (2N +1)×(2N +1) be anti-diagonal and diagonal matrices, respectively, such that B i,2N +2−i = sign(N + 1 − i) and Then, we set D n = C −1 (AB) n C so that D n h amounts to an approximation of the n th -derivative of h.We can now write down a nonlinear SDE with multiplicative noise which approximates the solution of the SPDE: where: are the drift vector and the square root of the diffusion matrix, respectively, with y ∈ R 2N +1 defined as y i = 1.
Having the approximating SDE, all that is left is to choose an appropriate time-integration scheme.Since the deterministic problem is stiff, it is more appropriate and reliable to use two separate timesteps ∆t s and ∆t d < ∆t s , for the stochastic and deterministic parts of the SDE, respectively.Consider the uniform partition t n ⊂ [0, T ] with spacing ∆t s .Then, for each interval, [t i , t i +∆t s ] we integrate over time the deterministic equation, dh = b(h)dt with h ti as the initial condition using MATLAB's ode15s stiff explicit multistep solver to obtain h * ti+∆ts , where the solver selects ∆t d adaptively.We then add the stochastic component using the Milstein scheme, i.e. h ti+∆ts = h * ti+∆ts + √ ∆t s σξ+ ∆ts 2 σσ T (ξ 2 −1), where ξ is a vector of independent and identically distributed (iid) Gaussian random variables.

Simulations
Here we perform numerical simulations of the deterministic and the stochastic thin-film equations.We know from our analysis in Sect.3.2 that the uniform solution is linear unstable.This raises the question of whether or not there exist other stationary solutions to the thin-film equation.Choosing the domain length to be an integer multiple of L max and time-integrating with the initial condition h 0 = 1 + sin(k max x) for 1, the system converges to a second stationary solution, as is evident in Fig. 4. The second stationary solution has a cluster-like structure and corresponds to the film rupturing and breaking up into droplets.This behaviour resembles that seen in systems of interacting particles in statistical mechanics where the uniform distribution destabilises and gives rise to molecular clusters (see, e.g., Refs.[11] and [42]) but also in systems with state transitions induced by thermal fluctuations (e.g.Ref. [52]).
To ensure that there is no transient stationary solution between the uniform state and the droplet state, we have also computed the values of the free energy functional H as a function of time (see Fig. 5.2) and the functional derivative δH δh near the rupture time (see Fig 5 .2).As can be seen in Figs.5.2 and 5.2, the functional derivative is strictly non-zero between the two states, which indeed implies there is no transient intermediate solution.This calculation is completely necessary to know whether the system has to overcome an energy "barrier" during the breakup process, which then would comprise a nucleation event.
The inclusion of noise does not seem to induce a qualitative difference onto the behaviour of the thin film (see Fig. 6).The transitions from fluctuations about the uniform state to fluctuations about the droplet state can be seen in Fig. 6.To study the effect of noise on the thin-film dynamics we vary the value of the noise intensity and obtain statistics for the rupture time t r , which we define as: The reason behind this choice lies in the fact that the minimum film height changes rapidly in the vicinity of the film rupture time as can be seen in Fig. 7. Using this definition of the film rupture time, Fig. 4 Deterministic evolution of a thin film under the influence of the disjoining potential, φ(h) = h −3 30 − h −2 2 .The sequence must be understood in the following order (A)→(B)→(C)→(D).The uniform solution destabilises and then converges to a second stationary solution which has a cluster-like structure, i.e. the film ruptures and breaks up into droplets.we can then compute its dependence of the noise intensity.For each value of T we set the number of sample paths.The variation of tr = E(t r ) with √ 2T can be seen in Fig. 8. From the same figure is also evident the noise has the effect of reducing the time it takes on average for the film to breakup, as expected from physical intuition.Nevertheless, as the noise intensity increases, the rupture time tends to saturate to a fixed value.This saturation time seems to reflect the fact that the system has a proper time scale to react to fluctuations, no matter how intense they are.

Concluding remarks
In this work we have introduced an alternative and rigorous derivation of the stochastic thin-film equation from first principles.We have also presented numerical simulations to study the stability of an initially homogeneous state and the response of the system to fluctuations.The starting point of the derivation is the stochastic Navier-Stokes equations for the velocity field of an incompresible fluid on a planar horizontal solid substrate.The relation of these equations with an underlying Hamiltonian dynamics of the constituent particles of the fluid is sketched in Fig. 1, where we also highlight our  contribution to the state-of the art and summarise all possible model equations obtained from the original Hamiltonian system.Having the fluctuating equations governing the relevant hydrodynamic fields for the system dynamics, we apply the widely-known long-wave approximation, which makes possible a considerable reduction of the dynamics in terms of the height of the liquid film.Despite such a simplification, the resultant SPDE describing the temporal evolution of the height of the film along the streamwise direction contains a fluctuating term which is not convenient for practical purposes.
Inspired by the work of Davidovitch et al. [10] and Grün et al. [30], we propose a tractable but general state-dependent noise term to replace the original one.At this point, we analyse the condition the noise term must fulfil in order to satisfy the detailed-balance condition.We show that, for the new equation to fulfil the detailed-balance condition the noise must have the same structure as the one proposed in previous works [10,30].We subsequently justify the structure of the alternative noise term proposed by using a Q-Wiener representation of the original noise.We then show that the original stochastic dynamics and alternative SPDE converge when the long correlation-length limit is imposed along the wall-normal direction.That is, both noise terms produce an equivalent statistics when the simpler is considered to be perfectly correlated along the wall-normal direction.We believe this limit is more physically meaningful than the uncorrelated noise originally derived by Grün et al. [30].
We also demonstrate the gradient-flow structure of the thin-film equation and define the associated energy functional, H.By studying the variation of such a functional, we show that an initially spatially homogeneous film is unconditionally linearly stable to sufficiently small perturbations in the case of a negligible interface potential.In the case of a general interface potential, φ(h) = α h −3 − β h −2 , which is the sum of a non-negative convex and a concave term, we find that the condition required for the film to become unstable is β > 2α.This result is crucial in that it gives the conditions under which the dewetting process can occur.To scrutinse the nonlinear dynamics of the film and resulting rupture process, we propose a numerical algorithm based on a spectral collocation method.We perform simulations of the dynamics of the thin film and study the evolution of the energy functional close to the rupture of the film.Finally, we study the effect of the noise intensity on the rupture time, and our results are in agreement with those by Grün et al. [30].As a remark, we observe that the rupture time seems to saturate as the intensity increases, resembling the saturation of the escape time in some thermally-activated processes, e.g.nucleation.

4 )Fig. 1
Fig. 1 Flow diagram of the approach used to derive the stochastic thin-film equation from the full underlying Hamiltonian dynamics.Arrows indicate the interconnectedness of the different approaches.Thick black boxes/arrows: this work.Thin boxes/arrows: previous works.Text on arrows gives brief descriptions of the approximations/manipulations made.

Fig. 2
Fig.2Linear stability of the uniform state for the interaction potential φ = αh −3 − βh −2 with α = 1 30 and β = 1 2 as a function of the surface tension coefficient, γ.Panel (a) shows the stability map in terms of the k and γ values, whereas panel (b) shows the same on the L-γ plane.Shaded red regions correspond to the parameter values for which the system is unstable.Dashed lines represent the critical stability curve given by Eq. (31b) and(32).Solid lines represent the values associated with the maximum growth rate, Eq. (33).

Fig. 3
Fig. 3 Distribution of the eigenvalues bn of the covariance operator Q for different values of the correlation length along the wall-normal direction.

Fig. 5
Fig.5Evolution of the energy functional H as a function of time (left), showing two minima: the uniform and droplet states, and of its functional derivative δH δh near the rupture time (right), where the two white regions represent values of x and t for which the functional derivative is zero.

Fig. 7
Fig. 7 Single-path variation of the minimum film height near the rupture time with √ 2T = 0.01 and lx = 0.

Fig. 8
Fig. 8 Mean rupture time tr as a function of the noise intensity √ 2T .Error bars show the standard deviation of tr.