Field emission microscope for a single fullerene molecule

Applying strong direct current (DC) electric fields on the apex of a sharp metallic tip, electrons can be radially emitted from the apex to vacuum. Subsequently, they magnify the nanoscopic information on the apex, which serves as a field emission microscope (FEM). When depositing molecules on such a tip, peculiar electron emission patterns such as clover leaves appear. These phenomena were first observed seventy years ago. However, the source of these emission patterns has not yet been identified owing to the limited experimental information about molecular configurations on a tip. Here, we used fullerene molecules and characterized the molecule-covered tip by an FEM. In addition to the experiments, simulations were performed to obtain optimized molecular configurations on a tip. Both results indicate that the molecules, the source of the peculiar emission patterns, appear on a molecule layer formed on the tip under strong DC electric fields. Furthermore, the simulations revealed that these molecules are mostly isolated single molecules forming single-molecule-terminated protrusions. Upon the excellent agreements in both results, we concluded that each emission pattern originates from a single molecule. Our work should pave the way to revive old-fashioned electron microscopy as a powerful tool for investigating a single molecule.

the top of the protrusions. However, it was not clear whether the molecules at the top are single or clustered [7][8][9][10][11][12][13] , unless one could undertake an exceptional study with an ultra-sharp tip 14 or a single molecule emitter 15 . For instance, the authors in reference 10 consider the protrusions to consist of a cluster of oxygen atoms, and each spot in the emission patterns representing a single atom. This interpretation, however, has not explained all the characteristics of the observed emission patterns.
Unlike an STM, an FEM shows only molecules with high electric fields at the top of protrusions. This limitation in the experimental information makes it difficult to reveal the detailed structures of the protrusions. In this study, together with our FEM experiments, we performed simulations to determine the optimized molecule configurations on a tip under strong DC electric fields, and thus, we aimed to learn about the sources of the emitted electrons. Here, we employed C 60 molecules because of their highly isotropic properties. Based on its isotropic aspect, we assumed the molecule as a simple sphere in the simulations 16 and could thus vastly reduce the computational cost. The experimental results indicated formations of protrusions on a buffer layer, as discussed in the previous work 10 . The simulations further revealed that each protrusion is primarily a single molecule, and the protrusion is the potential source of electrons.

