Unravelling the multilayer growth of the fullerene C60 in real time

Molecular semiconductors are increasingly used in devices, but understanding of elementary nanoscopic processes in molecular film growth is in its infancy. Here we use real-time in situ specular and diffuse X-ray scattering in combination with kinetic Monte Carlo simulations to study C60 nucleation and multilayer growth. We determine a self-consistent set of energy parameters describing both intra- and interlayer diffusion processes in C60 growth. This approach yields an effective Ehrlich–Schwoebel barrier of EES=110 meV, diffusion barrier of ED=540 meV and binding energy of EB=130 meV. Analysing the particle-resolved dynamics, we find that the lateral diffusion is similar to colloids, but characterized by an atom-like Schwoebel barrier. Our results contribute to a fundamental understanding of molecular growth processes in a system, which forms an important intermediate case between atoms and colloids.

is an extract from the (indeed very faint) colour image (main text, Fig. 2c). A clear peak splitting in the first monolayer is observed, which supports our analysis approach for the first monolayer. A clear peak splitting in the first monolayer is observed. (1)

Supplementary
Supplementary Note 1: Anti-Bragg intensity during the growth of the first monolayer of C 60 on mica We note the curious side-observation of the C 60 growth oscillations, where the reflectivity of the C 60 growth oscillations reaches almost zero after deposition of one C 60 ML (i.e. the reflected counts are reduced from to ). We emphasize that this has no particular implications for the growth itself, but we shall nevertheless briefly explain the origin of this scattering feature (see also reference [1]).
Using the formula for the Anti-Bragg oscillations given in the methods section, we fit the experimental growth oscillations and in particular the value of the minimal reflectivity at 1 ML coverage: Obviously, the scattering can vanish if the amplitudes involved are comparable and the phases are destructive (near ). If we insert the numbers appropriate for our system, we find indeed a nearly vanishing resulting intensity: Note that C 60 has a scattering amplitude that is a little smaller than the mica scattering amplitude ( ( ) ) and has a relative phase of (see Supplementary Fig. 1a for a graphical illustration of scattering amplitudes).
An alternative calculation of the absolute reflection intensities (with the same qualitative and quantitative result) uses the material densities of g/cm 3 for C 60 2 and g/cm 3 for mica 3 . From the density the mica and C 60 scattering length density (SLD) are calculated using with the bulk material density , the molar mass (taking into account the chemical composition of mica (KAl 3 Si 3 O 12 H 2 ) and the fullerene C 60 ) the number of electrons (for KAl 3 Si 3 O 12 H 2 and C 60 respectively) and the Thomson scattering length . This gives a mica scattering length density of in agreement with literature values 4 , and for C 60 a scattering length density of 1 is obtained. Using these values for the SLD and a C 60 monolayer thickness of (corresponding to the C 60 lattice constant) we arrive at SLD profiles shown in Supplementary Fig. 1b for the bare substrate (vacuum/ mica interface) and the one monolayer C 60 on mica structure (vacuum/ C 60 (1ML)/ mica). We use the Parratt formalism for the calculation of the reflectivity of stratified layers from the scattering density profile (see reference [5] for details of this recursive transfer matrix method).
Supplementary Fig. 1c shows the calculated x-ray reflectivity as a function of q for the bare substrate and 1ML on the substrate. For the C 60 on mica reflectivity a pronounced dip (Kiessig fringe) can be seen at the anti-Bragg point of C 60 ( which is equivalent to an incidence angle °) due to the C 60 layer reflectivity. From the simulated reflectivity we find a reduction by a factor of 31 at the anti-Bragg point ( for bare mica to a reflectivity of with one C 60 ML) neglecting the small contribution of surface roughness in the calculation. Comparison with the experimentally observed reduction by a factor of as obtained from the growth oscillations count rates for 0ML and 1ML in Supplementary   Fig. 1d shows good agreement between theory and experiment.

Supplementary Note 2: Epitaxial order of C 60 on mica
We have investigated the epitaxial order and thin film structure of C 60 on mica using Grazing Incidence X-ray

