Defect Sensitivity in Instability and Dewetting of Thin Liquid Films: Two

Aug 12, 2006 - Rajesh Khanna , Narendra Agnihotri , Manish Vashishtha , Ashutosh Sharma , Prabhat Jaiswal , Sanjay Puri. Physical Review E 2010,011601...
0 downloads 0 Views 671KB Size
3108

Ind. Eng. Chem. Res. 2007, 46, 3108-3118

Defect Sensitivity in Instability and Dewetting of Thin Liquid Films: Two Regimes of Spinodal Dewetting Ruhi Verma and Ashutosh Sharma* Department of Chemical Engineering, Indian Institute of Technology, Kanpur-208 016, India

The instability, dynamics, and dewetting engendered by the van der Waals forces in a thin liquid film sandwiched between a solid substrate and bulk fluid phase are investigated based on linear stability analysis and large-area, nonlinear 2D and 3D simulations. The effects of initial free-surface heterogeneities on the length scale of instability, morphology, order of the resulting structures, and their dynamic evolution are examined. The simulations clearly show two distinct regimes of spinodal instability: (a) deep inside the spinodal territory (DIST) and (b) a defect-sensitive spinodal regime (DSSR). The latter regime with increased sensitivity to the initial conditions and local ordering and clustering of holes is obtained toward the periphery of the spinodal region where spinodal destabilization becomes weak, or the force per unit volume ∂φ/∂h f 0. Thus, relatively thin and thick films close to the spinodal boundaries that do rupture by the spinodal mechanism may give the appearance of nucleation-induced rupture. Finally, the effects of simulation domain size (or number of holes) on the morphological and dynamic characteristics of spinodal instability are studied to assess the minimal sample size of simulations required to faithfully mirror some important aspects of thin-film rupture. 1. Introduction Thin polymer films find extensive applications in industrial processing and manufacturing.1-19 In many of the coating applications, homogeneity, uniformity, and durability of the thin films are essential,3-7 but in soft lithography-based patterning applications, the film instability needs to be controlled and guided for engineering of self-organized mesostructures.9-19 As films become thinner, they become increasingly susceptible to defects resulting from the coating process, conjoining pressure,20-32 or substrate nonuniformities.33-37 Understanding the mechanisms leading to the instability of a thin film is, therefore, of increasing importance as the coatings become progressively thinner. Rupture in thin films can be initiated by (i) spinodal dewetting, in which the initial thermal or mechanical fluctuations in film thickness experience a driving force for growth at a certain lateral length scale and grow exponentially with time, leading to dewetting as soon as the amplitude reaches the size of the film thickness;1-8,20-32 and (ii) heterogeneous nucleation, in which dewetting occurs from defects or impurities in the film or substrate.33-37 At a fundamental level, spinodal dewetting is a spontaneous self-organizing process occurring because of a lowering of the system free energy resulting from a spontaneous growth of initial fluctuations or defects. On the other hand, nucleation requires that an energy barrier be overcome either by thermal fluctuations or by taking a lower-energy passage offered by some preexisting nucleation sites.33-40 The problem of thin-film rupture and spontaneous pattern formation during dewetting also holds promise in the creation of desired microand nanoscale patterns for a host of bulk-nanotechnology applications, for example, in optical coating, texturing of superhydrophobic surfaces, and organic and polymeric optoelectronics, including light-emitting displays.36 Although destabilization of thin ( 13 nm are dominated by heterogeneous nucleation, while those thinner than 13 nm are dominated by spinodal dewetting or a spontaneous mechanism. Seemann et al.41 studied the dewetting of thin low molecular weight PS film on silicon substrate, and based on their experimental results, they reported the presence of three different rupture mechanisms within one system, namely: (i) spinodal dewetting in 3.9 nm thick PS film on coated silicon wafer, (ii) thermal nucleation in 4.1 nm thick film, and (iii) heterogeneous nucleation in 6.6 nm thick PS film. Similarly, Becker and co-workers42 also demonstrated both experimentally and theoretically that a PS film of 3.9 nm thickness dewets the silicon wafer by spinodal mechanism, whereas a slightly thicker film (4.9 nm) produces an ordered structure. Thiele et al.45 have also investigated the presence of nucleation in a spinodal regime having altered structural signatures around a localized depression. Bolline et al.50 have reported results on dewetting of polystyrene films on silicon substrates in which the behavior resembles spinodal dewetting but does not strictly follow all the expectations of spinodal dewetting The main purpose of this paper is to theoretically characterize the hallmarks of spinodal instability and dewetting based on large-scale nonlinear simulations that take into account realistic levels of preexisting free-surface defects and heterogeneities. This should aid in the design and interpretations of thin-film experiments. In particular, we investigate in more detail an interesting aspect of spinodal dewetting where its morphological and dynamical aspects seem to depart from classical descriptions

10.1021/ie060615q CCC: $37.00 © 2007 American Chemical Society Published on Web 08/12/2006

Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007 3109

used to model a polymer film supported on a silicon substrate with an associated oxide layer (SiO2) and polymer adsorbed layer. 2.1. Governing Equation. The following equation derived from the long-wave approximation governs the stability and spatio-temporal evolution of a thin film subjected to the excess intermolecular interactions.6-8,19-22,27,30,33-36,38,44 Figure 1. Schematic diagram of a thin liquid film supported on a solid substrate. ho, dc1, and dc2 represent the mean film thickness, the coating thickness, and the adsorbed layer thickness, respectively. γ and µ denote the surface tension and the viscosity of the film, respectively.

Figure 2. Variation of spinodal parameter, (∂φ/∂h), with the film thickness. The inset shows a magnified view near the upper critical thickness where the film is spinodally stable, (∂φ/∂h) > 0.

and start to resemble nucleation-induced dewetting in some aspects.42,44-48 The following questions are addressed: (a) Are there subregimes of spinodal instability and dewetting with distinct behaviors? (b) What are the unique dynamical and morphological signatures to differentiate between the different pathways of dewetting in the spinodal regime? (c) How do the initial free-surface heterogeneities affect the dynamics and morphology of spinodally unstable films? (d) What is the role of domain size chosen in numerical simulations for capturing different essential complexities of the evolution of instability? The last issue is of importance because there are, as yet, no large-scale 3D simulations to assess the role of simulation size on different morphological and dynamical features of the instability. The above issues are addressed systematically by looking at the morphological response of different thickness films to random initial disturbances and localized small defects at the free surface of the film, as well as a combination of the two in simulation domains of different sizes. The largest domain considered contains ∼100 holes, which allows for a clear resolution of suprahole arrangements. 2. Theory Figure 1shows a schematic representation of a rather general system studied that can be used as a model for a variety of thin-film systems. The system consists of a thin liquid film supported on a coated substrate and a pseudocoating of a thin adsorbed layer. Here, we chose a high-energy substrate which has repulsive (stabilizing) interactions with the film for a thick film, a low-energy coating that tends to destabilize the film for a range of intermediate thickness, and a repulsive adsorbed layer that again stabilizes the film at very low thickness. As discussed in Section 2.2, the variation of energy (per unit volume) for such a system has the qualitative form shown in Figure 2, where films in an intermediate range of thickness hc1-hc2 are spinodally unstable. For example, the system chosen can also be

3µ(∂h/∂t) + ∇‚[γh3∇(∇2h)] - ∇‚[h3∇φ] ) 0