Results
Experimental: summary of characteristic FEM images and our interpretations. We shall first examine the experimental results. After depositing molecules on a tip, different FEM images could be observed depending on the amount of deposited molecules and DC electric fields. These images are summarized in Fig. 2, together with the schematic diagrams of our interpretations. Deposition of molecules is done by resistively heating the evaporator boat, as shown in Fig. 1. In the experiments presented in this article, the temperature of the evaporator boat during the deposition was 420 °C and we gauged the amount of deposited molecules by a deposition time. A clean tungsten carbide sample was mainly used in our experiments (see the Method section for details on the experiments). The emission pattern from the tungsten carbide is shown in image I at the leftmost panel in Fig. 2. After depositing molecules for 2 min, the lower side of the tip became darker, as shown in image II. This is because of the formation of a molecular layer around the lower half of the apex, as schematically depicted below the FEM image. (We call this layer a buffer layer, which will be explained later.) Since the molecule source was installed below the tip, molecules were deposited mainly on the lower side of the tip. Note that we applied -500 V at the tip during the deposition. Under this condition, molecules are attracted towards the tip apex, and some molecules can be deposited even on the upper side 17 .
After depositing molecules for 4 min, peculiar emission patterns could be observed over the molecule layer, as shown in image III. It should be emphasized that the electron emission from the tungsten carbide surface could be slightly discerned from the upper half of the apex, and their signals got weaker in the lower half, indicating that there was a molecule layer mainly at the lower part of the apex. It should be also mentioned that the molecule layer was more firmly bonded to the surface than the molecules above the first layer, as we will discuss later. The molecules above the first layer are considered to be chemically pristine fullerenes (though, when the applied field is present, they may be partially charged and/or have a dipole moment) because previous studies have revealed that the electronic properties of molecules above the first layer are the same as that of the pristine C 60 , which was true for various metal substrates, including tungsten [21][22][23] . Hence, we call the first layer a buffer layer to distinguish it from the other layers. As depicted in the interpretation schematic for image III, those molecules on the buffer layer were well separated from each other, and thus they formed protrusions of single molecules. The peculiar emission patterns originated from the protrusions.
Under the same deposition time but with lower DC electric fields, the FEM image became different, as shown in image IV in Fig. 2. In this image, a flock of random and tiny spots on the lower half of the image can be observed. We believe that several molecule layers were formed on the tip, and electrons were emitted from some protrusions that appeared on the topmost layer, as depicted in the figure. Those molecules above the first layer were weakly bound to the surface, and they were evaporated under strong DC electric fields. This transition  Experimental: evaporation of molecules by DC electric fields. Figure 3a is an FEM image of clean tungsten carbide. For illustrative purposes, the emission areas of the tungsten carbide have been highlighted by white lines, and the same lines are indicated in Fig. 3b-e. After depositing molecules for 4 min, we saw a flock of tiny spots under a weak DC electric field, as shown in Fig. 3b. These spots are surrounded by dashed lines. With increasing electric fields, the flock of spots gradually disappeared from the top of the tip apex, and the dashedline area gradually moved downwards, as shown in Fig. 3c-e. As the flock of spots disappeared, the signals from the substrate became visible. The applied electric fields and the corresponding timestamp of related experiments are shown in Fig. 3f. At the highest fields in Fig. 3e, we could observe the peculiar emission patterns such as a large round spot and a two-leaf pattern surrounded by dashed circles. Additionally, the electron signals become weaker below the dashed line, which indicates that a buffer layer remains mainly below the line. This is the same situation as image III in Fig. 2. It should be mentioned that the dashed line in Fig. 3e does not indicate the clear boundary of the buffer layer but indicates that the density of molecules in the buffer layer becomes sparse above the line. The buffer layer molecules should exist even on upper side of the tip as explained in image II of Fig. 2.
We believe a molecule that is the source of the two-leaf pattern in Fig. 3e is situated on a buffer layer, though the pattern appeared above the dashed line. The reason for this is that we have rarely seen those peculiar patterns at the low coverages where the formation of the buffer layer cannot be observed. These results in Fig. 3 indicate that the flocks of spots were signals from the molecules on the buffer layers, and they were weakly bound to the surface. These weakly-bound molecules would be forming more than a monolayer on the buffer layer as the emission patterns were rather different from the ones on the buffer layer. Note that these molecules on the buffer layer consistently evaporated from the top of the apex, because the DC electric field was the strongest at the top of the tip apex and got weaker toward the shank side.
Experimental: evaporation of molecules by heating. The buffer layer could not be removed by the DC electric fields though we tried up to around 4.5 V/nm as in image II of Fig. 2. However, it could be removed The temperature of the evaporator boat was 420 °C. The value of the DC electric field at the top of the tip apex, F DC is denoted in each image. A schematic diagram of the interpretation of each image is also depicted. F DC was determined using the following equation: F DC = V tip kr , where V tip is the voltage applied to the tip, r is the tip radius, and k is a parameter that varies depending on geometrical factors, such as tip shape, geometries of surrounding electrodes and polar angles of the sphere of the tip apex 7 . In this particular case, we determined k top , which is defined at the very top of the tip apex. The tip radius was separately measured once using a scanning electron microscope (SEM). In addition, before each experiment, we recorded threshold tip voltages at which the electron signals from the field emission could be seen on a phosphor screen. Here, we assumed that the F DC for the threshold voltages is constant as long as the tip surface conditions, such as material and cleanliness, are the same. The threshold voltage is inversely proportional to the tip radius under this assumption. At each experiment, we determined the tip radius by multiplying the tip radius that was measured using the SEM by a proportion determined by the changes in threshold voltages. The parameter, k, was determined by fitting the theoretical FN curves provided in the references 18,19 onto the experimentally obtained FN plots. The signals for the FN plots were taken from the (310) type facets of a tungsten tip, as indicated by the rectangles in the inset of Fig. 4e. There are three fitting parameters: temperature, work function and k. The previously obtained value of 4.45 eV was used for the work function 20 . The temperature was set to room temperature or 300 K. The k parameter was determined by a curve fitting using the FN plots. Thus, the obtained k value for the (310) type facets was multiplied by a factor of 1.04 to obtain k top 7 . The uncertainty of F DC determined in this method is 20 percent at most 7 .  Fig. 4a, a molecule emission pattern marked by a dashed circle can be observed at room temperature, and the buffer layer can also be observed as the darker area spreads below the dashed line. The surface condition in this case was the same as that in image III of Fig. 2. With increasing heating power, the marked molecule pattern disappeared, as shown in Fig. 4b. Further increasing the heating power caused the substrate signal to increase, and the difference between the upper and lower sides also disappeared, as shown in Fig. 4c and d, which indicates the evaporation of the buffer layer from the apex. We also have estimated the temperature when the molecule patterns and the buffer layer disappeared. The temperature at the tip apex was extracted from Fowler-Nordheim (FN) plots in Fig. 4e. The signals of the FN plots were taken from the emission from (310) type facets of a clean tungsten surface as indicated by rectangles in the inset of Fig. 4e. The FN plots change with increasing heating power because of their dependency on temperature 18,19 . From the power-dependent FN plots, we successfully extracted the sample temperatures and plotted them in Fig. 4f. In addition, we recorded the heating power when the molecule patterns and the buffer layer disappeared. We repeated the same experiments a couple of times. The ranges of the resulting powers are indicated by bands in Fig. 4f, and we determined their temperature from the plots. The results show that the molecular patterns disappear at around 640-660 K (green bands), and the buffer layer disappears at around 1050-1120 K (red bands). The previous work showed that the evaporation temperature of the weakly-bound layers was around 570 K 24 . So the molecules, the sources of the peculiar emission patterns, are more or less as weakly bound to the buffer layer as the weakly-bound layers in image IV of Fig. 2. In contrast, the evaporation temperature of the buffer layer is roughly two times higher than that of the weakly-bound layers.

