Behavior of critical nuclei of pentacene formed on a substrate surface based on the results of molecular dynamics simulations

In the author’s last paper [Ikeda, Appl. Phys. Express 13, 015508 (2020)], it was suggested that the number of molecules that constitute a critical nucleus of pentacene is about ten based on the results of molecular dynamics (MD) simulations performed to investigate the stability of clusters comprising standing pentacene molecules on substrate surfaces. In this study, the author carried out additional MD simulations of clusters comprising ten pentacene molecules and found that these molecular clusters of critical size show stochastic behavior, which can be understood based on curves of free energy change. The discrepancy in the number of molecules that constitute a critical nucleus of pentacene between previous reports and this study is also discussed.


Introduction
The research field of organic semiconductors is continuously developing based on various factors, including the development of new compounds on the basis of a novel concept for chemical design [1][2][3] as well as the development of fabrication technology of devices. However, in spite of the accumulation of excellent research results, organic field-effect transistors (OFETs) 4,5) seem to be facing a barrier that impedes practical applications. One of the factors in this difficulty is considered to be related to the difficulty in obtaining high-quality thin films essential for fabricating organic thin film transistors (OTFTs), 6) which will probably become a major part of OFETs in practical devices. Thin films have several factors reducing charge carrier mobility compared to single crystals; for example, grain boundaries that trap the carriers. 7,8) Deep understanding of the mechanism of nucleation and growth of organic semiconductor thin films is indispensable to control and improve the quality of the films. In reality, fundamental questions remain unanswered with respect to the nucleation and growth process of organic semiconductor thin films. Even the question of when rod-shaped molecules stand up on substrates is still under debate, though such a "standing" arrangement of molecules is one of the well-known phenomena universally observed in organic thin films.
In the author's previous study, 9) classical molecular dynamics (MD) simulations were applied to this problem. Rod-shaped molecules in thin films tend to stand on inert substrate surfaces such as amorphous SiO 2 when the substrate temperature is relatively high and the growth rate is relatively low. 10,11) Pentacene, among the rest, tends to stand on inert substrates and this phenomenon is explained by the assumption that (001) plane has the lowest surface energy and the system leaves the lowest-energy surface exposed. 12) The mechanism can also be explained universally based on thermodynamics and kinetic analysis. 13,14) However, the dynamic process of this phenomenon is still unclear. The isolated molecules tend to lie on surfaces because of the interaction with the surfaces of the substrates, and most of the molecules adsorbed on the surfaces probably lie flat in the early stage of molecular deposition. 9) Therefore, we can hypothesize that some kind of transition from flat-lying to standing occurs at a certain point in time during the deposition process.
The author first approached this general tendency of the standing manner using MD simulations but found that it was very difficult to observe the process in which molecules spontaneously stood up on substrate surfaces using this approach. As shown in supplementary movie 1, available online at stacks.iop.org/JJAP/59/115506/mmedia and also suggested in Ref. 9, isolated molecules adsorbed on the substrate tend to lie flat, and other molecules deposited later are combined with the existing (flat-lying) molecules and the molecular clusters maintain the flat-lying manner (the deposition rate applied in the simulation shown in supplementary movie 1 was unrealistically high due to the limitation of computational capabilities, but this result can provide evidence that molecules tend to lie flat on surfaces in the early stage of deposition). Judging from the grain size of pentacene generally observed in experiments (∼1 μm), 6,8,11,12) a nucleus probably appears in every area of ∼10 −12 m 2 , which is several orders of magnitude larger than the area the author can treat in MD simulations (∼10 −16 m 2 ). It seemed hard for the author to deal with the nucleation problem (specifically, direct observation of the moment in which molecules stand up spontaneously) using MD simulations. Based on this preliminary result and consideration, the author changed the research strategy and focused on the stability of the clusters composed of standing molecules (molecules that have already stood up as the initial condition of the MD simulations). 9) In Ref. 9, the author showed that clusters comprising 15 or more pentacene molecules stably maintained the standing manner for 100 ps and could grow further by capturing additional molecules, whereas clusters comprising 8 or fewer molecules fell down and lay flat on the substrate surface; these results suggested that the sizes of critical nuclei of pentacene correspond to more than 8 and fewer than 15 molecules (on a hydrophobic surface). Subsequently, in Ref. 9, the author performed MD simulations of clusters comprising 12 standing molecules (three runs) and 10 standing molecules (three runs). In the case of 10 molecules, the cluster fell down in one run and maintained the standing manner in two runs, whereas the standing manner was maintained in all three runs in the case of 12 molecules. Based on these results, the author suggested that the number of molecules that constitute a critical nucleus is about ten. However, there remain uncertainties in this conclusion because nucleation is a stochastic process and three runs of MD simulations seemed to be insufficient to discuss it.
In this study, to discuss the behavior of a cluster comprising ten standing molecules in more detail, the author performed additional three runs of MD simulations using the same conditions as used in Ref. 9 (the basic technique of the MD simulations used in both Ref. 9 and this study was first demonstrated in Ref. 15).

