Direct imaging of structural disordering and heterogeneous dynamics of fullerene molecular liquid

Structural rearrangements govern the various properties of disordered systems and visualization of these dynamical processes can provide critical information on structural deformation and phase transformation of the systems. However, direct imaging of individual atoms or molecules in a disordered state is quite challenging. Here, we prepare a model molecular system of C70 molecules on graphene and directly visualize the structural and dynamical evolution using aberration-corrected transmission electron microscopy. E-beam irradiation stimulates dynamics of fullerene molecules, which results in the first-order like structural transformation from the molecular crystal to molecular liquid. The real-time tracking of individual molecules using an automatic molecular identification process elucidates the relaxation behavior of a stretched exponential functional form. Moreover, the directly observed heterogeneous dynamics bear similarity to the dynamical heterogeneity in supercooled liquids near the glass transition. Fullerenes on graphene can serve as a new model system, which allows investigation of molecular dynamics in disordered phases.

D isordered states of materials, such as supercooled liquid or glass, display peculiar non-equilibrium behavior and have been intensively investigated for decades [1][2][3][4] . Various glassy materials are widely used in industrial applications and consumer products, and a fundamental understanding of the disordered non-equilibrium structural phase may also facilitate advancement in material processing and fabrication. Structural rearrangements, which are often spatially heterogeneous, govern the various properties of glassy systems 2,4-6 . Yet, direct visualization of these dynamical processes has been difficult, and diffraction-based analysis such as structure factors or intermediate scattering functions has been mainly utilized only to yield the spatially averaged signals 1,2 . To complement this issue, researchers have utilized computer simulations 6,7 or other macroscopic model systems such as granular 8 or colloidal particles [9][10][11] , providing verification on spatially heterogeneous dynamics and increasing characteristic length scales near the glass transition.
The main experimental limitation in directly imaging atoms and molecules in disordered states can be overcome by modifying the sample geometry, to an atomically thin two-dimensional (2D) glassy system. Until now, 2D silica glass 12,13 , 2D carbon glass 14 , and Si atoms at the surface of amorphous Si 15 have been successfully visualized using transmission electron microscopy (TEM) or scanning tunneling microscopy (STM). Although direct imaging of atomic structure has been successfully performed in these atomic disordered systems, in-depth observation and analysis of both structure and dynamics at the atomic or molecular resolution are still mainly lacking. In particular, the strong covalent bond in the previously-studied atomic glass leads to relatively slow structural evolution and dynamics under experimental conditions, limiting systematic studies on the dynamical behavior.
Here, we prepare a C 70 molecular system on graphene, and directly visualize both the structural and dynamical evolution of the system at molecular resolution. The relatively weak van der Waals interaction between C 70 molecules can be perturbed using an electron-beam (e-beam) during aberration-corrected transmission electron microscopy (acTEM), emulating the melting process of the molecular crystal. Our computerized method precisely identifies molecular positions in the disordered state, and the pair correlation functions of molecules clearly show the short-range liquid-like ordering. Time-dependent relaxation behaviors of the molecular structure are studied in-depth by van Hove correlation functions, which clearly shows the relationship between the local structure ordering and the dynamical behavior of the system. Real-time tracking of individual molecules also allows us to extract the spatially heterogeneous dynamics during the melting process. Our study demonstrates that fullerenes on graphene can serve as a new model system for investigation of super-cooled liquid and glass at molecular resolution.