(1)

where h(x,y,t) is the local thickness as a function of time, t, and spatial coordinates, x and y; φ is the conjoining pressure (body force per unit volume); and γ and µ denote the surface tension and viscosity of the film, respectively. The above equation represents the sum of three forces: (i) unsteady force/pseudoviscous force (first term), (ii) surface tension force due to local curvature of the free surface (second term), and (iii) the sum of excess intermolecular forces (third term). The unsteady force, which includes viscous effects, controls the growth rate of instability. The surface tension effect may be stabilizing or destabilizing. Finally, a gradient of conjoining pressure, ∇φ, caused by local changes in the film thickness engenders flow of fluid in the direction of lower free energy. 2.2. Excess Intermolecular Interactions. A general potential for an unstable film is made up of antagonistic (attractive and repulsive) interactions with different rates of decay with distance. In the absence of attractive forces, surface instability is not possible, whereas no study of dewetting/contact angle is possible unless short-range repulsion near the surface is present. The excess free energy per unit area, ∆G, and the conjoining pressure, φ ()∂(∆G)/∂h), for a thin film on a substrate with a thin coating and an adsorbed polymer layer (Figure 1) are obtained by a pairwise summation of various van der Waals interactions: (i) interactions among the molecules of the film, (ii) interactions among the molecules of the film and the molecules of coating, (iii) interactions among the molecules of the film and those of the substrate, (iv) interactions among the molecules of the film and the molecules of the adsorbed layer, (v) interactions among the molecules of the coating and the molecules of the overlying bulk fluid, (vi) interactions among the molecules of the substrate and the molecules of the adsorbed layer, (vii) interactions among the molecules of the adsorbed layer and the molecules of the overlying bulk fluid, and (viii) interactions among the molecules of the film and the molecules of the overlying bulk fluid.28 The total energy per unit area is given by

12π∆G ) [-(Aff + Aab - Afb - Aaf)/h2 (Aaf + Acb - Acf - Aab)/(h + dc2)2 (Acf + Asb - Acb - Asf)/(h + dc1 + dc2)2] (2) where Aij refers to the Hamaker constant for the van der Waals interactions between molecules of type i and j. The subscripts f, a, c, s and b denote the film, the adsorbed layer, the coating, the substrate, and the bounding fluid, respectively. In terms of effective Hamaker constants, Ac, Aa, and As,