Experimental methods
2.1. Models of the substrate and a pentacene molecule Figure 1 shows the models of the amorphous SiO 2 substrate and a pentacene molecule used for MD simulations in this study. First, a unit cell (1.6 × 1.6 × 1.6 nm 3 ) of amorphous SiO 2 consisting of 90 Si atoms and 180 O atoms was prepared [ Fig. 1(a)]. The arrangement of Si atoms and O atoms in the cell were determined by MD simulations using the force field proposed in Ref. 16; the system was first completely molten (mixed) at 5000 K, and then relaxed at 300 K. The Si atoms on the outermost surface of amorphous Si substrates were replaced by Si-OSi(CH 3 ) 3  Gaussian 09 with the Hartree-Fock method using 6-31 G(d) basis sets. 15) By combining the 64 cubic unit cells shown in Fig. 1(b) [8 cells along the X-direction and 8 cells along the Y-direction as shown in Fig. 1(e)], a model of the substrate with dimensions of 12.8 nm (X) × 12.8 nm (Y) × 1.6 nm (Z) was prepared as shown in Fig. 1(f) [strictly speaking, a box-shaped vacuum space with dimensions of 12.8 nm (X) × 12.8 nm (Y) × 9.6 nm (Z) was attached on the substrate model as shown by the green lines in Figs. 1(e) and 1(f)]. Although another model of the substrate showing its hydrophilic nature was also used in Refs. 9 and 15, only the hydrophobic substrate was used in this study. This is because the molecular clusters were less stable on the hydrophilic surface than on the hydrophobic surface in the results of MD simulations shown in Ref. 9 and it was considered that the hydrophobic surface was more suitable for the discussion of critical nuclei.