Supplementary Note 3: Equivalence of real-space and reciprocal-space methods
Real-space images of the first layers using AFM are unfortunately not possible due to strong dewetting effects in the first monolayers for C 60 on mica. Nevertheless, despite post-growth dewetting effects in the first layers, we could image the morphology of stable, thicker films with atomic force microscopy, as shown in Supplementary Fig. 4. In the kinematic approximation, the Fourier transform (FT) of the real-space structure corresponds to the diffusely scattered intensity in our x-ray experiments 6 . Comparing the diffuse x-ray scattering, which unfortunately for the thick films is close to the resolution limit and the FT of the corresponding AFM image for an 14 ML thick film ( °C, high deposition rate) we find very similar peak splitting of nm -1 , respective nm -1 . The line graph of the FT was convoluted with the resolution function of the experiment (normalized Gaussian function with a FWHM of nm -1 ). This good agreement between reciprocal-and real-space experiments confirms that our analysis determines correct lateral length scales.

Supplementary Note 4: Island shapes in experiment and KMC simulation
In general, the island shape can be obtained from the shape of the diffuse x-ray scattering (GISAXS) as the scattered intensity is composed of the structure factor (average island distance) and the form factor (island shape) of the thin film morphology. By assuming a certain island shape and (size) distribution of islands one can calculate the diffusely scattered intensity from this film morphology using dedicated programs like IsGISAXS 7 . Here we do not have to assume island shapes, but can directly use our KMC data to calculate the expected shape of the diffusely scattered intensity. In Supplementary Fig. 5  The comparison of the experimental lateral information (island density) and vertical information (anti-Bragg growth oscillations and layer coverages) was performed for rates of 0.1 ML min -1 and 1 ML min -1 and temperatures of 40 °C, 60 °C, and 80 °C. In addition to the data for a temperature of 60°C and rate of 0.1 ML min -1 shown in the main text, we present in Supplementary Fig. 6a and b an additional example for a temperature of 40°C and a rate of 0.1 ML min -1 .
In accordance with growth theories 8 and also correctly predicted by our KMC simulations, we find experimentally that the island densities increases (that is, the island sizes decreases) with lower substrate temperature ( Supplementary Fig. 6a). The anti-Bragg intensity ( Supplementary Fig. 6b) for °C shows distinct growth oscillations indicating layer-by-layer growth. From fitting the anti-Bragg growth oscillations using an analytical growth model 9 one can extract the layer coverage. The results support the simulated KMC layer coverages for °C as can be seen in Supplementary Fig. 6c.

Supplementary Note 6: KMC parameter choice: the correlation between energy barriers
To adjust the parameters , and appearing in the energy barrier ( ) (see equation (1) in the main text) we start from initial values suggested in the literature 10, 11 . We then optimize the parameter set to match, as accurately as possible, the experimental data for the island density and the filling fraction at °C and ML min -1 . This experimental data set has the best lateral and temporal resolution. In performing such an optimization, we have to note that the influences of the different energy barrier contributions on our observables are strongly correlated. For example, determines the free diffusion time and is therefore of prime importance for the calculated island density. The latter, however, is also influenced by : If the neighbour energy is sufficiently small, island nuclei can dissociate, which effectively reduces the island density. Therefore the influence of and on the island is correlated. Another example occurs during step-edge crossing: The energy barrier of this process is given by . However, the effective crossing rate is also influenced by the island morphology (which determines the probability of reaching a step), and thus, by . Due to this mutual influence it is clear that optimization of the energy parameters is rather challenging, one danger being that the resulting set may not be completely unique. We note, however, that the final parameter set selected in our study (see Fig. 1 in the main text) yields a satisfying match of the experimental data not only for the case °C and ML min -1 , but also at the other temperatures and adsorption rates considered. This "robustness" strongly supports our predictions. For the estimation of the error bars of the extracted energies their mutual correlation as well as the experimental confidence interval was taken into account. Notably, even a small alteration of one of the energy contributions (changes of the order of 20 meV) results in significant changes of the calculated morphology, i.e., in a deviation from experimental results. The confidence interval of the experimental island density data shown in Fig. 3 in the main text is calculated from the experimental uncertainties of x-ray wavelength, sample-to-detector distance as well as the fit uncertainty of the double peak distance in the diffuse x-ray scattering experiments. When we take the experimental confidence interval of the island density and anti-Bragg intensity into account, the error bar of the diffusion barrier increases to 40 meV.