-12π∆G ) [(Aa)/h2 + (Ac - Aa)/(h + dc2)2 + (As - Ac)/(h + dc1 + dc2)2] (3) where the three different effective Hamaker constants in the above are defined as Ai ) Aff + Aib - Afb - Aif, where the subscript i refers to phases a, c, or s.

3110

Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007

∆G can also be expressed in terms of the macroscopic spreading coefficients as

∆G(h) ) d02[Sa/h2 + (Sc - Sa)/(h + dc2)2 + (Ss - Sc)/(h + dc1 + dc2)2] (4) where Si ) -(Ai/12πd02) denotes the spreading coefficient of the film fluid on phase i. In the context of thin films considered here, the excess free energy (∆G) is a combination of the following: (a) Ss > 0, signifying a long-range van der Waals repulsion between a highenergy wettable substrate (e.g., silicon) and the film (e.g., polystyrene), (b) Sc < 0, signifying an intermediate-range attraction (e.g., between a low-energy, nonwettable thin SiO2 coating of the substrate and the polystyrene film), and (c) Sa > 0, signifying a short-range repulsion close to the solid surface, for example, because of the chain adsorption or grafting, with dc2 being an effective adsorbed-layer thickness. Figure 2 shows the qualitative behavior of the spinodal parameter, ∂φ/∂h, for this generic potential. The system is unstable for film thicknesses where ∂φ/∂h < 0, namely, between the two critical thicknesses hc1 and hc2. The spinodal parameter ∂φ/∂h ) ∂2∆G/∂h2 is obtained as (based on eq 4)

∂φ/∂h ) 6d02 [Sa/h4 + (Sc - Sa)/(h + dc2)4 + (Ss - Sc)/(h + dc1 + dc2)4] (5) As is shown below, based on the pathways of dewetting and morphology, the spinodally unstable region (∂φ/∂h < 0) can be divided into two different regimes: (a) deep inside the spinodal territory and (b) the defect-sensitive spinodal regime. The first regime is obtained away from the two critical thicknesses hc1 and hc2. The defect-sensitive spinodal regime (DSSR) is manifested toward the periphery of the spinodal region where ∂φ/∂h f 0, namely, near to the critical thicknesses. From simulations it is shown below that, as films gets closer to hc2 (∂φ/∂h f 0, in the DSSR), the dewetted structures become increasingly more sensitive to the initial conditions and produce the appearance of nucleated clusters of holes. 2.3. Linear Stability Analysis. Linear stability analysis of eq 1 gives the dominant length scale and time scale of instability. Equation 1 is linearized by h ) ho +  exp(ωt) sin (kx), where k is the wavenumber, ω is the growth coefficient, and  is the amplitude (,ho) of the initial disturbance. This leads to the following linear dispersion relation: ω ) -(ho3k2(φh0 + γk2))/ 3µ, where φh0 ) ∂φ/∂h is evaluated at the mean thickness, h ) ho. The spinodal dominant wavelength (λm) of the fastestgrowing linear mode (corresponding to dω/dk ) 0) and the shortest time of rupture (tr) are given by the following: λm ) (-8π2γ /(∂2∆G/∂h2))1/2 and tr )12µγ[ho3((φh)o)2]-1 ln(ho/), respectively. Parts a and b of Figure 3 show the variation of the spinodal wavelength (λm) and the time of rupture (tr), respectively, with film thickness in the region of spinodal instability. It is evident from these figures that both the quantities (λm and tr) increase tremendously as the upper critical thickness is approached. 2.4. Nondimensional Equation. The following scales were used to nondimensionalize the governing equations and the boundary conditions for a compact representation and numerical solutions: H ) h/ho, X ) x/((γ/6|Ss|d02)1/2)(ho2)), T ) t/(µγho5/ 12Ss2d04), and Φ ) φ/((6|Ss|d02)/ho3). The nondimensional thin-film equation, obtained from its dimensional counterpart (eq 1), governs the stability and spatio-

Figure 3. Variations of (a) the spinodal wavelength and (b) the rupture time, with the film thickness in the defect-sensitive spinodal regime. The critical thicknesses for the system are hc1 ) 1.83 nm and hc2 ) 6.81 nm. The parameters used for the simulations are as follows: Ss ) 20.0 mJ/m2, Sc ) -12.0 mJ/m2, Sa ) 2.0 mJ/m2, dc1 ) 2.5 nm, dc2 ) 1.0 nm, µ ) 1 Pa s, and γ ) 38.0 mJ/m2

temporal evolution of the film subjected to the excess intermolecular interactions.6-8,19-22,27,30,33-36,38,44

∂H/∂T + ∇‚[H3∇(∇2H)] - ∇‚[H3∇Φ] ) 0

(6)