Force field and computational conditions
MD simulations have been applied to various problems regarding pentacene at the molecular level; for example, deposition of a pentacene molecule at the step edge of a pentacene crystal, 17) packing structures of monolayer pentacene thin films on amorphous silica, 18) crystal structure (polymorphs) of multi-layer pentacene films, 19) growth of pentacene films on the (001) surface of C 60 crystals, 20) and intermixing at the pentacene-C 60 fullerene bilayer interface. 21) The force fields (sets of potential equations and parameters) expressing intramolecular and intermolecular interactions used in those studies varied from each other, specifically MM3, DREIDING, AMBER, etc, and there seems to be no predominant force field used for the investigation of molecular-level phenomena of pentacene.
In this study, the DREIDING force field, 22) a simple generic set of force equations and parameter values, was used for calculating intramolecular and intermolecular interactions. There are two expressions, the Lennard-Jones potential and the Buckingham potential (the exponential-6 form), for describing intermolecular interactions (van der Waals force) commonly used in the DREIDING force field. In this study, the Buckingham potential was used for calculating intermolecular interactions between pentacene molecules based on the preliminary tests performed by the author, 15) while the Lennard-Jones potential was used for other intermolecular interactions, such as between a pentacene molecule and the functional group that modified the substrate surface. Periodic boundary conditions were applied to all simulations, but the Ewald method was not used for calculating Coulomb interactions. All interactions, including Coulomb interactions between atoms within a cut-off distance (R c ) of 3.0 nm, were calculated and used to solve the equation of motion. In the MD simulations of the deposition of pentacene molecules, R c was decreased to 2.0 nm to reduce the calculation time.
The atomic clusters of the functional group, the part of -OSi(CH 3 ) 3 , placed on the substrate surface were able to freely move during the simulations, whereas the Si and O atoms in the amorphous SiO 2 substrates were fixed (velocity was zero). The shape of each pentacene molecule was also flexibly changed by calculating the intramolecular interactions expressed by the four energy terms: bond stretch, angle bend, torsion, and inversion (out-of-plane) angle. Please refer to Ref. 15 for more details of the mathematical expressions of these interaction potentials (force field). All MD simulations were carried out using SCIGRESS ME, a software package for exploring materials produced by Fujitsu Limited, with parallel computing using two Xeon ® processor cores.
The NTV ensemble was applied to the simulations, and the temperature (related to the average velocity) of the pentacene molecules and the functional group was kept at 370 K. This temperature, 370 K (∼97°C), is a relatively high substrate temperature for pentacene thin film growth. One of the reasons for applying this temperature for MD simulations was that pentacene molecules tend to show a flat-lying arrangement at lower substrate temperatures. 11) Another reason (though related to the first reason) was that reevaporation is one of the important factors that influence nucleation of organic materials including pentacene, 10,13,14,23) and the author expected that reevaporation would effectively occur in MD simulations by slightly increasing the temperature. (As a result, there was no obvious effect of the increase of temperature on the behavior of molecules or molecular clusters in this study.) Since we need to consider the time scale of the motion of hydrogen atoms, such as stretching vibration, which is faster than that of other elements, 24,25) the time of each step to solve the equation of motion in MD simulations was set as 0.2 fs in this study. The results of the remaining three (Runs 4, 5, and 6) are new ones carried out in this study. The random sampling numbers used to determine the initial velocity of all movable atoms in the system in the MD simulations were changed for each run to investigate the stochastic behavior. By looking carefully at the results of these six runs, the behavior of the cluster can be classified roughly into the following three patterns: (a) keeps standing (Run 1 and Run 3) or barely keeps standing (Run 6) (b) broken into some pieces and subsequently falls down (Run 4) (c) falls down but keeps the shape of the cluster with herringbone molecular packing (Run 2 and Run 5). In the case of (a), if additional molecules are deposited near the cluster, the cluster can grow further by capturing the molecules and the "standing" cluster is increasingly stabilized as shown in Fig. 3 (and supplementary movie 3).
In the case of (c), the situation is not as simple. The cluster seems to have fallen, but has not fallen on the substrate surface completely. As shown in Fig. 4, the bottoms of the molecules seem to retain some attractive interaction with the substrate and the cluster shows a distorted shape.
3.2. Discussion 3.2.1. Stochastic behavior of critical nuclei discussed based on free energy change. These features can be theoretically understood using the free energy diagram shown in Fig. 5. According to the classical nucleation theory, 26,27) the change in free energy associated with the formation of a spherical nucleus comprising n atoms/molecules can be expressed as where Δμ is the chemical potential difference between the vapor phase (in the case of vacuum deposition) and stable solid phase (per atom/molecule), f is the factor required to derive surface area from volume, and γ is the surface energy density. The first and second terms originate from the volume and surface area of the nucleus (atomic/molecular cluster), respectively. Equation (1) represents a curve being upwardly convex with a value that starts from the origin, increases with increasing n (or radius r of the nucleus), and decreases after exceeding the maximum point (this vertex corresponds to the critical nucleus with the critical radius r c ). If the number of atoms/molecules is small, the effect of the second term (effect of the surface energy which prevents solidification) is predominant and such clusters tend to break into smaller pieces. However, if the radius r of the cluster happens to exceed r c , the first term of Eq. (1) (volume ∝ r 3 ) surpasses the second term (surface area ∝ r 2 ) and a decrease in free energy is caused. As a result, crystal growth begins by capturing further atoms/molecules so that the system can reduce its free energy and be stabilized.
In the case of thin film growth, the volume is not simply proportional to r 3 because films are not spheres. However, the tendency is almost similar to the case of a spherical nucleus and a peak (upwardly convex) appears in ΔG versus n diagrams as shown in the appendix (Fig. A·3). In Fig. 5, the free energy diagram obtained in the appendix under the assumption Δμ = 1.53 × 10 −20 J molecule −1 (corresponding to 3 k B T at 370 K) has been selected and shown. Although the peak position and absolute values are changeable based on the calculation assumptions, the free energy of the standing molecules (blue curve) is placed higher than that of the flatlying molecules (red curve) in the range of a small number of constituent molecules (less than several tens of molecules). A similar tendency has been reported in a previous theoretical study using more detailed models, 14) though the peak position relationship between the standing and flat-lying cases was slightly different from the relationship seen in Fig. 5.
Concerning the effect of Δμ, for example, Ref. 14 also suggested that the nucleation rate suddenly increases when Δμ exceeds 3.7-3.8 × 10 −20 J molecule −1 . In this study, the curves of ΔG calculated under the assumption Δμ = 1.53 × 10 −20 J molecule −1 was used for discussion about the behavior of the molecular clusters (Fig. 5), and this means that the discussion in this study is performed in the range for which Δμ is lower than the critical point causing the sudden rise of nucleation rate suggested by Ref. 14.
The cluster comprising ten standing pentacene molecules is considered to be situated at the peak (vertex) of the blue free energy curve shown in Fig. 5. Since the cluster at the peak is thermodynamically unstable, eventually the cluster will move, for example, to the right or left side of the peak or down to the lower (red) curve to decrease its free energy. Therefore, the three patterns of behavior (a), (b), and (c) observed in the MD simulations described in Sect. 3.1 are consistent with the pictures (arrows (a), (b), and (c) shown in Fig. 5) predicted based on thermodynamics. In the case of pentacene, the clusters that hold ten standing molecules (the state with the maximal free energy) are at the turning point in their lives; some of them will grow further by capturing additional molecules and become more stabilized along arrow (a), some will be broken into fragments and then fall down along arrow (b), and the rest will fall down and provisionally be stabilized in a flat-lying manner along arrow (c) in Fig. 5.
As mentioned before, path (c) in Fig. 5(c) is not as simple. The cluster seems to have fallen, but has not fallen down completely. As shown in Fig. 4, the bottoms of the molecules seem to retain some attractive interaction (bonds) with the substrate and the shape of the cluster is different from those lying completely flat, which can be seen in supplementary movie 1. There is probably a potential barrier between the state of "standing" (the peak of the blue curve in Fig. 5) and the state of "flat-lying" (the peak of the red curve in Fig. 5). The molecular cluster needs to break the bonds between the bottoms of the molecules and the substrate, and the cluster is trapped in the intermediate state as shown in Fig. 4 before completely overturning. The height of this potential barrier probably corresponds to Δγ 1 as shown in Fig. A·1 (in appendix), the advantage of energy that the system obtains by reducing the exposed surface, and the energy can be estimated as of the order of ∼1 × 10 −21 J per molecule (∼1 × 10 −20 J for the cluster comprising ten molecules) if the assumptions (parameters) used in the appendix are applied. This value is not large compared to the value of ΔG; however, this may act as an energy barrier to prevent the overturning of the molecular cluster. 3.2.2. Size of a critical nucleus of pentacene. In some previous experimental studies on the nucleation and growth of pentacene thin films, precise measurement of island size distribution or adsorption probability on oxidized Si surfaces were carried out, and the number of molecules that constituted a nucleus (the smallest stable cluster) was estimated    Fig. 2 at 100 ps). The cluster has almost fallen down on the substrate surface but still retains some attractive interaction between the bottoms of the molecules and the atoms of the outermost surface. It shows a distorted shape due to this remaining interaction with the surface. The white lines show the level of the surface of the substrate and the green lines show the vacuum space (one of the components of the cell for MD simulations) attached on the substrate model shown in Fig. 1. All the atoms of the substrate have been removed from view to make the pentacene molecules easy to see.
to be four or three by analyzing the experimental data using scaling functions. [28][29][30] The author also observed a "stable" cluster comprising four pentacene molecules that maintained herringbone-like packing rigidly in an MD simulation as shown in Fig. 6 (and supplementary movie 4; this simulation was performed as one part of the study in Ref. 9). However, the cluster was lying flat on the surface and it seemed difficult for it to stand up spontaneously and grow further into a crystalline thin film comprising standing molecules. According to the results of the MD simulations shown in this study, it seems natural for the critical nuclei of organic semiconductor thin films (composed of rod-shaped molecules like pentacene) to be defined as clusters comprising "standing" molecules. If we define the critical nuclei as stated above, the number of molecules constituting a critical nucleus of pentacene is probably about ten (possibly depending on the surface conditions), larger than the number estimated based on experiments. [28][29][30] In those experimental studies, the state of the critical nuclei, standing or flat-lying, was not identified clearly.