Supplementary Note 7: Discussion on the diffusion barrier E D
Both Körner et al. 12 and Liu et al. 11 report a free diffusion energy eV; however,they also use an attempt frequency of Hz. Moreover, both studies are based on a hexagonal lattice under consideration of interstitial sites. Here, we neglect these sites, yielding a somewhat coarse-grained approach.
We note that without the coarse-grained approach it would not be possible to simulate such a large system for minutes to hours of experimental time. In one diffusion step on our coarse-grained lattice a particle overcomes the barrier eV twice. Moreover, one has to note that there are three diffusion options from the interstitial site. Since only one option leads to our coarse-grained destination site, an additional geometric factor of 1/3 needs to be included in the diffusion rate. Furthermore, taking the difference in the attempt frequency into account, we obtain the following estimate of a coarse-grained free diffusion barrier from the values reported in 11,12 This value lies within the error margins of our value . This estimate was obtained using K. (4)

Supplementary Note 8: Attractive interaction range of C 60 compared to atoms and colloids
The effective, centre-of-mass interaction between two C 60 molecules decays significantly faster to zero with the (centre-of-mass) distance r than that between atoms. Specifically, the potential has the form 13,14 (as obtained based on a calculation of Girifalco 15 ) , where is the depth of the potential and the centre-of-mass separation at which the potential is zero. The distance dependence of the attractive part of the potential, -1/r 9 , results from an angle-average over all the vander-Waals interactions (~ 1/r 6 ) between the individual carbon sites. (5), the attractive interaction between atoms can be described by a conventional Lennard-Jones potential which decreases as -1/r 6 for large distances:

Contrary to Supplementary Equation
To illustrate the different attraction ranges, the interaction potentials for atoms (here Argon) and C 60 are depicted in Supplementary Fig. 7. Note the significantly shorter attractive interaction range of C 60 compared to atoms.
An even shorter range of attractive interactions occurs in colloidal systems. Here, the attractive interaction between two colloidal particles typically originates from depletion forces induced by solvent particles. The range of the resulting attractive interaction is determined by the size of the solvent molecules, which can be orders of magnitude smaller than the size of the colloidal particles themselves 16 Supplementary Note 9: Calculation of particle-resolved dynamics for atomic systems: parameters except the ratio R allows us a direct comparison of single-particle dynamics despite the smaller time-and length scales of growth in atomic systems relative to C 60 . Clearly, this strategy is not suited to make predictive simulations of the growth of the atomic systems, but it does enable us to determine the influence of the difference in range of attractive interactions.

Supplementary Note 10: Origin of Ehrlich-Schwoebel barrier: Calculation of geometric contribution
The subsequent argumentation closely follows that given by Ganapathy et al. 18 for colloidal systems. We consider a C 60 particle moving on a surface formed by other C 60 particles. As a consequence of the short range of interactions the travelling particle tries to be in constant contact with two other particles. This effectively reduces local transport to a one-dimensional (1D) motion along a straight path, along which the potential landscape can be assumed to be constant. We note that the length of the 1D path between binding sites ( ) on an island is smaller than that of a path crossing the step-edge ( ), see Supplementary Fig. 8 for an illustration of such paths. Therefore, the travel time for a step-edge jump can be up to a factor of 〈 〉 〈 〉 ( ) longer, where 〈 〉 with r as the Arrhenius-type rate describing the surface processes 19 . Associated with this increase in diffusion time along the 1D potential is an increase of the probability to return to the original site 20 . As a consequence, the step-edge crossing probability effectively decreases. This consideration leads to an effective, geometrical Ehrlich-Schwoebel barrier determined by (( ) ) .
To quantify this geometry-induced energy contribution, we recall that our KMC simulations are based on a triangular lattice. This yields for the ratio a value of 1.4. This value is based on the considerations made in the supplemental material of Ganapathy et al. 18 , where they found for a hexagonal lattice that the pathlength of local motion in a 1D potential along the path crossing a step-edge (see b in Supplementary Fig. 8) is 2.8 times as long as the path-length involving diffusion between neighbouring sites (see a in Supplementary   Fig. 8). Since our simulation is restricted to a triangular lattice, an in-plane diffusion process involves two nearest-neighbour steps. This results in the ratio . Combining this result with Supplementary Eq. (7) we conclude that the geometrical Ehrlich-Schwoebel barrier is less than 50 meV for a temperature equal or less than 80°C.