2.5. Numerical Methods. The above nondimensional eq 6 is discretized in space using central difference technique combined with half node interpolation. The resultant set of stiff coupled ordinary equations (ODEs) was solved as an initial value problem with a volume preserving initial random perturbation, as well as other forms of initial perturbations, including localized defects on the film surface. These equations were integrated in time by using the Gear’s algorithm for stiff equations. Periodic boundary conditions were used. Numerical accuracy and convergence were checked by varying the number of grid points. 3. Results and Discussions Figure 3 depicts the predictions of the linear analysis regarding the dominant wavelength, λm, and the linear time of rupture of the film by the dominant wave. Both of these diverge as the film thickness approaches the stable boundary, hc2. Further, the linear theory predicts the formation of one hole or droplet in an area of λm2, regardless of the details of the initial disturbances. The effect of initial surface roughness (determined by the preparation conditions and materials), on the evolution of liquid-air interface, morphology, number density of structures formed, and the time taken for rupture, is critically examined by nonlinear simulations. From simulations, it is found that, for relatively thicker films in the “defect-sensitive spinodal regime, DSSR” defined by ho f hc2; ∂φ/∂h f 0, the morphology of dewetted structures formed crucially depends on the details of the initial surface heterogeneities. The effects of initial conditions are studied by considering the following three different cases of initial surface heterogeneities at the free surface of the film: (i) a random perturbation of the free surface over the entire simulation domain, (ii) a localized small depression on the free surface surrounded by a random perturbation, and (iii) periodic perturbations of the free surface. Parameters such as the mean film thickness and the amplitude of initial free-surface heterogeneity, localized defects, and random undulations are varied to observe the change in the film morphology. In this paper, the term “surface heterogeneity” will be used to denote the localized or randomly distributed preexisting defects on the free surface of the film or departures that occur from a flat free surface during the film preparation.

Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007 3111

Figure 4. Evolution of a 2.0 nm thick film subjected to a random perturbation of small amplitude (1% of film thickness). The length of simulation domain is 3λm, where λm is the characteristic length scale from the linear theory. Continuous shading between the maximum and minimum thickness in each image is used. Black represents the maximum thickness, and white represents the minimum thickness in all the figures. Images 1-6 are at times T ) 6.83 × 104, 1.47 × 105, 1.81 × 105, 2.25 × 105, 5.07 × 105, and 6.73 × 105, respectively. The parameters used are as follows: Ss ) 20.0 mJ/m2, Sc ) -12.0 mJ/m2, Sa ) 2.0 mJ/m2, dc1 ) 2.5 nm, dc2 ) 1.0 nm, µ ) 1 Pa s, and γ ) 38.0 mJ/m2.