Results
Preparation of C 70 crystals on graphene. Fullerenes have relatively high robustness to e-beam under TEM imaging conditions compared to other organic molecules and can serve as model molecules [16][17][18][19] . During e-beam irradiation, energetic electrons can occasionally transfer significant momentum and energy to the fullerene molecules via elastic collisions, generating structural displacements and dynamical motions 20 . In particular, the weak van der Waals interaction between fullerene molecules can be easily perturbed with an e-beam, inducing strong dynamical behavior. Indeed, researchers have previously observed an e-beam induced molecular dynamics of C 60 in one-dimensional confinement of carbon nanotube 21 . In our study, we chose C 70 deposited on graphene as an ultra-thin model molecular system. Graphene has high electrical/thermal conductivity and can further reduce the structural damage to samples because of fast energy transfers [22][23][24] . The atomically flat surface of graphene facilitates the dynamical behavior of the molecular liquid under imaging conditions, allowing systematic study of its molecular structure and dynamics. Compared to C 60 fullerene 16,18,25 , C 70 has a slight anisotropy in its molecular shape, which may also contribute to the more pronounced molecular movement.
We deposited C 70 molecules onto a graphene membrane by thermal evaporation to make a thin film (see methods section for detailed information) 18 . Using electron diffraction, we find that the nearest neighbour distance between molecules in the initial C 70 film was 1.07 nm without any sign of anisotropy in terms of C 70 orientation ( Supplementary Fig. 1). The observed thin-film (10 nm thick) C 70 crystal structure is consistent with the hightemperature crystal phase of a face centered cubic (fcc) 16 . The film also exhibits a uniform and well-ordered structure ( Supplementary Fig. 1). For the main study on the molecular structure and dynamics of liquid-like state, C 70 samples with 3~5 nm thickness (3~5 molecular layers) were used.
Melting of C 70 crystals by e-beam irradiation. A well-ordered C 70 packing structure can be perturbed with e-beam irradiation (Fig. 1a). Figure 1b-d show TEM images of C 70 film on graphene, revealing the structural transition from ordered to disordered configurations. E-beam irradiation drives the molecular crystal to a more disordered state, emulating the melting process of the molecular crystals. We note that the process is driven mainly through decreasing the molecular ordering by the random energy transfer of e-beam to molecules, not through actual heating of the system 13,14 . This raises the possibility that the effects of e-beam irradiation may not be the same as thermal effects. For example, in the case of single-layer graphene, it was shown that the population of e-beam induced defects deviates from a thermally induced Boltzmann distribution 26 . Nonetheless, e-beam irradiation can drive the system into high-energy states, similar to what can be expected in a thermally activated system. Such similarities have been studied theoretically in 2D silica 27 , and experimentally observed during the transformation of various graphitic nanostructures [28][29][30] .
As the C 70 molecular long-range order diminishes under prolonged e-beam irradiation, we observed a change in the local C 70 film thickness. The bare graphene surface, which can be regarded as a pore in the C 70 film, was also observed locally as shown in Fig. 1d. The non-uniform film thickness under e-beam irradiation mainly results from the molecular rearrangements on the graphene surface rather than the molecular ejection because the energy barrier for the molecular migration is much lower than that for the molecular desorption from the surface. 20,31 The molecular rearrangements seem to be driven by the stronger C 70 -C 70 interactions than the C 70 -graphene interactions 31 . Around the pores in the C 70 film, the local thickness of the film was sub-monolayer, and isolated C 70 molecules with circular shapes could be clearly observed. The observed molecular shape (Fig. 1g) and the line intensity profile along the molecule (Fig. 1h) were consistent with the TEM simulation data (Fig. 1e, f). Figure 2 is a series of TEM images showing the behavior of C 70 molecular crystal at the very early stage of melting, namely stage 0. TEM images of C 70 initially show highly-ordered molecular positions together with clearly visible circular molecular shape. Later, molecular movements were induced by e-beam irradiation, which results in local heterogeneous disordering of the molecules. The area bounded by red lines in Fig. 2 indicates the disordered molecular regions. We find that the disordered regions nucleated locally and were growing with fluctuations in their shape and eventually merging with adjacent regions over time as shown in Fig. 2 and Supplementary Movie 1. The nucleation and growth of a disordered phase out of a homogeneous crystalline phase and the observed microstructural phase coexistence with clear phase boundaries strongly suggest that the observed phase transition has the characteristics of a thermodynamic first-order phase transition 32,33 . This observation supports that the e-beam induced crystal-to-liquid melting in our study has similarity to the conventional crystal-to-liquid melting induced by thermal activation.
Identification of C 70 positions and pair correlation. The welldefined nearly circular shape of the C 70 molecules, as shown by the dark circular line at our imaged condition (defocus value at −13 nm), can be used to efficiently identify molecular positions, even in an area where the direct recognition of molecules is difficult due to multiple molecular overlaps. To identify molecular positions reliably, we devised an image processing scheme whereby circles with a predesignated radius r were automatically identified. To achieve this, we processed TEM images through a two-phase Mumford-Shah (MS 2 ) model 34 and calculated the probability density function at each image pixel, PDF(x, y), which is a parameter indicating the probability of molecular presence at a certain pixel location (x, y). (See Supplementary Note 1 and Supplementary Fig. 2 for details). With a proper choice of PDF threshold, we could assign molecular positions in twodimensional image space.
To validate our method, we first applied our image processing to simulated TEM images using a model with disordered molecules on graphene, as shown in Fig. 3a, b. The molecular model with~3-nm thick fullerene film was constructed by Monte Carlo simulations, in which a reasonable three-dimensional (3D) molecular pair correlation was used as shown in Fig. 3c correlation function (PCF) is given by where r ij is the distance from the ith to the jth molecules, N is the total number of molecules, and ρ is the number density of molecules. The identification process was verified using simulation images with different noise levels ( Supplementary Fig. 3). We found that the process identifies the centers of molecules with precision higher than 92% for a model with three molecular layers and the main inaccuracy of the identification process is originated from uncounted molecules due to significant molecular overlaps ( Supplementary Fig. 3). The main effect of the uncounted molecules can be seen from the undercounted data points at r < 0.3 nm in 2D-projected PCF (2D g 2 (r)) as shown in the inset of Fig. 3d. Nevertheless, the general features in 2D g 2 (r), including the peak intensity and position, were well-captured by our identification process. The identification process was also validated with control images of amorphous carbon and computer-generated random noise (Supplementary Figs. 4 and 5). The well-known Hough transform 35,36 and another scheme (simple circumference transform) were also tested, but we find that MS 2 model yielded the most reliable identification of C 70 molecules.
With the application of the MS 2 model and PDF calculations, experimentally obtained TEM images were processed (Supplementary Movies 2-4) and the centers of molecules were identified, as shown in Fig. 3e, f. Figure 3g presents PCF at different stages of imaging (stage from 1 to 3), where the structural ordering was diminished as the molecular crystals underwent melting under ebeam irradiation. The PCF at its fully disordered state exhibited peaks around r = 0.5 and 1.0 nm and converged to unity in the range r > 2.0 nm. This observation clearly shows that the C 70 molecular system hosts a liquid-like short-range ordering under e-beam irradiation. It is noticed that the peak positions in 2D PCF were slightly different from the expected positions (red lines) calculated from projection image using close-packed C 70 fcc molecular assembly. In particular, the peak positions were slightly down-shifted, showing an apparent lattice contraction compared to the original C 70 crystal due to 2D projection of 3D counterpart 37 . We note that the obtained PCF in Fig. 3g was calculated from 2D-projected TEM images, which shows some differences from the 3D PCF as shown in Fig. 3c, d. Tomography-based TEM imaging 38 is potentially available to obtain a static molecular 3D PCF, but it is currently challenging to apply the technique to our dynamic molecular systems due to possible e-beam induced motions during extended experimental imaging time.