Conclusions
In summary, the behavior of the cluster comprising ten pentacene molecules on the surface-modified (hydrophobic) SiO 2 substrate has been investigated in detail by classical MD simulations, and stochastic behavior that can be associated with free energy curves was visually shown. The cluster situated at the peak of free energy sometimes remains standing and has a chance to grow by capturing additional molecules to be stabilized further, sometimes falls down after breaking into some fragments, or sometimes falls down without breaking (but an intermediate state may exist before reaching the state of complete rollover). The number of molecules that constitute a critical nucleus of pentacene is about ten if we define the critical nuclei as the "standing" ones. These visual results of MD simulations help us to understand the realistic figure of "critical nuclei" of organic semiconductors, though there still remains a significant question of how such clusters comprising standing molecules can appear on surfaces spontaneously. The author believes that more powerful computation capabilities for MD simulations, which can deal with greater time and size scales, are necessary to visualize the stochastic phenomenon in which clusters comprising more than ten standing molecules appear by chance on a computer display.
If the nucleation process becomes clear, the subsequent growth process can be observed in situ using microscopy systems appropriate for organic materials, 31) and the whole process of nucleation and growth is expected to be fully understood by the combination of theoretical and experimental studies including simulations. The author believes that these efforts to elucidate the microscopic mechanism of nucleation and growth of organic semiconductor materials will contribute to developing organic electronic devices such as high-performance OTFTs through the fabrication of highquality thin films.