Unless otherwise specified, the parameter values used for the simulations are (as in Figure 2) as follows: Ss ) 20.0 mJ/m2, Sc ) -12.0 mJ/m2, Sa ) 2.0 mJ/m2, dc1 ) 2.5 nm, dc2 ) 1.0 nm, µ ) 1 Pa s, and γ ) 38.0 mJ/m2. 3.1. Evolution of Random Surface Heterogeneities. Figures 4-9 depict the various pathways of dewetting and morphological structures formed in thin liquid films of different thickness during the growth of initially random free-surface heterogeneities. Previous studies19,24,27,30,32-36 have predicted that pathways of dewetting and the various morphological structures formed crucially depend on the location of the mean film thickness ho vis-a`-vis its deviation from the location of minima hm in the spinodal parameter (φh) vs thickness (h) curve (Figure 2). Figure 4 depicts the evolution of a relatively thin 2.0 nm film, which is thinner than the location of the minimum in the φh vs thickness curve (ho < hm, Figure 2). Initial evolution thus indeed starts by the reorganization of random surface heterogeneity of small amplitudes (1% of initial film thickness) on the length scale predicted by the linear theory (image 1 of Figure 4). Further evolution forms a bicontinuous undulating structure consisting of hills and valleys (image 2 of Figure 4), which directly disintegrates to form liquid droplets without first forming distinct circular holes. The number density of drops is identical to the predictions made by linear stability analysis (images 3 and 4 of Figure 4). Later, the smaller droplets start to merge into bigger droplets because of the Laplace and disjoining pressures induced flows from the smaller to bigger drops (images 5 and 6). The final equilibrium morphology (not shown) consists of a single drop surrounded by an equilibrium flat film. This simulation depicts that the liquid films on the left of the minimum in the spinodal parameter evolve by the formation of liquid droplets which further coalesce and finally form a single droplet. Figure 5 shows the evolution of a slightly thicker film (ho > hm), which is to the right of the minimum in the spinodal curve, subjected to an initial random surface heterogeneity of very small amplitude (1% of initial film thickness). Simulations are shown for two different domain sizes. The dewetting now initiates by the formation of holes (images 2 and 3 of Figure 5a) rather than by drops. The initial holes grow (image 3) and coalesce (image 4) to produce long liquid threads (image 5), which in turn fragment because of the Rayleigh instability and

Figure 5. Evolution of a 4.0 nm thick film subjected to random perturbation of small amplitude (1% of initial film thickness). (a) The length of simulation domain is 3λm, and images 1-6 are at times T ) 4.46 × 103, 6.50 × 103, 6.96 × 103, 7.99 × 103, 1.10 × 104, and 3.07 × 104, respectively. (b) Length of simulation domain is 10λm, and images 1-6 are at times T ) 6.84 × 103, 7.49 × 103, 8.50 × 103, 9.43 × 103, 1.48 × 104, and 3.87 × 104, respectively.

form droplets (image 6). The maximum number density of holes again closely matches with the linear stability analysis (8 holes on 3λm × 3λm simulation domain and 92 holes in 10λm × 10λm). Figure 5b represents the same simulation but with larger simulation domain length (10λm × 10λm) to gain a better understanding of the distribution of the structures formed. The large-scale simulation shows a rather even distribution of holes over the entire area and similarity of hole sizes (image 3 of Figure 5). The large scale simulation also better illustrates the late time ridge structure resulting from the hole coalescence (image 5). As shown below for thicker films, the resolution of larger structures on the scale of clusters of several holes becomes increasingly more important for statistical characterization. Figure 6 illustrates the 2D simulation results of dewetting for thicker films, ho > hm, close to critical thickness, hc2 ()6.81 nm), where ∂φ/∂h f 0. Figure 6a depicts the time evolution of instability of a 5.5 nm (>hm) thick film on 10λm domain, subjected to initial random surface heterogeneity of very small amplitude (1% of initial film thickness). The random surface heterogeneity rapidly reorganizes close to the characteristic length scale of instability (profile 1 of Figure 6a), and eventually, the film dewets the substrate. The maximum number density of holes (9 on 10λm domain, profile 2 of Figure 6a) is remarkably close to the linear theory predictions. However, as the film thickness is further increased, the number density of full thickness holes decreases, and for films close to the critical (h g 6.3 nm), where ∂φ/∂h is vanishingly small, the number density of holes decreases rather drastically (parts b and c of Figure 6). Thus, for the films in very close vicinity of the metastable or spinodally stable region, the number density of

3112

Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007

Figure 6. Evolution of instability in films of different thicknesses when subjected to a random perturbation of small amplitude (1% of initial film thickness) over a domain of 10λm. Profile 1 is at the instant of rupture, and profile 2 is at the time when the maximum number of holes have formed. (a) Time evolution for 5.5 nm thick film; profiles 1 and 2 are at times T ) 14 899 and 19 027, respectively, and the maximum number of holes formed is Nmax ) 9. (b) For a 6.5 nm thick film, profiles 1 and 2 are at times T ) 159 312 and 213 051, respectively; the maximum number of holes formed is Nmax ) 4. (c) For a 6.6 nm thick film, profiles 1 and 2 are at times T ) 292 550 and 277 930 365, respectively; the maximum number of holes formed is Nmax ) 3.

holes may gradually lose correlation with the spinodal scale. In this case, divergence of the linear time scale makes the evolution very sensitive to the local amplitude of the heterogeneity. The few holes that are able to form and grow initially in this case suppress the formation of other potential holes. Thus, based on the morphology and the number density of dewetting features, one can divide the instability zone into two regimes, namely: (i) deep inside the spinodal territory or DIST (ho e 5 nm in this example) and (ii) defect-sensitive spinodal regime, DSSR (ho f hc2 ) 6.8 nm). In the deep inside spinodal territory regime, the predictions of linear theory (λm ≈ h2 and tr ≈ h5) are in complete agreement with the direct numerical result of thin-film evolution equation (eq 1). However, in the DSSR, even extremely small variations of the initial heterogeneity structure are greatly amplified to produce a faster rupture by a smaller number of holes initially, which grow to suppress the formation of other potential holes. Thus, even though the length scale of the initial short-time (small-amplitude) evolution may be similar to the linear theory prediction both in the DIST and DSSR, the final number of fully grown holes are drastically different. This issue is important in the interpretation of thin-film experiments where often only the more accessible hole-formation stage is resolved. A clear differentiation of the nucleative dewetting from the DSSR should require observations of both the smallamplitude initial growth and the late holes. The spinodal correlations visible in the very early stages may be completely lost later in the spacings and number count of fully grown holes, thus mimicking the nucleation-induced dewetting. Another interesting morphological aspect observed in thicker films (in the DSSR) is a wider distribution of hole sizes, because not all holes appear at about the same time. These morphological

Figure 7. Evolution of a 5.5 nm thick film subjected to a random perturbation of small amplitude (1% of initial film thickness). (a) The length of simulation domain is 3λm, and images 1-6 are at times T ) 4.35 × 103, 1.58 × 104, 1.69 × 104, 1.83 × 104, 2.09 × 104, and 3.10 × 104, respectively. (b) The length of simulation domain is 5λm, and images 1-4 are at times T ) 1.89 × 104, 2.05 × 104, 2.17 × 104, and 2.64 × 104, respectively. (c) The length of simulation domain is 10λm, and images 1-4 are at times T ) 1.91 × 104, 2.10 × 104, 2.24 × 104, and 2.54 × 104, respectively.

aspects are best illustrated with the help of large-domain 3D simulations as discussed in Figures 7 and 9. Figure 7 depicts the 3D evolution of a thicker 5.5 nm thick film, which is toward the defect-sensitive spinodal regime. Parts a, b, and c of Figure 7 present the evolution of morphology in simulation domains of area 9 λm2, 25 λm2, and 100 λm2, respectively. Dewetting is initiated by the formation of a few widely separated primary holes (image 1 of Figure 7), which develop rims during their growth (image 2 of Figure 7). Interestingly, the second generation of holes clusters around the initial holes, rather than forming evenly in the empty spaces (image 2). The early holes continue to grow even as the new ones appear in their surroundings. Finally, the holes coalesce

Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007 3113

Figure 8. (a) Kinetics of hole growth and (b) kinetics of increase in dewetted area for different simulation domain lengths for a 5.5 nm thick film.

Figure 9. Evolution of a 6.0 nm thick film subjected to a random perturbation of small amplitude (1% of initial film thickness). The length of simulation domain is 10λm, and images 1-4 are at times T ) 4.11 × 104, 4.27 × 104, 4.43 × 104, and 4.53 × 104, respectively.

to form liquid ridges that fragment (image 4) and then ripen at long times (image 6 of Figure 7a). The morphological evolution of this intermediate case has signatures of both the regimes, in that the length scale of instability and the maximum number density of holes are quite close to the linear theory predictions as shown in Figure 8a, but the structure evolves by the formation of a small number of primary holes around which are seeded a large number of satellite holes that grow continuously to fill the spaces available. Further, even small-area simulations can adequately resolve the fractional dewetted area as shown in Figure 8b. However, the small-area simulations (Figure 7a) do not clearly resolve the seeded-cluster structures, which are depicted in the large-area simulation of Figure 7c and Figure 9 for a thicker film close to the spinodal boundary in the DSSR. Figure 9 shows the evolution of 6.0 nm (fhc2 ) 6.81) thick film, where dewetting is initiated by a smaller number of primary holes that grow in size and, later, smaller satellite holes start to open up around the primary holes (images 1 and 2 of Figure 9). The number density of holes is found to be much smaller than as compared to the linear theory predictions. This is because the initially formed primary holes (which evolve much before the appearance of secondary holes) grow and occupy a large fraction of available area and a very small space is left behind for other depressions to open up. A wider distribution of hole sizes becomes even more pronounced compared to the case of the 5.5 nm thick film. The size distribution of holes and other statistics in various stages of evolution for different thickness films could be extracted from the large-area simulations of 100 λm2, where the linear theory predicts 100 holes. These simulation results are summarized in Figures 10A and 10B. For thinner films, the

Figure 10. Size distribution of holes, number of holes (N) vs hole diameter (D), at an intermediate stage of dewetting when the dewetted area is ∼(13 ( 2)%. Parts a-c correspond to film thicknesses 4, 5.5, and 6 nm, respectively. Figure 10B. Size distribution of holes at a later stage of instability when nearly maximum number of holes are present. Parts a-c correspond to film thicknesses 4, 5.5, and 6 nm, respectively.

maximum number of holes are close to the linear theory expectations (Figure 10B (part a)), and a majority of them form within a narrow window of time to give a near-normal distribution of hole sizes at an early stage of dewetting (Figure 10A (part a)). Increased film thickness (Figures 10A (part b) and 10B (part b)) results in a delayed formation of greater proportion of holes leading to a wider distribution of hole sizes and a propensity toward bimodality. The distribution becomes almost uniform for the thickest film investigated in the DSSR (Figures 10A (part c) and 10B (part c)), indicating that newer holes continue to appear even as the growth of early holes continues. Compared to the linear theory predictions, fewer holes in thicker films also translate into a larger mean diameter of the holes (in terms of λm) and larger fractional dewetted area. Thus, a wide and more uniform hole size distribution together with clustering of holes are hallmarks of the DSSR. It may be noted that these important structural phenomena are not statistically resolved in small-scale simulations. 3.2. Localized Depressions Combined with Random Surface Heterogeneity. During the film preparation, small localized defects on the background of statistically small amplitude heterogeneities can result. In this section, the effects of a localized depression of slightly higher amplitude surrounded by random heterogeneities of comparatively smaller amplitude are investigated. Figures 11-18 depict the structural evolution of the same films as in the previous section, but with slightly different initial surface profile. Figure 11 shows that evolution of a thinner film deep inside spinodal territory proceeds with a

3114

Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007

Figure 11. Evolution of a 4.0 nm thick film starting with a localized depression of 2% amplitude and a 1% random roughness on its surroundings, over a domain size of 10λm. The nondimensional times of evolution for profiles 1-5 are T ) 0, 3.31 × 103, 5.64 × 103, 2.53 × 105, and 5.01 × 107, respectively.

cascade of holes that grows like a ripple progressively around the initial defect (profiles 1-3). Interestingly, the hole spacing still respects the linear theory prediction of λm. Thus, at a later stage when the maximum number of holes have formed, the signatures of the initial defects are obliterated (profile 4). The maximum number of holes/drops matches exactly with the predictions of linear theory. Further, the ripening and coalescence of drops leads to the formation of a single drop (profile 5). However, it will be shown below that, as film thickness approaches the critical thickness, the details of the initial surface heterogeneity critically affect the density of structures and their morphology. Figure 12 depicts the evolution of a localized surface heterogeneity in the liquid film that is a few nanometers thinner than critical thickness (ho < hc2). The localized depression quickly engenders a primary hole (profile 2) that continues to expand for a large distance (>λm, profile 3) before the formation of secondary satellite holes. Thus, the number density of holes decreases significantly compared to the predictions of the linear theory (profile 4, 6 holes on 10λm domain length) because of the difference between the various time scales of evolution involved in the formation of primary hole, the hole growth, the formation of satellite holes, and the spinodal time scale. The primary hole size can grow with the film thickness because of a greatly increased spinodal time scale (Figure 3b) on which the secondary holes make their appearance. For very thick films, close to a spinodal boundary such that ∂φ/∂h f 0 (in the DSSR), the number density of structures formed and their morphology shows even greater sensitivity to the amplitude and wavenumber of initial surface heterogeneity. Figure 13 shows the evolution of a film close to critical thickness (ho f hc2 ≈ 6.81 nm) with identical initial profile as in Figure 12. The evolution now results in a single hole in the 10λm domain (!) and any correlation with the linear theory is completely lost. Thus, films in the DSSR can give the appearance of nucleative dewetting despite the fact that these films are spinodally unstable. Figure 14 better illustrates the extreme structural sensitivity of the films nearing the DSSR when extremely small localized heterogeneities are present. Parts a-c of Figure 14 depict the

Figure 12. Evolution of a 5.5 nm thick film starting with a localized depression of 2% amplitude and a 1% random roughness on its surroundings, over a domain size of 10λm. The nondimensional times of evolution for profiles 1-5 are T ) 3.25 × 103, 5.43 × 103, 1.28 × 104, 1.84 × 104, and 3.89 × 108, respectively.

Figure 13. Evolution of a 6.5 nm thick film starting with 2% amplitude of localized depression and a 1% random roughness on its surroundings, over a domain size of 10λm. The nondimensional times of evolution for profiles 1-4 are T ) 9.9 × 103, 2.40 × 104, 4.90 × 104, and 2.64 × 105, respectively.

effects of the amplitude and diameter of a localized depression in a 5.5 nm thick film. The area of each simulation is identical, 25λm2. As seen in parts a and c of Figure 14, a single localized defect of negligible amplitude (0.1 nm) and diameter hm) close to the periphery of the spinodal region (∂φ/∂h f 0) starting with a localized depression of slightly larger amplitude as compared to the surrounding random heterogeneity. Parts

Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007 3115

Figure 14. Ordered structures in a 5.5 nm thick film: (a) Starting with a 0.1 nm amplitude depression of radius 0.7 µm and random roughness 0.05 nm on its surroundings. The nondimensional times of evolution for images 1-6 are T ) 1.70, 1.55 × 104, 1.70 × 104, 1.80 × 104, 1.87 × 104, and 1.96 × 104, respectively. (b) Starting with a 0.5 nm amplitude depression of radius 0.7 µm and random roughness of 0.05 nm on its surroundings. The nondimensional times of evolution for images 1-4 are T ) 9.22, 1.09 × 104, 1.63 × 104, and 1.81 × 104, respectively. (c) Starting with a 0.1 nm amplitude depression of radius 0.2 µm and random roughness of 0.05 nm on its surroundings. The nondimensional times of evolution for images 1-8 are T ) 1.62 × 102, 8.03 × 103, 1.42 × 104, 1.65 × 104, 1.76 × 104, 1.93 × 104, 2.03 × 104, and 2.27 × 104, respectively. The area of simulation domain in each case is 25λm2.

a-c of Figure 15 illustrate the effect of increasing the initial depth of the depression on the number density and the morphology while keeping the amplitude of the random perturbations unchanged. The local depression engenders a primary hole, which expands until the secondary holes surrounding the primary hole appear (profile 4 of parts a and b of Figure 15). The size of the primary hole increases with an increase of the amplitude of depression (profile 4), and for a 10% amplitude depression, only a single hole forms over 10λm length (profile 4 of Figure 15c). For thicker films, the small amplitude linear time scale for the instability is very high because of ∂φ/∂h f 0 and is, thus, very sensitive to the initial amplitude. Hence, the local small depressions quickly grow into holes and expand to occupy most of the space where the other

Figure 15. Evolution of a 6.5 nm thick film: (a) Starting with a 2% amplitude depression and random roughness of 1% amplitude on its surroundings, over a domain size of 10λm. The nondimensional times for profiles 1-4 are T ) 1.13 × 105, 1.32 × 105, 1.48 × 105, and 1.88 × 105, respectively. (b) Starting with a 5% amplitude depression and a 1% random roughness on its surroundings. The nondimensional times of evolution for profiles 1-4 are T ) 6.22 × 104, 7.72 × 104, 1.06 × 105, and 1.60 × 105, respectively. (c) Starting with a 10% amplitude depression and a 1% random roughness on its surroundings. The nondimensional times of evolution for profiles 1-4 are T ) 2.27 × 104, 3.68 × 104, 6.66 × 104, and 1.44 × 105, respectively.

potential holes of the linear theory could have eventually evolved. As the depth of local depression increases, this imbalance between the time scales for the formation of primary and secondary holes also increases and the hole-size distribution becomes bimodal and eventually acquires more broadband characteristics. Figures 16 and 17 depict some examples of the structural evolution when several localized interacting defects are present. In these examples, dimensional amplitude of local depression, amplitude of random surface heterogeneity, and dimensional separation distance between the two successive depressions are kept constant, but the film thickness is changed from deep inside spinodal territory to the defect-sensitive spinodal regime to analyze the sensitivity of structures to the initial conditions. For a lower film thickness (h ) 4 nm) in the deep inside spinodal regime, the local depressions initiate a cascade of holes that eventually resolve in a structure that is virtually indistinguishable from that produced without localized defects (profile 4 of Figure 17). A signature of local defects is visible only during the initial evolution, as reflected in the clustering of holes around the primary holes. However, increased film thickness under the same initial conditions engenders a distinct bimodal distribution of a smaller number of holes. Thus, the number density and size

3116

Ind. Eng. Chem. Res., Vol. 46, No. 10, 2007

Figure 16. Evolution of a 4.0 nm thick film, starting with a depression of 0.5 nm amplitude and a random roughness of 0.05 nm amplitude on its surroundings, over a domain size of 20λm (42 µm). The separation distance between two successive depressions is 10λm (21 µm). The nondimensional times for profiles 1-4 are T ) 0.0, 2.91 × 103, 5.76 × 103, and 7.76 × 103, respectively. Figure 18. Evolution of instability in films of different thicknesses when subjected to a single sine wave of amplitude 1% of initial film thickness, over a domain size of 10λm. Profile 1 is at the instant of rupture, and profile 2 is at the time when maximum number of holes formed at the interface. (a) Time evolution for 5.5 nm thick film, profiles 1 and 2 are at times T ) 33 551 and 304 367, respectively; the maximum number of holes formed is Nmax ) 8. (b) For 6.0 nm thick film, profiles 1 and 2 are at times T ) 77 007 and 112 574, respectively; the maximum number of holes formed is Nmax ) 4. (c) For 6.5 nm thick film, profiles 1 and 2 are at times T ) 404 776 and 646 116, respectively; the maximum number of holes formed is Nmax ) 2.

Figure 17. Evolution of a 5.5 nm thick film, starting with a depression of 0.5 nm amplitude and a random roughness of 0.05 nm amplitude on its surroundings, over a domain size of 20λm. The separation distance between two successive depressions is ∼21 µm. The nondimensional times for profiles 1-5 are T ) 6.24, 4.69 × 103, 8.42 × 103, 1.23 × 104, 1.56 × 104, and 9.10 × 108, respectively.

distribution of holes can be greatly affected by the distance between the local defects in the DSSR. 3.3. Periodic Surface Heterogeneity. Here we show that the results related to the extreme sensitivity of the DSSR to very small defects is not limited to random and localized defects but is equally shared even by initially periodic deformations of the film surface. This situation is of special interest in the context of mesopatterning of thin polymer films by controlled dewetting.9-11,19,33-36 Figure 18 summarizes the simulation results for the evolution of a single Fourier mode in different film thicknesses. As in the cases discussed previously, the number

density of holes is close to the linear theory prediction for relatively thinner films in the DIST but decreases with increasing mean thickness. This example illustrates that the regularity and length scale of an initially periodic structure are gradually lost in the DSSR, but not in the DIST, even by the introduction of extremely small numerical noise. This numerical experiment thus again points to the extreme sensitivity of the DSSR to small defects, which in this case are automatically introduced by the finite numerical representation of a periodic function. In summary, as one moves closer to the critical thickness (DSSR), the density of structures and their morphology show far greater sensitivity to the amplitude and wavenumber of the initial surface heterogeneity, unlike the structures in the deep inside the spinodal territory, where small defects do not produce significant deviations from the linear spinodal theory. 4. Conclusions On the basis of the large-domain nonlinear simulations, we have probed the structural response of spinodally unstable thin (