Time-dependent correlation function of molecular structures.
The detailed dynamical behaviors of C 70 molecules can be studied using time-dependent correlation functions. Van Hove correlation functions (vHCF) were obtained using Figure 4a shows g(r,t) = G(r,t)/ρ at the state where the disordered state was developed from the C 70 crystal (stage 2). We note that the proper calculation of vHCF involves the overall drift correction of time-series TEM images, which is described in detail in Supplementary Fig. 6. The peaks at r~0.5 nm and r~1.0 nm show rapid decaying behaviors and converge to unity (Fig. 4b).
The dynamical structural relaxation characteristics can be further investigated by analyzing the shape of the relaxation function. Figure 4c shows the decaying vHCF plots of the first peak (at r~0.5 nm) at the different stages of melting. The relaxation behavior can be fitted using g t fitting parameters τ and β. The fitting functional form does a good job of capturing the experimental relaxation, resulting in approximately τ~10 s and β~0.6. The observed non-exponential decaying function with β < 1 (stretched exponential decaying) is a strong indicator of the heterogeneous dynamical behavior of the system [3][4][5]39 . The observed dynamical relaxation characteristics exhibited correlations with local structure ordering. As shown in Fig. 4d, the relaxation time τ becomes smaller as the local structural order parameter, PCF 1 st peak, diminishes. In situations where the structural order parameter is high, the relaxation dynamics can be suppressed due to the higher energy barrier originated from the larger molecular coordination number, which can be obtained from R r¼1:1nm On the other hand, as the molecular structure undergoes disordering, the energy barrier for molecular movement can be lowered due to the smaller molecular coordination number. This result implies that the relaxation process is highly sensitive to local structural ordering and the observed heterogeneous dynamical behaviors can be originated from the heterogeneous structural ordering.
Spatially resolved heterogeneous dynamics in C 70 liquid. The heterogeneous dynamical molecular behaviors can be further directly accessed by visualization of molecular movement during melting. With the ability to track dynamics at molecular resolution, we visualized the 2D diffusional behaviors of C 70 molecules. Figure 5a shows 2D molecular trajectories of C 70 molecules at a relatively early stage of melting (stage 1), demonstrating that the structural and dynamical behaviors of molecular disordering display spatially heterogeneous process. The central area of Fig. 5a shows a more ordered structure together with suppressed dynamical behavior, whereas the non-central area shows a more disordered structure with enhanced molecular dynamics. The positions of molecules with pronounced dynamics also display this heterogeneity, as shown in Supplementary Fig. 7. The zoomed-in molecular trajectories (Fig. 5b) show that the molecular diffusion is associated mainly with wiggling movements and sporadic large jumps. This non-trivial behavior is reminiscent of the cage rearrangement observed in model colloidal systems near the glass transition 10,11 .
Detailed diffusional behaviors of molecules can be examined using the mean-squared displacement of molecules defined by Fig. 5c. The meansquared displacement at two different stages during the melting process clearly shows that the molecular diffusion became more pronounced at the later stage of observation. This is consistent with our analysis and conclusion from vHCF in Fig. 4. Interestingly, in Fig. 5c, the slope of the plots (red dashed line) begins with a lower value and then approaches the slope of the blue lines over time. The blue dashed lines indicate the diffusional behavior obtained from free diffusional form 6 . The time scale of the transition from sub-free to free diffusional behavior is~10 s, which is consistent with the relaxation time scale observed in Fig. 4. Heterogeneous dynamical behavior can also be studied using the non-Gaussian parameter, α 2 ¼ r 4 t ð Þ h i= 3r 2 t ð Þ h iÀ 1; which quantifies deviations from a Gaussian distribution 5,10,11 . Calculated α 2 at the stage 1 shows a peak value of around 2 at t~10 s as shown in Fig. 5d. The similar α 2 feature showing a peak at the characteristic relaxation time scale was previously observed for the alpha-relaxation of various glass systems 5  The vertical error bars represent the ranges of τ and β obtained from curve fitting to the data points in Fig. 4c with a stretched exponential function heterogeneous dynamics, and our molecular system bears some similarity to the dynamics observed in supercooled liquids near the glass transition 3,6 . We note that the effect of overall drift was eliminated by the drift correction during the post-imaging process as shown in Supplementary Fig. 6.

