Asymptotic regimes in elastohydrodynamic and stochastic leveling on a viscous film

An elastic sheet that deforms near a solid substrate in a viscous fluid is a situation relevant to various dynamical processes in biology, geophysics and engineering. Here, we study the relaxation dynamics of an elastic plate resting on a thin viscous film that is supported by a solid substrate. By combining scaling analysis, numerical simulations and experiments, we identify asymptotic regimes for the elastohydrodynamic leveling of a surface perturbation of the form of a bump, when the flow is driven by either the elastic bending of the plate or thermal fluctuations. In both cases, two distinct regimes are identified when the bump height is either much larger or much smaller than the thickness of the pre-wetted viscous film. Our analysis reveals a distinct crossover between the similarity exponents with the ratio of the perturbation height to the film height.


I. INTRODUCTION
The motion of an elastic sheet supported by a thin layer of viscous fluid is a phenomenon that manifests itself in processes spanning wide ranges of time and length scales, from e.g. magmatic intrusion in the Earth's crust [1,2], to fracturing and crack formation in glaciers [3], to pumping in the digestive and arterial systems [4][5][6], or the construction of 2D crystals for electronic engineering [7]. Elastohydrodynamic flows have been studied in model geometries in order to understand their generic features and the inherent coupling between the driving force from the elastic deformations of the material and the viscous friction force resisting motion [8][9][10][11][12][13][14][15][16].
The investigation of an initially flat elastic membrane that is subsequently subjected to an applied deformation has helped disclose how system size, magnitude and direction of elastic deformations and spatial confinement affect the membrane dynamics. Whilst stretching [17,18] or compressing [19,20] of the membrane produces wrinkling patterns [21], it has been shown that slowly deforming the membrane, by means of injecting additional fluid into the supporting layer, leads to a dynamics where the fluid pressure is solely balanced by elastic bending forces [22,23]. As the out-of-plane deflection increases, a change in the physical mechanism that dictates the dynamics occurs as elastic stretching becomes the dominant driving mechanism until the system reaches a critical size for which gravity starts to dominate the dynamics. For small systems, an analog to such a predicted transition is observed for a thin perturbed elastic plate resting on a nanoscopic fluid layer, where the restoring elastic bending force is opposed by van der Waals forces leading to an elastohydrodynamic touchdown [24] similar to capillary film dewetting [25]. To describe the dynamics theoretically, one can solve the full Navier-Stokes equation in the fluid phase using dynamic boundary conditions at the elastic interface given by the solution of the Föppl-von-Kármán equation [26], using e.g. the immersed-boundary method [27]. Viscous flow in thin films can be described by the lubrication theory [28] that has been widely used to study different elastohydrodynamic flow phenomena [8,[22][23][24]26]. However, not much is known about how elastohydrodynamic flows are affected by the ratio between the geometric parameters that characterize the system as it undergoes large changes while the driving force remain the same.
For instance, when an elastic sheet deforms onto a wall pre-wetted by a thin viscous film, the dynamics of the advancing front are dictated by the local curvature of the interface [16,22]. This elastohydrodynamic relaxation is reminiscent of the capillary spreading of a viscous drop onto a solid substrate [29][30][31][32]. Similar to capillary flows, elastohydrodynamic relaxation processes are not only limited to very thin pre-wetted films. In fact, an elastic sheet with zero spontaneous curvature but with an initial shape of a bump (Fig. 1) with a height much larger than the pre-wetted viscous film will relax towards a flat equilibrium state. Inevitably, the system must then crossover from a situation where the bump height is larger than the pre-wetted film height, to a situation where instead the pre-wetted film becomes thicker than the bump. Here we investigate how the elastohydrodynamic leveling changes with the ratio between the bump height and the pre-wetted film thickness. In particular, are there different asymptotic regimes, and Initially, the overall profile presents a localized bump, whose profile is invariant in the y-direction, i.e., quasi two-dimensional. Far away from the perturbation, the viscous film has a constant thickness . In the bump region, the height profilesĥ(x, t) and h(x, t) =ĥ(x, t) − of the viscous film and the bump, respectively, vary with the horizontal position x and time t, and remain symmetric about x = 0. At x = 0, we define the characteristic height h(x = 0, t) = h0(t) of the bump and its typical radius R(t), with initial values given as h0(t = 0) = hi and R(t = 0) = Ri. how do the system transition from one to another? At the nanoscale, thermal fluctuations are expected to contribute and may even dominate the dynamics [33][34][35][36][37], which we quantify in the leveling dynamics. To answer these questions, we combine numerical solutions of a mathematical model based on the lubrication theory [28] with scaling analysis and experiments.

