Post-Deposition Wetting and Instabilities in Organic Thin Films by Supersonic Molecular Beam Deposition.

We discuss the formation and post-deposition instability of nanodrop-like structures in thin films of PDIF-CN2 (a perylene derivative) deposited via supersonic molecular beam deposition technique on highly hydrophobic substrates at room temperature. The role of the deposition rate on the characteristic lengths of the organic nanodrops has been investigated by a systematic analysis of atomic force microscope images of the thin films and through the use of the height-height correlation function. The nanodrops appear to be a metastable configuration for the freshly-deposited films. For this reason, post-deposition wetting effect has been examined with unprecedented accuracy throughout a year of experimental observations. The observed time scales, from few hours to months, are related to the growth rate, and characterize the thin films morphological reordering from three-dimensional nanodrops to a well-connected terraced film. While the interplay between adhesion and cohesion energies favors the formation of 3D-mounted structures during the growth, wetting phenomenon following the switching off of the molecular flux is found to be driven by an instability. A slow rate downhill process survives at the molecular flux shutdown and it is accompanied and maybe favored by the formation of a precursor layer composed of more lying molecules. These results are supported by simulations based on a non-linear stochastic model. The instability has been simulated, for both the growth and the post-growth evolution. To better reproduce the experimental data it is needed to introduce a surface equalizer term characterized by a relaxation time taking into account the presence of a local mechanism of molecular correlation.

The focus of the theoretical analysis reported in this supplemental material will be on the step barrier diffusion effect, known also as Ehrlich-Schwoebel (ER) barrier effect. This effect does not allow atoms to diffuse over the edge of a step on the surface creating an overall uphill current of diffusive particle flux. The step barrier diffusion effect creates mounds on the surface, introducing a new length scale λ, which measures the average mound separation. In this study, we will analyze the ER effect not only on the growth process when the particle flux impinges on the substrate, but also on the consequent relaxation dynamics after the flux switching off. In section I, the model and the time-dependent self-consistent approach are reviewed; in section II, the growth in the presence of the flux is theoretically analyzed; in section III, the relaxation dynamics after the flux switching off is discussed.