Scientific
Simulations: quick overview of characteristics of dipoles. In order to have a clearer vision of the molecule configurations on a tip, we have simulated optimized molecule configurations on a tip under strong DC electric fields. The details of these interactions and parameters used in the simulations are given in the Method section. In the simulations, we ignored the atomic structures of C 60 and regarded it as a simple sphere. Under DC electric fields, C 60 will be polarized and generate a dipole. The interactions associated with the dipoles govern the molecule configurations under strong electric fields. We shall first examine a couple of characteristics of dipoles. First, as shown in Fig. 5a, if two dipoles are horizontally aligned, they will undergo repulsive forces from each other. In contrast, if they are aligned vertically, they will become attracted to each other. These aspects indicate that molecules at the same height will tend to stay apart from each other, and they prefer to be aligned  The horizontal axis is 1/V tip . The temperature of the tip was extracted by fitting theoretical FN curves given in the references 18,19 onto experimentally obtained FN plots at each heating power. There are three fitting parameters: a temperature, a work function and a voltage conversion factor C L , where the applied voltage V tip can be converted to a local electric field F L via F L = C L V tip . The unit of F L is V/nm. As described in the caption in Fig. 2, C L is expressed by C L = 1 kr . In this series of experiments, the tip radius r was 165 nm. The parameter k was determined by a curve fitting using the FN plots at room temperature. The resulting value of k was 6.91 and hence the value of C L was 8.77 × 10 −4 . For the work function, the previously obtained value of 4.45 eV was used 20 . Then the temperature was determined for the FN plots at each elevated heating power via the curve fittings. The resulting FN curves are shown as solid lines in figure (e) by using the same colour code as those of the corresponding FN plots. (f) Plots of temperature versus heating power. The temperatures were extracted from the FN plots from both rectangle areas in (e) and were calculated by taking an average of the two cases. The width of the colour bands indicate the ranges of heating powers where the molecule patterns (green bands) and the buffer layer (red bands) disappeared. www.nature.com/scientificreports/ vertically if dipoles are induced in them. Another important factor is the total energy of a dipole 25,26 . This energy draws a dipole, or a molecule, to a position with higher electric fields 27 . We shall see the characteristics of the electric field distributions. Around the tip, strong DC electric fields radially spread from the apex along their lines of field force, as illustrated by the green lines in Fig. 5b. The induced dipoles will modify the electric fields around the apex. For instance, we shall assume four molecules on a tip apex as drawn in Fig. 5b. The electric fields at position A will be reduced by the dipole fields in the two molecules beside it because the signs of the dipole fields are opposite to those of the DC electric fields. In contrast, the electric fields are enhanced at position B because the sign of the dipole fields is the same as that of the DC fields. This means that the single molecule on a buffer layer has stronger electric fields than the molecules in the layer beneath. Because of the characteristics mentioned above, molecules would tend to form a protrusion under strong DC electric fields, and the protrusion becomes the potential source of electrons as it will have high electric fields. It should be noted that the normal explanation of field enhancement above a surface protrusion involves charging of the protrusion 28 . In this simulation, since the charging effect was ignored, the enhanced electric fields at a protrusion are supposed to be underestimated. To our minds, further electric-field enhancement at a protrusion due to the charging effect would further stabilize the protrusion by enhancing a dipole moment at the protrusion molecule. We believe that the further stabilization would not significantly affect our conclusion of the formation of a single molecule protrusion below inferred from our simple dipole model. In order to achieve more quantitatively accurate simulations, the development of simulation models based on density functional theory is required.
Simulations: optimized molecule structures under DC electric fields. In the simulations, we have randomly deposited 500 fullerenes on a tip apex with its curvature at a radius of 160 nm and formed three layers. Here, we assumed the atomic structures of tungsten for the model tip, and the tip apex is oriented towards a [001] direction. An optimized structure under weak DC electric fields, F DC of 0.5 V/nm, is shown in Fig. 6a. The colour of the upper panel represents the height of molecular centres (sphere centres) from the surface of the sphere of the tip apex, Ds. The colour of the middle panel represents the strength of the electric fields at the centres of the molecular spheres. In the lower panel, these electric fields are plotted against Ds. One can quickly confirm that the molecules are piled up to the third layer in the lower panel. We have simulated optimized structures with increasing DC electric fields step by step. At an F DC of 2.6 V/nm, in Fig. 6b, four molecules jumped up to the fourth layer and formed protrusions. As can be seen in the middle and lower panels, the protrusions had stronger electric fields than other molecules beneath. Those molecules with high electric fields would emit electrons, which would be the same situation as image IV of Fig. 2. At an F DC of 3 V/nm, in Fig. 6c, most of the molecules above the first layer evaporated due to the wind force, which is a force caused by the field emission from a molecule. This is the evaporation of the weakly-bound layers, which was observed in the experiments. Note that the weakly-bound layers evaporated around 3.0 V/nm in the experiments, as shown in Fig. 3. The agreement is remarkable, though the wind force is a very simple model 29,30 . For a more accurate descriptions of the physics, one has to simulate using a density functional theory 31 .
After the evaporation of the weakly-bound layer, some molecules remained on the first layer, as shown in Fig. 6c. Those molecules on the first layer were separated from each other and formed single-molecule-terminated protrusions. At even higher electric fields, in Fig. 6d, most of the molecules evaporated from the second layer, and only one molecule remained. Those protrusions in Fig. 6c and d have higher electric fields compared to others