II. MATHEMATICAL MODEL AND NUMERICAL PROCEDURE
We consider the system depicted in Fig. 1, where we focus on a system where any influence of gravity can be neglected, i.e., the bump height is smaller than the elasto-gravity length [8]. Only situations where the bump height h(x, t) is small compared with its horizontal extent and where the film slopes are small, i.e., ∂ĥ(x, t)/∂x 1, are considered. We describe the viscous flow between the plate and the solid substrate using lubrication theory [28]. When the initial deflection h i of the elastic plate is small compared to its thickness d, we can neglect stretching and the pressure reduces to p(x, t) = B∂ 4ĥ (x, t)/∂x 4 , where B = Ed 3 /[12(1 − ν 2 )] is the bending rigidity of the plate, E is the Young modulus and ν is the Poisson's ratio [38]. In addition, the system is a spatially unconfined elastic sheet with the two lateral boundaries being free to move relative to the underlying fluid. Thus, the in-plane compression is suppressed, and bending stresses dominate the relaxation process regardless of the ratio d/h i . By assuming incompressible flow and imposing no-slip conditions at the two solid substrates, and considering a onedimensional geometry as there are no variations along the y-direction, one obtains the governing equation for the evolution of the height profile (see e.g. Ref. [8]) where µ is the fluid's dynamic viscosity. At small scales, thermal fluctuations can also influence the dynamics, which is described by the last term of Eq. (1). This term mimics the stress generated by thermal fluctuations, originates from an additional symmetric random stress term in the Navier-Stokes equations and is obtained by an integration in the z-direction (for details see [36,[39][40][41]). The noise term η(x, t), is multiplied by a prefactor Γ = k B T A /(6µb) where k B is the Boltzmann constant, T A is the ambient temperature, b is the width of the plate along the y-direction, and η(x, t) is modelled as a spatiotemporal Gaussian white noise such that η(x, t) = 0 and where the symbols indicate average quantities. We non-dimensionalize Eq. (1) When Γ = 0, this non-dimensionalization procedure gives us a parameter-free partial differential equation forĤ(X, T ). When Γ > 0, the non-dimensional number N = 2k B T A R 3 i /(Bh 2 i b) appears as a pre-factor in front of the stochastic term, and N 2 measures the ratio between thermal and bending energies. For the macroscopic system provided in our experiment and described in detail below, i.e., T A = 300 K, h i = 2.5 µm, R i = 20 µm, µ = 10 4 Pa s and B = 1.3 · 10 −12 Nm we get the noise prefactor Γ = 2.5 · 10 −13 ms −1/2 and the energy ratio N = 1.75 · 10 −6 which is well within the elastic bending dominated regime. However, a transition from a dominant elastohydrodynamic leveling to a dominant stochastic leveling would occur for a system with temperature T A = 300 K, membrane perturbation height h i = 10 nm and radius R i = 5 µm for a bending modulus B in the range of 10 − 100 k B T A where k B T A = 4 × 10 −21 Nm which corresponds to N in the range of 0 − 8 [36].
We solve the dimensionless version of Eq. (1) numerically by using a finite element method, and we split it into three coupled equations for: the bump profile H(X, T ) =Ĥ(X, T ) − /h i , the linearized curvature ∂ 2 H(X, T )/∂X 2 , and the bending pressure ∂ 4 H(X, T )/∂X 4 . These fields are discretized with linear elements and solved by using Newton's method from the FEniCS library [42]. For the deterministic case N = 0, an adaptive time stepping routine has been used with an upper time step limit of ∆T = 0.001 and a discretization in space ∆X ∈ [0.001; 0.01]. For the stochastic case N > 0, we have used a constant time step ∆T = 0.001, together with a discretization in space ∆X = 0.0025. At T = 0 we impose the initial condition: H(X, T = 0) = 1 − tanh(50X 2 ). We further impose the following boundary conditions at the boundary ∂Ω of the numerical domain: H X ∈ ∂Ω, T ) = H X ∈ ∂Ω, 0), ∂ 2 H(X ∈ ∂Ω, T )/∂X 2 = 0, and ∂ 4 H(X ∈ ∂Ω, T )/∂X 4 = 0. The noise Θ(X, T ) is introduced independently at each discrete position and time step using the "random" class with the "randn" Gaussian subclass from the Numpy library [43], with zero mean and a variance 1/(∆X∆T ). We avoid negative values ofĤ(X, T ) (that might occur in the stochastic case due to the fluctuations), by imposing that whenĤ(X, T ) < 10 −6 , it is put back to 10 −6 as in [33,35]. To verify the predictions of Eq.(1), we construct an experimental setup which is described in the following section.

III. EXPERIMENTAL PROCEDURE
The experimental setup is composed of a fiber of polystyrene (PS) with a glass-transition temperature T g, PS ≈ 100 • C deposited on a film of the same polymer supported on a silicon (Si) substrate. These samples are capped by a thin sheet of polysulfone (PSU) with T g, PSU ≈ 180 • C. Sample preparation is carried out as follows: PS fibers (with number-averaged molecular weight M n = 15.8 kg/mol, and polydispersity index PDI = 1.05, Polymer Source Inc., Canada) are pulled from the melt at 175 • C using a glass rod. Thin PS films are spin cast from a toluene solution onto 10×10 mm 2 Si substrates, leading to a thickness of 25 to 380 nm, measured using ellipsometry (Accurion, EP3). The films are annealed at 110 • C for at least 12 hours in vacuum to remove residual solvent and relax residual stresses. The PS fibers are then transferred onto the PS films and the ensemble is heated briefly above T g, PS . The heating allows the PS to flow, thereby resulting in a bump. Thin PSU films (M n ≈ 22 kg/mol, Sigma-Aldrich) are prepared by spin casting from a cyclohexanone solution onto freshly cleaved mica substrates (Ted Pella, USA). The PSU films have a thickness of ≈ 160 nm, measured using ellipsometry, and are annealed in vacuum at 200 • C for at least 12 hours. The PSU films are floated onto water and then transferred onto a supporting apparatus (described previously [44]), held only by the film edges. These freestanding films can be relaxed to an unstrained state ensuring no in-plane tension. The PSU films are finally transferred onto the PS sample. The part of PSU film at the edges of the Si wafer was then removed using a scalpel blade prior to annealing. This was done to ensure slippage at the boundary between the PSU film and liquid PS layer, thus rendering the relaxation bending dominated as discussed above.
After preparation, the samples were annealed on a hotstage (Linkam, UK) at 130 • C, which is above T g, PS but below T g, PSU . Hence, the PS becomes a viscous liquid while the capping PSU film remains an elastic solid, thus realizing the system illustrated in Fig. 1. The height profile is imaged during annealing using optical microscopy with a red laser line filter (λ = 632.8 nm, Newport, USA), which creates interference fringes in the region of the bump, as shown in Fig. 2(a), due to the light that is reflected from the Si substrate. It is clear from these fringes that the initial fiber and resulting flow are one-dimensional over length scales that are many times the width of the perturbation itself. Each interference fringe corresponds to a change in height of λ/(4n), where n ≈ 1.57 is the average index of refraction of the two polymers that make up the sample (n PS = 1.53 and n PSU =1.61). This allows the bump profile h(x, t) to be reconstructed by fitting a polynomial to the fringe data, as shown in Fig. 2(b). Such profiles can then be used to track the leveling dynamics, and to extract in particular the evolution of the height h 0 (t) of the bump with time for various initial geometries.

A. Elastohydrodynamic leveling
We first start by investigating the elastohydrodynamic leveling in the absence of thermal fluctuations (N = 0). In Fig. 3 we show the numerical solutions of the dimensionless version of Eq. (1) for N = 0 and we can see that the aspect ratio h i / controls both the time scale for leveling and the detailed features of height profile. The smaller h i / , the faster the dimensionless leveling process. Also, the dip created near the advancing front of the perturbation is enhanced both in magnitude and lateral extent for smaller h i / . We remark that for each initial aspect ratio there is a transition period of a few numerical time steps preceding the onset of the leveling process. This part of the data is not included in Fig. 4 as it is considered to depend on the initial condition, but does not influence the later dynamics. For h i / = 43, the numerical height profiles are further compared with our experiments, which are found to be in good agreement. We recall here that the elastic plate is floating on the liquid film and has edges that are free to move. Therefore, the pressure contribution from bending still largely dominates any contribution from stretching and Eq. (1) is still valid.
We now turn to a scaling analysis of Eq. (1) for N = 0. When h 0 (t)/ 1, the equation can be linearized and reduces to 12µ∂h/∂t = B 3 ∂ 6 h/∂x 6 and we deduce the long-term scaling for the temporal evolution of the horizontal length of the bump: R(t) ∼ [B 3 t/(12µ)] 1/6 . Since there is area conservation in the (x, z)-plane, we assume R(t)h 0 (t) to be constant, that is evaluated to R i h i at t = 0. By combining these scaling relations we get for h 0 (t)/ 1 where τ = 12µh 6 i R 6 i /(B 9 ) is the characteristic time scale for the bending-driven leveling dynamics. As we operate within the regime where bending dominates over stretching, a similar result is obtained by considering the force balance between the viscous and bending forces [19]. Also if we include isotropic stretching due to clamped boundaries a similar scaling law appear, but now with an additional logarithmic term, R(t) ∼ (t/log(t)) 1/6 [20]. However  2) and (3) are indicated with triangles. The inset provides a zoom in the region containing the experimental data for the three samples, with initial aspect ratios hi/ = 30, 43, 56; corresponding respectively to = 50, 67, 26 nm; hi = 1.52, 2.9, 1.48 µm; and Ri = 9.6, 16.5, 9.9 µm. The uncertainties in all experimental length scales are about 5%. To compensate for those, the characteristic time τ for each sample is multiplied by a free fitting factor α = 0.7, 0.13, and 1.3, respectively.
quasi-static solution to obtain the correct scaling [22], i.e., constant pressure in the bump, leading to [36] h 0 (t) ∼ τ t By balancing the two asymptotic predictions above, we expect the crossover between them to occur around t/τ ≈ 1.
In addition, these asymptotic regimes suggest that h 0 (t)/ is essentially a function of t/τ only, independent of the value of h i / in particular. In order to test our scaling predictions, we compute numerical solutions of the dimensionless version of Eq. (1) for N = 0, with h i / ∈ [10 −2 , 10 3 ] and extract h 0 (t)/ as a function of t/τ . These numerical results are plotted in Fig. 4 and compared to the experimental data. For each sample, the experimental data is matched to the numerical data through one fitting parameter α in front of the time scale τ . The values of the fluid viscosity and elastic Young modulus are highly sensitive to the temperature in the experiments, and we estimate them to be µ ≈ 10 4 Pa s [45,46] and E ≈ 2.6 GPa [47], respectively. Since all experiments were carried out at the same temperature and with the same polymer, sample-to-sample variations in τ result only from uncertainties in the geometrical parameters h i , R i , d, and . The obtained α values are 0.13, 0.7, and 1.3 for the three samples and each of these values are reasonably close to unity. More importantly, the sample-to-sample variations in α do not exceed a factor of 10, which is well within the expected relative error arising from the high sensitivity of τ to the geometrical parameters. The general agreement between the experimental data and the numerical predictions is good, over about 5 orders of magnitude in t/τ . The systematic early-time tail in the experimental data might be attributed to the initial compressive thermal stresses in the elastic layer, which arise due to the rapid heating of the samples from room temperature to T = 130 • C, which relax prior to leveling and the time needed for the initial shape to enter the asymptotic regime.
The master curve in Fig. 4 confirms that h 0 (t)/ is a function of t/τ . Furthermore, the two scaling regimes predicted above are indeed present, with pre-factors close to unity, and the crossover between the two being located near t/τ ≈ 1 as predicted. Any bump that initially starts in a thin pre-wetted film regime h 0 (t)/ 1 will eventually crossover to a thick-film regime h 0 (t)/ 1, with the corresponding power laws in time. As a final remark, a similar combination (not included here) of numerical simulations and scaling analysis can been performed for an axisymmetric geometry, leading to h 0 (t) ∼ t −2/11 for h i / 1, and h 0 (t) ∼ t −1/3 for h i / 1.