Discussion
In conclusion, we observed an e-beam induced first-order transition-like crystal melting of C 70 molecular system and analysed its molecular dynamics with single-molecule sensitivity. The pair correlations for C 70 molecules were calculated using an automated molecular position identification process. The spatially heterogeneous dynamical behavior of these molecules bore the similarity to the dynamical heterogeneity observed in supercooled liquids near the glass transition. Considering the possibility of using various molecules with different anisotropy (C 60 , C 72 , and C 82 ) and modifying interactions between molecules through surface/internal functionalization (C 60 F 48 and M@C 82 ) 16 , fullerenes on graphene can serve as a new model system for investigation of molecular glass and supercooled liquid, providing unprecedented real-space imaging of dynamical heterogeneity. Future work combining the e-beam irradiation with in situ heating will make it possible to reveal the details on the thermodynamic equilibrium and non-equilibrium properties of fullerene model systems.

Methods
Sample preparation. Graphene was synthesized using chemical vapor deposition (CVD) 40 . Twenty-five micrometer-thick copper foil was used as the synthesis substrate. CVD graphene was transferred to Quantifoil holey carbon grids via direct transfer 18 . C 70 films were deposited onto graphene TEM grids by thermal evaporation. Before the thermal evaporation process, TEM grids were pre-annealed in the air at 200°C for 30 min to minimize surface adsorbate on graphene. The C 70 films with thicknesses ranging from 0.5 to 10 nm were deposited at the deposition rate of 0.05 Å s −1 under a vacuum pressure of 2 × 10 −6 Torr. The graphene substrate was held at 110°C during the deposition.
TEM imaging and simulation. TEM imaging was performed using a FEI Titan equipped with an image aberration corrector operated at 80 kV and a JEOL ARM 200 F equipped with image and probe aberration correctors operated at 80 kV. Selected area electron diffraction was performed with a FEI Tecnai. The electron dose during TEM imaging was~2 × 10 4 e nm −2 s −1 . TEM videos were recorded using the CamStudio program with a temporal resolution of around 0.2 s per frame. TEM image simulations were performed using MacTempas software with experimental imaging conditions. The simulation images were obtained at a defocus value of −13 nm.
TEM image analysis. Molecular position identification was performed using the two-phase Mumford-Shah (MS 2 ) model and the calculation of the probability density function (PDF) for molecular presence at each pixel. First, the MS 2 model calculates the support function of an experimental TEM image. Based on the dark circular shape of fullerenes on the image, the annulus support with a given radius was applied to each image pixel. Second, the PDF was computed using the calculated support function and the chosen reference molecule, which described the probable molecular positions. To find positions that otherwise could be missed due to the aggregation of molecules, bright artificial molecules were inserted to the positions obtained from PDF, after which the PDF was re-calculated. Doublecounted positions were determined if the PDF value was more than 98% of the PDF maximum. Drift correction of time-series TEM images was performed using a custom ImageJ macro. To quantify overall drift for each frame, tracking of molecular positions was executed for all frames, and the average displacement of the tracked molecules at a given frame was calculated along both the x and y directions. The overall drift was compensated using the calculated average displacement for a given frame. These processes are described in greater detail in the Supplementary Information section.

Data availability
The authors declare that the data supporting the findings of this study are available within the Supplementary Information files and from the corresponding authors upon reasonable request.

Code availability
The codes used for data are available from the corresponding authors on reasonable request.