Discussion
We would like to point out configurational differences around single-molecule protrusions. As can be seen in Fig. 7a, some single molecules of protrusions evaporated around 3.5 V/nm and some remained at up to 5 V/nm. The difference in the electric fields for evaporation originates from the molecule configurations. For instance, one of the molecules in Fig. 6c, which remained even at the higher electric fields in Fig. 6d, is marked by a red circle and shown in its front view in Fig. 7b. The difference between the marked molecule and the others is the alignment of molecules in the first layer, as schematically drawn in Fig. 7c. Those molecules which evaporated around 3.5 V/nm are among case I, where the molecules in the first layer are closely located to each other. Case I includes configurations such as a molecule sitting on a hollow site, bridge site and on-top site. In contrast, the molecules that remained at the higher electric fields are in case II, where the molecules in the first layer had some spaces. Note that this configuration resembles the previously observed C 60 island with apexes 32 . These two cases are gauged by Ds, since Ds in case I is higher than that of case II. This tendency can be clearly seen in Fig. 7d. The average height of the topmost molecule is decreasing with increasing F DC . Thus, the value of F DC needed to induce molecule evaporation depends on the height of the topmost molecule above the surface, with low height corresponding to high F DC , and vice-versa. The molecules in case I evaporated first because (at all applied voltages) the local electric fields at the topmost molecules were higher in case I than in case II, due to the molecule configurations. When the local electric field is sufficiently high, then the local field emission current is sufficiently high, the local wind force is sufficiently high, and evaporation occurs. It should be emphasized that there are exceptional cases where protrusion of the case I group remains even at higher electric fields. As marked by a red circle in the upper panel of Fig. 7e, such a protrusion was situated at the edge of the molecule layer. In this protrusion, a molecule was situated at an on-top site, forming a vertically aligned dimer. As shown in Fig. 5b, if a molecule is surrounded by other molecules at the same level, the electric fields at the centre molecule will be reduced due to the dipole fields from surrounding molecules. This effect will be mitigated at the edge of a layer. This aspect can be clearly seen in the lower panel of Fig. 7e, where it shows plots of electric fields versus Ds. The two molecules of the protrusion are marked by red circles. The electric field of the marked molecule in the first layer was higher than the others in the first layer. This high electric field created a stronger attraction force in dipole-dipole interaction. As a result, it could remain even in higher electric fields.
In the end, we would like to comment on the emission patterns. As discussed above, the protrusions had different configurations. All such protrusions are considered to be the potential source of electrons. Here, a question arises: Are the molecule emission patterns influenced by the difference in the molecule configuration in the first layer? We believe that the answer is Yes. These emitted electrons were originally coming from the substrate metal, and they passed through molecules before emission. The emission patterns would reflect molecular orbitals that electrons pass through. Therefore, for instance, if a topmost molecule was at an on-top site, the emission pattern from that topmost molecule would reflect the molecular orbitals of two molecules in the first and the second layers. In contrast, the emission from the topmost molecule in case II in Fig. 7c would be less affected by the www.nature.com/scientificreports/ molecule in the first layer because there was no molecule just beneath the topmost molecule. In this case, the molecular orbitals in a single fullerene could be investigated. In this article, we have characterized the apex of a molecule-covered tip at different amounts of molecules and DC electric fields from an FEM. At the high amount of molecule deposition and under DC electric fields below 3 V/nm, weakly-bound molecules covered the buffer layer. At the same deposition amount and under stronger DC electric fields, the weakly-bound molecules were evaporated, and we confirmed that the molecules, which were the source of the peculiar emission patterns, appeared on a molecule buffer layer. The simulations on optimized molecule configurations on a tip reproduced the same phenomena. Also, they revealed that molecules on the buffer layers were mostly spatially separated from each other, forming single-molecule-terminated protrusions under strong DC electric fields. The simulations further showed the possible configurations of protrusions. These results are helpful to construct a theoretical model to understand the mechanism of electron emission from a single molecule on a tip. FEMs have been applied for investigating molecules many times over the past seven decades. Our study finally showed that an FEM is a simple and useful tool to extract electron signals from a single molecule. One does not need to create a special condition like an ultrasharp tip to realize an electron emission from a single molecule. This conclusion can be applied only to fullerene molecules in this study, but the same conclusion could be drawn for other molecules as well because the appearance of the emission patterns is the same among other molecules. Our results would be useful in reveal the interpretation of peculiar emission patterns and could pave the way to rekindling the FEM as a powerful tool to study a single molecule.