B. Stochastic leveling
Next we investigate the leveling process when it is dominated by thermal fluctuations (N > 0). As shown in Fig. 5, the numerical solutions suggest that the aspect ratio h i / is again essential, as it sets the time scale for leveling where the smaller h i / , the faster the dimensionless leveling process. Moreover, by comparison with the deterministic (N = 0) case in Fig. 3, the stochastic (N > 0) profiles exhibit spatiotemporal fluctuations and adopt different average shapes and leveling dynamics.
To go further, we propose a scaling analysis of Eq. (1), inspired by Ref. [31]. We consider specifically the N 1 limit, for which the thermal fluctuations are the dominant driving contribution to the dynamics and we assume that we can neglect the bending term so that Eq. (1) reduces to ∂h/∂t = Γ∂[( + h) 3/2 η]/∂x. We consider the average quantities h 0 (t) and R(t) , where we invoke the ∼ (tx) −1/2 scaling [33] for the root mean square value of the averaged noise over a space interval x and a time interval t. By assuming that the average area conservation in the (x, z)-plan can be expressed as is the characteristic time scale for the stochastic leveling dynamics. Interestingly, Eq. (4) describes a complete crossover between two asymptotic regimes in the stochastic leveling dynamics: for h 0 (t) / 1, we obtain h 0 (t) / ∼ (τ Γ /t) 1/4 , and thus we recover h 0 (t) ∼ t −1/4 [33]; while for h 0 (t) / 1, we get h 0 (t) / ∼ τ Γ /t. We expect the crossover between the two asymptotic regimes to occur around h 0 (t) / ≈ 1, i.e., around t/τ ≈ 1/8.
In order to test the prediction in Eq. (4), we compute the numerical solution of the dimensionless version of Eq. (1) for 5 ≤ N ≤ 8, with h i / ∈ [10 −1 , 10 2 ]. By averaging over a minimum of 30 realizations, we can extract h 0 (t) / as a function of t/τ Γ and the results are plotted in Fig. 6. The data from the numerical solutions is in good agreement with Eq. (4) for all h 0 (t) / , with no adjustable parameter. Our results highlight that Eq. (4) gives an accurate prediction of the stochastic leveling dynamics and show that the missing pre-factor is close to unity. Finally, in order to further highlight the underlying self-similarity associated with each of the two asymptotic regimes, the insets of Fig. 6 show the corresponding bump height profiles rescaled according to Eq. (4). In each asymptotic regime the height profiles collapse onto a universal shape which confirms the overall self-similarity in the leveling dynamics.

V. CONCLUSION
We have described the elastohydrodynamic and stochastic leveling of an elastic plate placed atop a viscous film. By combining numerical solutions, scaling analysis, and experiments, we identified various canonical regimes. Our results highlight the importance of the driving mechanism, either by elastic bending of the plate or thermal fluctuations, and the influence of the aspect ratio of the bump height to the film height. For each of these two driving mechanisms, a crossover between two distinct asymptotic regimes is controlled by the aspect ratio. These findings can be helpful to explain elastohydrodynamic leveling dynamics found in biological, engineering, or geological processes.