II. MODEL AND TIME-DEPENDENT SELF-CONSISTENT APPROACH
We will introduce a model for the modeling of the growth dynamics in supersonic molecular beam deposition of PDIFCN2. We will use a local stochastic continuum equation which includes the ER effect. The continuum equation involves the height profile h(x, t) and its derivatives (with respect to the time t and to the space represented by a bi-dimensional vector x) focusing on the local nature of diffusion and the step barrier diffusion effect. Actually, the process of diffusion is localized since one atom/molecule diffuses affecting only its neighbor. Indeed, the ER effect involves the diffusion of particles on the surface, therefore it gives rise to a characteristically local growth.
In this study, the growth process starts from t = 0, finishes at t c , switching off time of the flux, after which a relaxation dynamics takes place. The overall dynamics is modeled by the following equation 1,2 : where the first term describes the uphill growth due to the ER effect (∇ is the gradient operator), the second one is the Mullins diffusion term modeling the overall effect of surface diffusion. We remark that in the first term there is the sign −, therefore the ER effect introduces an instability in the growth process favoring the formation of mounds. The stochastic term present during the growth is modeled through the noise η(x, t), which is, as usual, assumed Gaussian: (2) where f (t) = θ(t−t c ), with θ(t) Heaviside theta function. In Eq. (1), the quantity F (t) models the particle flux: where R is the flux rate. The last term of Eq.
(1) provides a characteristic relaxation time τ : with G(t) indicating a primitive of the flux term This last term will be important after the switching off of the flux impinging on the substrate. The dynamics described by Eq. (1) is nonlinear due to the ER first term. Actually, the non-linear contribution in the denominator of the first term allows to provide a limit to the instability related to the mounding. We derive the temporal behavior of all the relevant dynamical quantities by using a time-dependent selfconsistent approach. The starting point is the observation that where m(t) is the surface slope. In order to solve the stochastic equation, we retain all the temporal fluctuations of the non-linear contribution in the first term of Eq.
(1), approximating its spatial fluctuations: Therefore, in the dynamical equation, one gets a term which depends itself on the equation solution. Therefore, we will solve the equation by using a self-consistent approach: the solution h(x, t) is accepted only when the slope function m(t) in the equation is, within a fixed tolerance, equal to that directly calculated through the height profile. In order to manage the inhomogeneous terms of Eq. (1), we divide the height profile into two contributions: where h (x, t) is the height profile affected by the stochastic noise, while G(t) = Rt c is the temporal contribution due to the external flux. We point out that G(t) = Rt c represents the average thickness of the organic film. The profile h (x, t) satisfies the following stochastic equation: with < |∇h| 2 >=< |∇h | 2 >= m 2 (t). In addition to the slope m(t), all the relevant quantities can be calculated through the height profile h : the interface width w(t) defined as the autocorrelation function R(r, t) as so that the height-height correlation function H(r, t) is (12) We solve Eq. (9) in the bi-dimensional momentum space q, therefore, we consider the spatial Fourier transformĥ(q, t) of h (x, t). One gets the following temporal equation forĥ(q, t): with q = |q| is the modulus of the vector q, andθ(q, t) is the spatial Fourier transform of the Gaussian noise term η(x, t), therefore The first-order temporal non-linear equation can be analytically solved determining the power spectral density function P (q, t), defined as We stress that the density function P (q, t) is proportional to D, the strength of the noise. Moreover, this density function allows to calculate all the relevant quantities of the dynamics. Actually, the function P (q, t) is the spatial Fourier transform of the autocorrelation function R(r, t) given in Eq. (11). Furthermore, it allows to directly calculate the interface width w(t) being Finally, we point out that the power spectral density function P (q, t) depends on the slope m(t) being derived from the solution of Eq. (13), but, at the same time, it allows to derive the slope: This turns out to be the equation that we have selfconsistently solved to determine the actual density function. Indeed, at the zero order, we have fixed m 0 (t) = 0, then, form Eq. (17), we get m 1 (t), and so on. A few steps of iteration are necessary to get a function m i (t) very close to m i+1 (t) (a tolerance is given in the norm of the difference between m i (t) and m i+1 (t)). It is the behavior at long times which is progressively improved with increasing the iterative steps. We stress that the self-consistent approach allows to carefully analyze the solutions for very long times. Most organic films analyzed in the experiments have an average thickness T around 20 nm. Since the flux rate R used in the experiments is known, one can estimate the switching off time t c = T /R. The rates R used in the experiments vary from around 0.04 nm/min to 0.22 nm/min. As a consequence, the times t c derived from experiments go from around 90 min to 500 min.
With the parameters ν and k present in Eq. (13), one can obtain the length l 0 = 1/q 0 = 2k/ν and the time t 0 = l 3 0 / √ νk = 2 √ 2k/ν 2 . The quantities l 0 and t 0 provide a term D 0 = l 4 0 /t 0 = √ 2k, which has the same dimensions as the term D, which measures the strength of noise. For the dynamics analyzed in this paper, we set the following values of l 0 and t 0 as unit length and time, respectively: l 0 = 77 nm and t 0 = 1.38 days. Then, the free parameters of the theory are R (or t c ), D, given in terms of D 0 , τ , expressed in terms of t 0 . For the dynamics analyzed in this paper, D is of the order of 0.001 − 0.01 D 0 , while the relaxation time τ is assumed to be quite large being of the order of 6 − 7 t 0 .