Methods
Experimental. Field emission experiments were performed in an ultra-high vacuum chamber. The base pressure was 1 × 10 −10 mbar. A tip and a phosphor screen were installed in the chamber to perform FEM experiments, as shown in Fig. 1. A metallic mesh was installed in front of the phosphor screen. In order to induce field emission, the mesh and phosphor screen were grounded, and the tip was negatively biased except in the experiments depicted in Fig. 4. In the series of experiments shown in Fig. 4, the tip was heated by heating a hairpin wire on which the tip was placed, and the phosphor screen and mesh were positively biased to generate field emission from the tip. A chevron-type microchannel plate was installed between the mesh and the phosphor screen for the low DC-electric-field experiments depicted in image IV in Fig. 2 and the series of experiments shown in Fig. 3. As shown in Fig. 1, fullerene molecules were enclosed in an evaporator boat with a hole, which was installed beneath the tip. Molecules were evaporated from the hole by resistively heating the boat. A thermocouple was used to monitor the temperature of the evaporator boat. A voltage of − 500 V was applied to the tip during the deposition of molecules. Note that we have also attempted the deposition at a positive voltage. With positive voltages, we observed fewer molecular patterns in an FEM image. This tendency was pronounced with increasing voltages. For instance, after deposition at a positive voltage of 1500 V, even the formation of a buffer layer could not be observed, while we could observe many molecular patterns after deposition at − 1500 V. According to the theory of a gas supply process in the field ion microscope (FIM), the polarized gas molecules are attracted to the shank region near the apex and then move towards the tip apex by some sort of hopping motion 27 . In our case with positive voltages, we believe that molecules could be deposited on the shank region, but they could not move towards the tip apex because of the weak electric fields around the tip apex compared to those in FIM.
In the experiments, we used a tungsten tip with a crystallised apex oriented towards the [011] direction. The surface condition of the apex can be assessed using FEM 33 . A clean tungsten tip surface can be obtained by heating the sample. The inset of Fig. 3e shows a typical FEM image of the clean tungsten tip. In the experiments, we frequently repeated molecule deposition on the tip and removed the molecules by heating the tip. The process converted the tungsten surface to a tungsten carbide surface, resulting in a typical FEM image, as shown in image I in Fig. 2. The tungsten carbide could be removed from the surface by annealing the tip under an oxygen atmosphere. However, because there were no remarkable differences in the appearance of molecular patterns between both surfaces, we mainly performed experiments with a tungsten carbide tip. The radius of curvature of the tips used in this study ranged from 150 to 250 nm.
Theoretical. Optimized molecular structures on a tip under strong DC electric fields are simulated by minimizing the energy of a molecular system. Here we included five different kinds of energies, along with a force. They are: (1) interaction energy between C 60 and C 60 ; (2) interaction energy between C 60 and a tip surface; (3) dipole-dipole interaction energy among molecules; (4) interaction energy between image dipoles in the metal and dipoles in the molecules; (5) total energy of a dipole; and (6) a wind force. The third to sixth terms appear under DC electric fields. The last term is a force exerted on a molecule when electron currents are induced at the molecule due to field emission 29,30 . First and second energy terms. The first two terms are calculated based on a density-functional theory using Quantum Espresso 34,35 . In the simulation, the non-local van der Waals functional (vdW-DF2-B86r) was applied 36 . In Fig. 8a, the resulting interaction energy of C 60 -C 60 is plotted as a function of their distance D. Each C 60 has its dimer face towards the other as schematically drawn in Fig. 8a. In the optimization simulation, we determined the interaction energy of C 60 -C 60 by interpolating the resulting plots in Fig. 8a. For the interaction energy of the C 60 -tungsten surface, we first simulated energy variations with different positions of C 60 with respect to the substrate atomic configurations to introduce the surface atomic corrugations into the simulations. The actual geometry of tungsten atoms is shown in the inset of Fig. 8b. The tungsten surface is orienting toward [001]. In Fig. 8b,  , where x is the distance between a C 60 sphere and a tungsten atom. In the experiments, the tungsten carbide surface was mainly used instead of tungsten. So we modified the obtained potential formula to mimic the experimental situation. From the heating evaporation experiments in Fig. 4, we have learned that the evaporation temperature of the buffer layer is approximately two times higher than that of weakly-bound molecules. The interaction energy for the weakly-bound molecules is roughly estimated by the case of a tetramer, as shown in Fig. 8a. The solid black lines indicate the interaction energy of the tetramer as a function of D. Then, the blue energy curve obtained for the hollow site in Fig. 8b was scaled down in such a way that the potential depth of the curve becomes double that of the tetramer curve, which is shown by the black curve in Fig. 8b. The scaling factor to the original blue curve was 0.2. Hence, we used the following formula for the interaction energy between C 60 and a substrate atom: V LJ2 = 0.8ǫ 0 The strength of the interaction energy between C 60 and a substrate atom was not important to our conclusion as long as it is significantly stronger than those of weakly-bound molecules so that all molecules do not evaporate at the low DC electric fields. Note that here we assumed the atomic structures of the tungsten (100) surface for the tungsten carbide surface and assumed all the substrate atoms were the same kind because the detailed structure of the tungsten carbide surface is not known. This factor would not very important as long as the atomic corrugation was introduced. Without the atomic corrugation of a substrate, molecules in the buffer layer can be slid in directions parallel to the substrate surface, no matter how strong the interactions are. The introduction of the surface atomic corrugation was to avoid such an unrealistic behaviour.
Third, fourth and fifth energy terms. The dipole-dipole interaction energy, V dip−dip for the third term is calculated using the following formula 37 : , where ǫ 0 is the dielectric constant, r i and r j are the position vectors for the two dipoles, r ij is given by r ij = � r i − � r j and µ i is given by � µ i = α � F i . F i is the electric field at the position r i , and α is a molecular polarizability. We employed the molecular polarizability of 106.3 meV/(V/nm) 2 (i.e. α/4πǫ 0 = 153.1 Å 3 ), which is calculated in the reference 38 . Note that meV/(V/nm) 2 is a polarizability unit compatible with the modern ǫ 0 -based equation system used to define SI units, whereas Å 3 is a unit of Gaussian polarizability, as used in the Gaussian-cgs equation system 39 . V dip−dip was calculated among the dipoles in molecules around the tip apex.
For the fourth term, namely the interaction energy between the image dipoles in the metal and dipoles in the molecules, we simulated the positions and amounts of classical image dipole charges for a metallic sphere. The image of an electric point dipole consists of not only an image dipole but also an image charge 40 . Suppose that there is a dipole with a dipole moment of µ i at r i on the tip surface, the image dipole and charge in the metallic sphere are located at � p im = R 2 surf r 2 i � r i , where r i is − → r i and R surf is the distance between the tangential image plane and the center of the sphere. R surf is given by R ′ surf + 1.37 Å in our simulations, where R ′ surf is the distance between the centre of the metallic sphere and the nearest substrate atom of the dipole at r i . Here, we assume the same arrangement of the substrate atoms as mentioned in the previous section. The small factor of 1.37 Å is added to account for the fact that the effective image plane is slightly away from the outermost substrate atom. The value of 1.37 Å was taken from a previous study, wherein it was determined for a tungsten substrate 41 . The induced dipole moment is given by � , and the amount of the induced charge is given by www.nature.com/scientificreports/ q im = R surf r 3 i � r i · � µ i . The interaction energy between the original dipole and the image dipole was calculated by using the abovementioned formula for V dip−dip . The interaction energy between the original dipole and the image charge is calculated using the following formula 37 : V charge−dip = q im 4πǫ 0 r 3 cd � r cd · � µ i . r cd is defined as � r cd = � p im − � r i , and r cd = |� r cd | . Using V dip−dip and V charge−dip , we calculated the interaction energies between the original dipole and image charges/dipoles that include the images induced by other dipoles at the tip apex. It must be noted that the interaction between the image dipoles and dipoles at the molecules in the second layer is negligibly weak and does not affect the formation of the protrusion molecules. Hence, this fourth term representing the interaction energy is not crucial to our conclusion.
It should also be mentioned that the image-related dipole-dipole term is around two orders of magnitude stronger than the image-related charge-dipole term for our simulation conditions. Now, we compare our imagerelated dipole-dipole term with the formula presented in previous works 25,42,43 . A dipole on the surface and its image in the surface are necessarily on the same vertical line, so our image-related dipole-dipole formula appears to reduce to the following: V dip−dip = − µ 2 16πǫ 0 z 3 , where z [ = r ij /2 ] is the distance of the dipole centre from the image plane (which is assumed to be half the distance between the dipoles for consistency), and μ is the common dipole moment. However, the formula that appears in references 25,42,43 is equivalent to V dip−dip = − µ 2 32πǫ 0 z 3 . This implies that the predicted local bonding energy is smaller by a factor of around 2. Despite this discrepancy, we highlight that the forces calculated by both formulae are the same. For convenience, we call our formula V1 and the formula obtained in the previous studies V2. If we move a dipole with a small displacement of + d along an axis perpendicular to the surface, r ij in the abovementioned V dip−dip formula becomes approximately 2z + 2d. Similarly, for a small displacement of -d, we obtain r ij ≈ 2z − 2d. Hence, the differentiation of V 1 can be approximately written as V1(z+d)−V1(z−d) 4d . In contrast, the differentiation of V2 is V2(z+d)−V2(z−d)