Acknowledgments
This research was supported by JSPS KAKENHI Grant Number JP15K04674 and JP18K05253, Japan, and the World Premier International Research Center Initiative (WPI), MEXT, Japan. The author would like to thank an anonymous referee and the associate editor for their helpful and constructive comments which greatly improved the manuscript.

Appendix
In this appendix, simple models to derive the free energy change associated with the formation of a molecular cluster (nucleus) composed of rod-shaped molecules on a substrate surface are proposed and the calculated results based on the models are shown. The contents are basically reproduced, with some modifications, from the supplementary data attached to "open access" articles written by the author (Ref. 9) which may be reused under the terms of the CC BY 4.0 license.
As described in the main text (Sect. 3.2.1), the free energy change associated with nucleation is generally treated based on classical nucleation theory assuming spherical nuclei. 26,27) However, in the case of thin films of organic semiconductors such as pentacene, nucleation generally occurs as a monolayer island, and therefore the increase of the volume is not simply proportional to the cube of the size (radius) of the island. Some modification is needed to discuss the nucleation of molecular semiconductors. Here, the author assumes that the shape of a molecule can be expressed as a square prism for which the length of the side of the base is l and the height is al (a > 1) as shown in Figs. A·1 and A·2. The author also assumes that the molecular cluster (nucleus) has the shape of a square prism comprising n molecules. These models for standing (normal) and flat-lying (lateral) manners can be thought of as those simplifying the models Kubono and Akiyama 14) formulated in their theoretical study in which contact angle was introduced as one of the parameters to express the shape of nuclei comprising flat-lying molecules.
In the case of the standing manner shown in Fig. A·1, the change in free energy of the formation of a monolayer nucleus can be expressed as ( ) ( · ) G n n nl a n l nl 4 A 1 where n is the number of molecules, Δμ is the chemical potential difference (per molecule) between the vapor phase (in the case of vacuum deposition) and the stabilized solid phase, γ 1 and γ 2 are the densities of surface energy of the individual faces, and Δγ 1 is the decrease of surface energy density caused by the interaction (chemical/physical bonding) with the substrate surface (the advantage of energy that the system obtains by reducing the exposed surface). The first term is proportional to the volume of the nucleus and stabilizes the solid state. The second, third, and fourth terms are related to the surface energy. Since the energy density of the side face of the rod-shaped molecule should be larger than that of the upper face because the side faces have bare π electrons, γ 2 can be expressed as bγ 1 (b > 1). The fourth term Δγ 1 is probably not greater than γ 1 , and Δγ 1 can thus be expressed as c 1 γ 1 (0 < c 1 < 1). Therefore, Eq. (A·1) can be transformed into In the case of the flat-lying manner (lateral) shown in Fig. A·2, it is assumed that the shape of the cluster is the same as that of the standing one (a square prism as shown in Fig. A·1) but the cluster is placed on the substrate surface rotated by 90°. In this case, the change in free energy of the formation of a nucleus can be expressed as follows: a = 5.0, b = 1.5, c 1 = 0.3, c 2 = 0.3, l = 4.0 × 10 −10 (m), and γ 1 = 2.0 × 10 −2 (J m -2 ). Figure A·3 show the results of the calculation performed using Eqs. (A·3) and (A·6) and the parameters shown above. Although the shapes of the energy curves largely depend on the assumed Δμ, the free energy is larger in the standing manner (blue curve) than in the flat-lying one (red curve) within the range of a small number of molecules, whereas this relationship is reversed when the number of molecules exceeds several tens (50-60). In the case in which small Δμ is assumed [ Fig. A·3(a)], the free energy continues to increase with increasing number of molecules (at least up to 100 molecules). In contrast, in the case in which extremely large Δμ is assumed [ Fig. A·3(e)], the free energy monotonously decreases without any peak. However, when Δμ is within the range of reasonable values [Figs. A·3(b)-A·3(d)], every curve shows a maximum point (a peak in the upwardly convex shape). These fundamental tendencies are similar to those reported in Ref. 14, except that the peak positions of the curves seem different in the standing case and in the flat-lying case in Ref. 14. Regardless, the fact that similar tendencies can be obtained even in the simpler model probably ensures its validity and universality in this kind of theoretical consideration.