III. GROWTH IN THE PRESENCE OF THE FLUX
In this section, we analyze the growth process in the presence of a flux impinging on the substrate, therefore in the case of times t < t c . We will analyze the growth features for t c larger than 500 min. Since the values of the relaxation τ are assumed large, the effects of this parameter on the growth dynamics are completely negligible for times t < t c .
We will investigate the interface width w(t) focusing in particular on the regime of power-law behavior and the related exponent β in the presence of the flux: Another important quantity is the lateral correlation length ξ(t) defined through the autocorrelation function of Eq. (11) such that R(ξ(t), t) = 1/e 0.36787.
We will study the correlation length ξ(t) focusing in particular on the exponent 1/z which characterizes the regime of power-law behavior in the presence of the flux: In the case of mounding, another important length is the average mound separation λ(t), which plays the role of the characteristic wave-length. At a fixed time t, given the width w(t), the length λ(t) is determined assuming, in the limit of large distance (|r| >> ξ(t)), the following behavior of the height-height correlation function of Eq. (12): where J 0 (x) is the zero-order Bessel function. We will analyze λ(t) focusing in particular on the exponent p which characterizes the regime of power-law behavior in the presence of the flux: When the mound instability is absent, the self-affine behavior takes place and the dynamic scaling is present, that is p = 1/z. In general, dynamic scaling does not hold for mounded surfaces. In the following, we will find slight deviations from the situation of dynamic scaling. In Fig. 1, we show the average mound separation λ, the lateral correlation length ξ, and the interface width w as a function of time. At first, we notice the different orders of magnitude of these lengths. Actually, the width w is of the order of a few nm, the correlation length ξ is of the order of ten nm, the wave-length λ of the order of one hundred nm. We stress that, in case of mounding instability, the correlation length ξ provides the average size of a mound. Therefore, one expects that ξ is smaller than λ.
In the presence of a flux, after a short transient, all the three lengths follow a power-law t γ . We have fitted all the three time functions in the experimentally relevant range of t c , therefore, as shown in Fig. 1, we consider a time interval between 1 hour and 6.5 hours. The exponent β of the interface width is about 0.29, a value perfectly compatible with that found in the literature neglecting the correction m 2 (t) in the ER term of the dynamical equation 3 . This result is expected since, up to times of the order of a few hours, the values of m 2 (t) are always smaller than 1. Therefore, the self-consistent procedure provides weak perturbations when the flux is present. The other two lengths increase more quickly as a function of times. Actually, the exponent p of λ is about 0.433, while the exponent 1/z of the correlation length ξ is around 0.405. Therefore, the dynamics scaling is only weakly violated. In the main text, we plot the correlation length ξ and the mound separation λ as a function of the flux rate making a comparison with experimental data.

IV. RELAXATION AFTER THE FLUX SWITCHING OFF
In this section, we analyze the relaxation dynamics from the flux switching off up to times larger than 20 days. Our self-consistent procedure is able to cover this very large time range providing accurate numerical results. We stress that the slope correction m 2 (t) in the ER term of the dynamical equation strongly affects the relaxation. First, we will consider the case of an infinite τ (1/τ 0), then the effects of a finite τ . We will show that a large but finite τ plays an important role in order to interpret the dynamical relaxation.
In Fig. 2, we show the lateral correlation length ξ as a function of time for different flux rates in the case 1/τ 0. We theoretically analyze the cases of intermediate (R = 0.1 nm/min), very large (R = 1 nm/min) and very small (R = 0.01 nm/min) flux rate. Focus is on the relaxation dynamics, therefore we plot the behavior of ξ only for t ≥ t c . We find that the asymptotic behavior of ξ at very large times is independent of the rate. Moreover, at very large times, ξ increases slowly reaching a nearly constant asymptotic value close to 200 nm. Therefore, as a result of the relaxation, the average size of a mound can be a stable quantity as a function of time. In Fig. 3, we show the average mound separation λ as a function of time for different flux rates in the case 1/τ 0. As in the case of the correlation length, we theoretically analyze the cases of intermediate (R = 0.1 nm/min), very large (R = 1 nm/min) and very small (R = 0.01 nm/min) flux rate for t ≥ t c . In analogy with ξ, the asymptotic behavior of λ is poorly independent of the rate at very large times. In contrast with ξ, at very large times, λ increases as a function of time even if the increase is quite slow. Therefore, during the relaxation, the average mound separation is a slow function of time.
One can expect such a behavior since the mound instability is always active in the system even in the absence of the flux.
We point out that the values of λ for times around 20 days are of the order of 650 nm, a value that is not much larger than 2π/q 0 = 2πl 0 = 2π 2k/ν 480 nm. This last value represents the natural average mound separation in a linear equation model without the corrective term m 2 (t) in the ER term 3 . We can interpret the numerical results considering that, due the slope corrective term, there is an effective ν ef f ν/(1 + m 2 ef f ), where m ef f is an effective value determined by the type and duration of the dynamics. Therefore, one expects a smaller ν ef f , therefore a larger λ ef f ∝ 1/ν ef f .
The increase of the values of λ as a function of time is obtained by analyzing the behavior of the height-height correlation function in the limit of large spatial separation. In Fig. 4, we show the height-height correlation It is apparent that the oscillatory behavior for large distances is also more marked with increasing time. Moreover, we notice that, just after the flux switching off time, the values of the correlation function get reduced. Only after a transient time of the order of a few days, the values of H starts increasing. We emphasize that this behavior is in good agreement with experimental data. This trend can be almost entirely ascribed to the time behavior the width w(t) (shown in the main text) being the correlation function H proportional to w 2 (t). As shown in the main text, the effect of a finite τ reduces w(t) at large times, therefore it will further decrease the values of the correlation function H at intermediate flux rate, getting a better agreement with experimental data.