2d
. Therefore, upon differentiation, our formula ends up providing the same results as those provided by the formula in references 25,42,43 . As explained below, since we evaluated mean forces exerted on molecules in our simulations, both formulae should provide the same results.
For the fifth term, the total energy of a dipole is given by the following formula: Sixth term: the wind force. The wind force is given by the following formula: � F wind = en e ρ a n a � J , where e is the elementary charge, n e is the electron density, n a the density of impurity atoms, ρ a the impurity resistivity, and J the current density of field emission from a molecule. For n e , we calculated the value for tungsten with its Fermi energy of 9.2 eV 44 , which is 1.27 × 10 23 /cm 3 . For n a , we took the value for a solid C 60 , which is 1.44 × 10 21 /cm 3 45 . For ρ a , we took an experimentally obtained value from the reference 46 , which was 200 cm . To calculate the current density of the field emission, we employed the work function of 5.0 eV. This is because the previous work revealed that the work function of C 60 -covered metal surfaces was about this value for various metal surfaces 21 .
Explanation of the algorithm of an optimization process. Using the six terms above, we calculated optimized molecule configurations. First, 500 molecules were randomly placed on the apex of a tip with its curvature at a radius of 160 nm. The area of molecules was limited within a polar angle of 3°, where the origin was situated at the centre of the sphere of the tip apex and the polar angle was measured from the tip axis. First, molecules were placed randomly one by one at a height of 5 Å from the tip surface. The distances of any two molecules were set to be more than 10 Å at the initial positions. If a molecule could not find a vacant space in the first layer, then it was placed in the second layer at a height of 15 Å from the surface. The same procedure was repeated until 500 molecules were placed. Thus we could make three layers. Then the molecule positions were relaxed and settled into optimized structures.
In the optimization process, we first, calculated the DC electric fields around the tip apex at each molecule position based on the method used in our previous work 20 . The electric fields induce dipole fields at molecules, which modulate the electric fields around the tip apex. In order to accurately calculate the modulated fields, one has to solve 3 × N simultaneous equations 47 , where N is the number of molecules in our case. Here, in order to accelerate the simulations, we simply integrated the induced and the original fields at each molecule. The resulting fields were similar to those calculated by solving simultaneous equations within a difference of around 10%. We confirmed that there were no apparent differences in the resulting optimized structures between the two methods. The updated electric fields were used to calculate the dipole terms and the wind force.
By calculating all the six terms discussed above and using a conjugate gradient method based on the Polak-Ribiere algorithm 48,49 , the directions in which to move the molecules were determined. The lengths of the displacements were determined by a line minimization process 50 . All the molecules were displaced to new positions given by the calculated directions and lengths. At the new positions, the same procedures were applied to determine new positions, and the process was repeated until the mean forces exerted on molecules decreased below a certain threshold. The threshold was set around 100 times lower than the initial mean forces. Further reduction of the threshold did not dramatically change the resulting molecule configuration. Here, we defined the electric field F DC as the value of the electric fields at the very top of the tip apex. The optimized structures were first simulated under F DC of 0.5 V/nm. Then we simulated the optimized structures by gradually ramping up the F DC such as 2 V/nm, 2.2 V/nm, 2.4 V/nm, … 3 V/nm, 3.1 V/nm, 3.2 V/nm … , 5 V/nm.
In the end, we would like to comment on two parameters whose values are scattering in some ranges among different pieces of literature. One is molecular polarizability. The molecular polarizability for a fullerene is scattering from 52.2 meV/(V/nm) 2 (i.e. α/4πǫ 0 = 75.1 Å 3 ) to 106.3 meV/(V/nm) 2 (153.1Å 3 ) in different works 38,51 . Here, we employed 106.3 meV/(V/nm) 2 (153.1Å 3 ), but we also tested with molecular polarizability of 60.1 meV/ (V/nm) 2 (86.5Å 3 ) 52 . With this lower value, molecules do not change their positions, and the formations of protrusions on a weakly-bound layer, such as the ones discussed in Fig. 6b www.nature.com/scientificreports/ contrast, we observed some spots in the emission patterns below 3 V/nm in our experiments, as shown in Figs. 2 and 3, which indicate formations of some protrusions on the weakly-bound layer. Since the lower molecular polarizability did not reproduce our observations, we employed the highest value. Additionally, the impurity resistivity ρ a of 3.8 cm was calculated for a fullerene molecule in different literature 36 , which is far lower than the value of 200 cm in our simulations. We have also tested this value. In this case, the weakly-bound molecules evaporate at around 3.7 V/nm. This change did not affect our conclusion.