Abstract Title of Dissertation: MODELING AND SIMULATION OF ORGANIC MOLECULAR CLUSTERS AND OVERLAYERS ON SOLID SURFACES Qiang Liu, Doctor of Philosophy, 2011 Dissertation directed by: Professor John D. Weeks Institute for Physical Science and Technology and Department of Chemistry and Biochemistry Driven by the rapid development of experimental methods and technology, nano scale physics and chemistry has become more and more important and practical to study. Monolayers of organic molecules have been studied a lot recently because of many potential applications, such as organic photovoltaic devices (OPV) or organic Liquid Electric Diodes (OLED). It is important to understand and interpret these new experimental advances. At molecular scales, Monte Carlo (MC) simulations and molecular dynamics (MD) are two important methods in computational chemistry and materials science. This dissertation will use these simulation methods along with statistical mechanical theory to study the dynamical properties of single monolayers of organic molecules on solid surfaces. First I give a brief introduction to two-dimensional (2D) molecular systems. Different from bulk system or single molecules, 2D systems have many unique properties, and attract much experimental and theoretical research attention. Some common methods in experimental and theoretical studies are reviewed. After introducing the properties and experimental results of ACA/Ag(111), we build a lattice gas model and run Monte Carlo simulations to help interpret the experiments. The pair approximation, a generalization of mean-field theory, is used to calculate the global phase diagrams and put our model into the more general class of spin-1 Ising models. The pair approximation can be used for modeling various monolayer organic molecular systems which correspond to different regions of the parameter space. Then I studied the C60/ZnPc/Ag(111) system, using molecular dynamic simulations. The C60 molecules form unusual chain structures instead of the close packed islands seen on metal surfaces, and we try to provide a theoretical explanation. Finally I use a density functional theory software to calculate the electronic structures of the C60/ZnPc/Ag(111) systems. This calculation predicts a 0.4e charge transfer from substrate to C60 molecule, believed important for the C60 interactions on these surfaces. In general this thesis studies the behavior of organic monolayers and bilayers on metal substrates. This basic work could help us to understand general 2-D system dynamics and electronic properties, and may help us to find new interesting systems with special properties and applications. MODELING AND SIMULATION OF ORGANIC MOLECULAR CLUSTERS AND OVERLAYERS ON SOLID SURFACES by Qiang Liu Dissertation submitted to the Faculty of the Graduate School of the University of Maryland, College Park in partial fulfillment of the requirements for the degree of Doctor of Philosophy 2011 Advisory Committee: Professor John D. Weeks Professor Ted Einstein Professor Wolfgang Losert Professor Dave Thirumalai Professor Janice Reutt-Robey ii Acknowledgements First and foremost, I would like to thank my thesis advisor, Prof. John D. Weeks, who offered me the opportunity to do my graduate research in his laboratory. His instruction, encouragement and support are important driving forces for my PhD study. It is my honor and pleasure to work under his supervision. I will always remember his patience and kindness. Another special thank goes to Prof. Ellen Williams for being my academic advisor and guidance with my first project. She spent a lot of time going over my report and giving me suggestions and insights which are very valuable for my research. To Prof. Janice Reutt-Robey for her generous help and significant discussions, thank you very much. Without her effort, this thesis would not be possible. Much gratitude goes to two experimentalists, Dr. Chenggang Tao and Dr. Wei Jin. They spent a lot of time doing their experiments, which are the most original inspiration of my research work. It?s always a pleasure to cooperate and discuss with them. I would like to thank my labmates Dr. Jeng-Da Chai, Dr. Jocelyn Rodgers, Dr. Natalia Denesyuk, and Dr. Zhonghan Hu, Mr. Shule Liu and Mr. Richard Charles Remsing. It is my pleasure to work with all of you. I also have to say thanks to all the MRSEC members. I have learned a lot from iii all the presentations and discussions with you. I am grateful to my committee members, Prof. Ted Einstein, Prof. Dave Thirumalai, Prof. Wolfgang Losert, and Prof. Janice Reutt-Robey, for their generous time and helps with my thesis. Last, but not least, I want to say thanks to my parents, Chunjing Liu, Jizhen Zheng and my brother Jian Liu, for their understandings, support and love through all my life. This research is supported by the UMD NSF-MRSEC and its shared facilities under grant # DMR 05-20471. iv Table of Contents Acknowledgement??.. ?????? ? ? ?????????????? ? ii Table of Contents??? ..????????????????????. iv List of Tables????????? ???????????????... ..vi List of Figures????????????????????????? vii List of Abbreviations??????????????????????. x List of Notations???????????????????????? .xi Chapter 1 Introduction??????? .???? .?????? .??. ? 1 1.1 Motivation??????????? ???.?. ????? .???1 1.2 Introduction to surface and steps?? ??????? ?? .???.2 1.3 Organic Molecular films??????? ????????.. ?? 11 Chapter 2 ACA molecules on Ag (1 1 1)?? ?????????.? ??.14 2.1 Introduction??????????? ????? ? ??????14 2.2 Experimental Results and Analysis?? ? ? ?? ? ? ...? ..15 2.3 Theoretical Model of Molecular interactions?? ???.. ..? .??. 25 2.4 Monte Carlo Simulations????????? ?????.?? ?.. 32 2.4.1 Introduction to Monte Carlo Simulation? ???.???.. .33 2.4.2 Simulation details and results??? ????????. ?. 33 2.5 Discussion and Conclusions????????? ???????. ..40 Chapter 3 Pair Approximation: Theory for competition between substrate and intermolecular interactions? ?????????????.. ? 46 3.1 Introduction? ?????????????????????? ..46 3.2 Pair Approximation?? ?????????????????.. ? 50 3.2.1 Pair approximation for two species????.?????.. ? 51 3.2.2 Pair approximation for three species?????? ?. ?..54 3.3 Calculation results?? ?????????????????? ?. 61 3.4 Conclusions and discussions?? ?????????? .???.. 64 3.4.1 Anisotropic situation???????????? ..???? 65 3.4.2 Spin Representation????????????? ..???.67 v 3.4.3 Another method with fixed surroundings???.?????69 Chapter 4 C60 molecules on ZnPc??????????????????74 4.1 Introduction????????? ?????????????. ?.. 74 4.2 Theoretical model of molecular Interactions??????? ??...?79 4.2.1 Van der Waals interactions?????????????..?80 4.2.2 Intermolecular interaction VG: Girifalco Potential???..?.83 4.2.3 Intermolecular interaction Vd: dipole-dipole repulsion??....86 4.2.4 Charge-Charge repulsion???????????...??? .89 4.2.5 Substrate potential???????????????..??.94 4.3 Simulation method: Langevin Dynamics?? ?? ??.. ??.. 96 4.4 Simulation results???????????????????..??97 4.4.1 The substrate effect??????????????.???97 4.4.2 The dipole-dipole repulsion effect??????? ???..?97 4.4.3 Both dipole repulsions and C60-substrate interactions??.?99 4.5 Conclusions?????????????????????...??101 Chapter 5 Quantum Chemistry Calculation (DFT)??????????105 5.1 Comments on density functional theory???? ?????. ??? 106 5.2 C60/ZnPc system??. ???? ?????????????? ?. 108 5.2.1 C60/ZnPc?????????????????????..108 5.2.2 C60/ZnPc/Ag(111)?????????????????....110 5.2.3 C60/ZnPc/Au(111)????????????????...?.117 5.2.4 C60/Ag(111)???????????????? .???.?118 5.3 Discussion and Conclusions????????? ?????? .?? 118 Appendix 1: C60/ZnPc/Ag(111) Quantum Espresso input files?? ..? 119 Chapter 6 Conclusions??????????????????????134 List of References????????????????????? ?? ...136 vi List of Tables Table 2.1 Synopsis of values (1/n, w, G0 and ?c) directly determined from the distribution and correlation functions, and the derived values (??, ?a, ?a) obtained using the limiting values of L=60nm and L=120nm. All values correspond to T =300K????????????????????????????... ?..26 Table 2.2 Parameters determined for MC simulations. A reasonable value of 0.37 eV was chosen for the H-bond energy ???? and the ratio of ???? to ??? = ??? was set to an assumed value of 5. The symmetry balance and step stiffness are mainly determined by this ratio and the difference ?? ??? as given in Eq. 2.10. We constrained our final parameter set to satisfy Eq. 2.10 for square T=0 islands, as suggested by experiment. The density in the gas phase increases with ??, which was chosen to give a high density gas phase at experimental temperatures? ???? .34 Table 3.1 Schematic representation of the effective contribution factors for three species isotropic 2-D square lattice gas. The species circled in the center are species we try to put in, while the surroundings are effective interactions. White sites and shaded sites may have different preference. Z is the normalization factor ??? ??????????????????????????. ?.57 Table 3.2 Schematic representation of the pair approximation method for a three species isotropic 2D square lattice gas. W is the normalization factor?? ??.. ?59 Table 3.3 Schematic representation of the effective contribution factors for three species anisotropic 2-D square lattice gas. The species circled in the center are species we try to put in, while the surroundings are effective interactions. White sites and shaded sites may have different preference. Z and Z? are the normalization factor??? ??????????????????????????.. ?.66 Table 4.1 Parameters for Lennard-Jones interactions? ???????.. ??.82 vii List of Figures Figure 1.1 STM images of a thermally equilibrated step on the Ag(110)???.. 3 Figure 1.2 Island shape with anisotropic interactions in two directions????. 10 Figure 2.1 Structural models and images of ACA on Ag(111)...???????16 Figure 2.2 STM image of ordered islands and surrounding disordered phase? ...18 Figure 2.3 Measurements of boundary fluctuations?? ??????.. ???.20 Figure 2.4 Gaussian distribution of the step position x(t)?? ?????.. ??.21 Figure 2.5 Typical correlation function G(t) and autocorrelation function C(t)? .22 Figure 2.6 Schematic illustration of molecular system with variable substrate interactions as well as intermolecular interactions???? ?.. ??.28 Figure 2.7 Illustration of a rectangular island of ACA?s?? ?? ??????.31 Figure 2.8 Images of Monte Carlo simulations of ACA/Ag(111) system???.36 Figure 2.9 Simulation of G(y) and fits for ????? ????. ???????38 Figure 2.10 Time correlation function G(t) for simulation?? ?????. ??..39 Figure 2.11 Auto correlation function for simulation? ??????.. ????.40 Figure 3.1 Effective interactions for two species 2-D square lattice gas? ?.. ?..52 Figure 3.2 Pair approximation for two species 2-D square lattice gas? ?. ??..53 Figure 3.4 Sub-lattices formed by the U and D molecules??? ?????.. ..55 Figure 3.5 Grand potentials ?G vary with chemical potential ??? ???.. ?..62 Figure 3.6 Phase diagrams with different parameters?? ???????.. ?..63 Figure 3.7 Phase diagram scaled with critical temperature T0? ???? ???64 Figure 3.8 Anisotropic 2D square lattice gas???? ???????.. ???65 Figure 3.9 Fixed surrounding field illustration?? ?????.. ??????.69 Figure 3.10 Fixed surrounding average spins??? ?????.. ??????70 Figure 4.1 C60, ZnPc, and Pentacene structures??? ????. ??????75 viii Figure 4.2 C60/ZnPc/Ag(111) System illustration? ?????. ???? ..?..76 Figure 4.3 ZnPc molecules on Ag (111) surface???? ?? ????? ..?.77 Figure 4.4 C60 meandering chains formed on different organic surfaces? ??.. 78 Figure 4.5 C60 close packed structure on Ag(111)???? ??????. ??78 Figure 4.6 DL_POLY simulation of C60 molecules??? ???????.. ?..83 Figure 4.7 Girifalco Potential and the C-C LJ potential? ??????? ??.85 Figure 4.8 Comparison of Girifalco potential and LJ potential? ??????. .86 Figure 4.9 Stripes formed by dipoles with vdW interactions???? ????. 87 Figure 4.10 Effective C60-C60 interaction: Girifalco potential plus dipole-dipole repulsion? ???????????????????????8 9 Figure 4.11 Error function and complementary error function???? ???... .92 Figure 4.12 Long-ranged and short-ranged Coulomb interactions with ? = 5? ??????????????????????????????. ?.93 Figure 4.13 Forces of short-ranged and whole Coulomb interactions with ? = 5? ???????????? ??????????????????.. ?94 Figure 4.14 Part of the manipulated substrate potential ??????? ?????.. 95 Figure 4.15 Simulation results with substrate potential????? ?????? 97 Figure 4.16 Final configurations with different dipole moments??? ????. 98 Figure 4.17 Simulation results of C60s with dipole d=19.6D, ?=0.16eV? ?. ?... .99 Figure 4.18 Simulation results of C60?s with 0.4e charge, Gaussian truncation ?=50?.Substrate interaction strength ?=0.16eV??? ? ????. 100 Figure 5.1 Initial configurations for DFT calculation of C60/ZnPc?? ???. .109 Figure 5.2 Final configurations of DFT calculation of C60/ZnPc?????. ?109 Figure 5.3 Unit cell for calculation of C60/ZnPc/Ag(111)?? ?????.. ?..112 Figure 5.4 Different positions of C60 on top of ZnPc???? ????? ??113 Figure 5.5 Differential charge density isosurfaces for C60/ZnPc/Ag(111)??..114 Figure 5.6 Differential charge distribution along Z directions for C60/ZnPc/Ag(111) ix system????????????????????????.. 115 Figure 5.7 Differential charge distribution isosurfaces of ??60 ???? ???? ???60 ? ????? ??? ? ??????????????????????. ..116 Figure 5.8 C60/Ag(111) system differential charge distribution isosurfaces?? ????????????????????. ..118 x List of Abbreviations 2D two dimensional ACA acridine-9-carboxylic acid AFM Atomic Force Microscope BEG Blume, Emery, Griffiths DFT Density Functional Thoery EC evaporation-condensation h-BN hexagonal boron nitride H-bond hydrogen-bond LJ Lennard Jones MC Mont Carlo MD molecular dynamics MFT Mean Field Theory ML monolayer NTCDA Naphthalin-tetracarboxylicacid-dianhydride OLED organic light emitting diode OPV organic Photovoltaic OTF organic thin films PWSCF Plane Wave Self-Consistent Field ESPRESSO open-Source Package for Research in Electronic Structure, Simulation, and Optimization rms root mean square STM Scanning Tunneling Microscope vdW van der Waals ZnPc Zinc Phthalocyanine xi List of Notations ? Species in lattice gas model. Same as: ?, A, U, D, G, V ?0 Energy coefficient ?? Effective interactions, Same as: ??, ???, ?? ? Step free energy P4 1 ???? P52 ?? Step stiffness ?0 Energy cost per unit length of a step ? Energy, Lennard-Jones potential constant ?x, ?y Energy bond ? Species in lattice gas model. Same as: ?, A, U, D, G, V P49 Friction constant P96 ??? Effective interactions, same as: ??, ??, ?? ?? Free energy per unit area on solid surface ? Delta function ?? 0 or 1 ? Potential function ? Lennard Jones potential constant ?? Correlation time ? Chemical potential ?(?) Random thermo noise ? External field in lattice gas, same as H ? Hamiltonian ? Area occupied by one atom xii ?( ,),?( ) Gamma function ?? Mobility of the step edge ??(r) Wave function ??, ?? Molecular dimensions A Species in lattice gas model, same as: ?, U, D, G, V ?? Effective interactions, same as: ?? ???, ?? ?(?) Auto correlation function D Species in lattice gas model. Same as: ?, A, U, G, V ?? Effective interactions, same as: ??, ???, ?? erf(?) Error function erfc(?) Complementary error function G Species in lattice gas model, same as ?, A, U, D, V ??(?) Gaussian function ?(?) Space correlation ?(?) Time correlation function H External field in lattice gas, same as ? J Interaction energy, same as: K ?? Boltzmann constant K Interaction energy, same as: J P67 Kelvin P26 L Step length M Number of sites in lattice gas xiii ??,??? 0 or 1 N Number of molecules P Probability r, ?? ???, ?? Distance S Spin ??,?+???,????? Effective spins T Temperature TR Roughing temperature V Species in lattice gas model. Same as: ?, A, U, D, G P50 Volume P33 Potential P84 W Rms of step fluctuation W Normalization factor, Same as Z ?? Momentum space step mode Z Normalization factor, same as: W 1 Chapter 1 Introduction 1.1 Motivation In the past half-century surface and interface phenomena has been studied extensively. Many new and rapidly growing technologies like organic light emitting diode (OLED) or Organic Photovoltaic (OPV) devices are mostly based on recent progress in the study of surface and interfaces. Surface nanostructures are particularly suitable for applications involving sensing and, optical electronics due to their accessibility to environment chemical species and light. To understand the structure and properties of surface is a great challenge and motivated by both technical applications and scientific interest. Recent development of experimental techniques such as scanning tunneling microscopy (STM) and Atomic Force Microscope (AFM) has dramatically enhanced our insight into the surface structures and their structural, chemical and electrical properties. In this thesis, we will try to understand and explain some of these experimental results. Along with the development of experiments, theoretical study of organic molecular films on solid surfaces has become more and more important. Understanding the underlying physics and chemistry at the interface can help us to explore new materials and develop experimental methods and perhaps build new devices. There are several aspects of theoretical study we will use in this thesis: classic molecular dynamics and Monte Carlo simulations, theoretical calculations 2 like mean field theory, and density functional theory (DFT) calculations. The different methods have their own advantages and limitations as will be explained later Computer simulation is a very powerful tool to study variety of processes, in chemical, biological and surface science applications. Simulation methods can be easily adapted in surface systems to study equilibrium and non-equilibrium processes. In the two examples in this thesis, we will use two kinds of simulations, molecular dynamics (MD) and Monte Carlo (MC) simulations, to solve different problems. We will also introduce an analytic method, the pair approximation, as a complement to the computer simulations. This approximation can help us understand the physics more clearly and expand the range of parameters we are studying. In addition we need to know the most basic electronic structures and properties of the interacting molecules and atoms, and here we use density functional theory to study essential aspects of it. Although we can only examine relatively small systems with current computers, the result might give us ideas to apply to larger system models. 1.2 Introduction to surface and steps First, we will introduce some basic concepts and theory for the evolution of surface morphology. A surface is a boundary between two macroscopic regions with 3 different phases. On the surfaces, there usually are a macroscopic number of mobile particles which are the components of one phase. At a temperature below the roughing temperature TR, a crystal surface mainly consists of terraces and steps, and some islands because of excitation on the terraces. Surface morphological changes arise mostly from the motion of these steps. There is a very good review article [1], which uses the well-known continuum step model to study surface evolution. Here we only cite some of the results we will use later for the single monolayer systems. Detailed discussions of most concepts discussed below can be found in this article. In the step model evolution of surface morphology is described in terms of motion of steps, which is directly connected to the attachment and detachment of atoms on the surface. The step motions can be illustrated below in Fig. 1.1, and described using function x(y,t). The surface steps are treated as the fundamental objects. Fig. 1.1 Four STM images of a thermally equilibrated step on the Ag(110) surface. The size of the image is 600 ?? 5000 ?. Each image was acquired in 18 s, and the time increment between images was 30 min. The schematic 4 illustration at the left shows the description of the step as a wandering line which is used in the continuum step model (provided by J. Reutt-Robey). We will define a local chemical potential and use its differences as the driving force to describe the surface dynamics. This is a unified thermodynamic approach that can be applied to most surfaces under near-equilibrium conditions. We suppose that the surfaces we are studying will eventually reach their equilibrium state where the local chemical potential is the same except for fluctuations. The energy cost per unit length of a step for rectangular lattice can be written as ?0 = ??|????|? ? +??|????|? ? (1.1) where ?x and ?y are the energy costs to break a bond, ax, ay are the lattice constants in the x, y directions, and ? is the angle between step tangent direction and y direction. With finite temperature, the steps fluctuate and have a configuration entropy S so that we can calculate the step free energy ? as ? = ?0 ??? (1.2) For a small angle change ? of the step, the change of free energy per unit length due to this distortion is ?? = 12?2(?(0)+???(0)) [2]. We define the step stiffness ??, which is related to the step free energy by ??(?,?) = ?(?,?)+? 2? ??2 (1.3) Then the Hamiltonian of an isolated step in a continuum model can be written as 5 ? = ?? ? 2( ?? ??) 2 ?? (1.4) From this Hamiltonian we can take a functional derivative to get the chemical potential of an isolated step without step-step interactions: ?(?) = ???? 2? ??2 (1.5) where ? is the area occupied by one atom. If we assume the evaporation-condensation (EC) case of mass transportation, then the local step exchanges adatoms with the vapor which serves as a reservoir of constant chemical potential. The local step moves according to this chemical potential difference between the step and the reservoir and the x(y, t) satisfies the equation below [3,4]: ??(?,?) ?? = ???? ??? ?2?(?,?) ??2 +?(?) (1.6) where ?? is the mobility of the step edge, and ?(?) is random white noise term: ??(?)?(??)? = 2??? ?? ?(? ???) (1.7) Often it will be convenient to use the diagonalized form of the Hamiltonian by using Fourier transform ? = 12??? ? 2? 2|??|2?? (1.8) where ?? = 12???(?)???????. The equipartition theorem gives us the relation 6 between different frequency modes: ??????? = ??????2 ?(? +??) (1.9) where the bracket denotes an ensemble average. If we define a mean square displacement along the step edge as, ?(?) ? ?[?(?)??(0)]2? (1.10) to characterize the spacial correlations of step edge positions, then this quantity should increase linearly with y ?(?) = 12??????????????(???? ?1)(????? ?1) = ??2????2(1?cos (??))???2 (???????? ???? 0 ?? ?)? ???? |?| (1.11) We can also study the time dependent fluctuations of a single position along the step edge by defing the time correlation function as: ?(?) ? ?[?(?,?) ??(?,0)]2? (1.12) which is related to the auto correlation function: ?(?) ? ??(?,?)?(?,0)? (1.13) with the relation ?(?) = 2?2 ?2?(?), where ?2 = ??(?,?)2? = ?(0). One simple theoretical assumption is to use periodic boundary condition for a finite positive range of 00 (? > 0) (1.16) which doesn?t depend on y. And the time correlation function is: ?(?) = ?4???????2 (1??? ?????2 ??? ?) ?>0 (1.17) First, the ?2 = ?(0) can be easily calculated as an integral or summation from 2?/? to infinity: ?(0) = ? ?? ? 2? ? 2??? ????2 ? 2? = ???? 2?2?? ?(0) = ? 2???????2 ? 2? ?? = ????12?? (1.18) 8 At long time ? ? ??, the decay of C(t) is governed by the mode with smallest ? = 2?? . This leads to the long time behavior of the auto correlation function [6]: ?(?) = ????2?2?????/?? ?(?) = ?????2?? (1??? ? ??) (1.19) with correlation time: ?? = ???? 2 4?2???? (1.20) Without the ? ? ?? condition, we can take the summation over q into the integral over q. Take this integral from 2?/? to infinity give us the result: ?(?) = ?(0)[?? ? ?? ??(1 2, ? ??)( ? ??) 1 2] (1.21) where ?(0) = ?2 = ????2?2??, and ?? = ????24?2? ??? . At short time ? ? ??, the above equation reduces to Euler gamma function: ?(?) = ?(0)[1??(12)(?? ? ) 1 2] (1.22) And finally we can get ?(?) = 2?? 1 2?(12) ? ( ??? ?? ) 1 2 ?12 = ?0?12 (1.23) For another boundary condition which fixed the two end of the finite step 00 ??????? = 2??????2 ?(? ???) ?? = 1???(?)sin (??) ? ,? = 0,1,?,? ?(?) = ???sin (??) ?>0 ,? = 0,??,2?? ,?,? ?(?) = ???? ?(1???) ?(?,?) = ??(?,?)?(?,0)? = ?2???????2 ???2(??)?? ?????2 ??? ? ?>0 ?(?,?) = ?4???????2 ???2(??)(1??? ?????2 ??? ?) ?>0 (1.24) Here the C(y,t) and G(y,t) do depend on y, because the two ends are fixed. The surface free energy depends on the orientation of the surface and the nature of the crystalline lattice. At zero temperature, the surface energy can be easily calculated by counting the number of neighbor bonds that are broken, because there will be no entropy at zero temperature. At equilibrium, the surface energy should be minimized by changing the area of different surface orientations. The Wulff [7] construction gives the solution and tells us that the equilibrium shape is the envelope of the surfaces perpendicular to the radius of a polar plot of the surface free energy, and the distance of a facet from the center of the equilibrium crystal shape is proportional to its surface energy. 10 ?? = ???????? ??? (1.25) where ?? is the free energy per unit area of the ith facet and ?? is the distance of this surface to the center of the crystal. This result can also extended to two-dimensional island situations: at zero temperature, the distance of one step to the center of the island, should be proportional to the free energy per unit length of that step (Fig. 1.2): Fig. 1.2 Island shape with anisotropic interactions in two directions Thus a material with very anisotropic intermolecular interactions would be expected to form islands with an anisotropic shape at low temperature that would minimize the number of ?broken bonds? involving the strongest intermolecular interactions. This idea plays a key role in our discussion of the ACA/Ag(111) system in the next chapter. 11 1.3 Organic molecules Organic materials that allow electron and hole transport and recombination have extensive materials and device applications [8,9]. Understanding the relationships between the transport characteristics of organic thin films (OTF) and how molecules are structurally arranged within the film is a fundamental and challenging issue. For OTFs, complex molecular symmetries and interactions often cause the formation of multiple density-dependent phases, with pattern size on the nanometer or micron scale. The boundaries of such local ordered regions may dramatically affect the stability, transport and charge recombination properties of the OTFs. To study the organic molecule monolayer structures and dynamics will be the main purpose of this thesis. The adsorption of organic molecules on solid surface is more complex than atomic species, because molecules have irregular shapes and anisotropic intermolecular interactions. Moreover, adsorbed molecules may not be commensurate with the substrate, so they may have different substrate interactions, based on different orientations of molecules themselves and also the position with respect to the substrate. Another complexity comes from the molecular self-assembly. Self-assembly describes the processes in which a disordered system of pre-existing components forms an organized structure or pattern as a consequence of specific, local interactions among the components themselves, without external direction [10]. These self-assembled molecules may have different configurations, which give 12 clusters of the same molecules different properties. Self-assembly is usually through relatively weak interactions, e.g. van der Waals interactions, hydrogen bonds, dipole-dipole interactions. Although they are weak interactions, they play an important role in material synthesis. Charge exchange between two kinds of materials and the corresponding dipoles formation, is one of the essential keys for the self-assembly processes, and also governs the electronic properties of OTF systems, forming the basis for semiconductor devices [11]. Chapter 4 and 5 discuss systems where we argue that charge transfer plays a key role in surface patterns seen in experiment. References: 1 H.-C. Jeong and Ellen D. Williams, ?Step on Surfaces: Experiment and Theory?, Surface Science Report, 34, 171-294 (1999) 2 M.P.A. Fisher, D.S. Fisher and J.D. Weeks, ?Agreement of Capillary-Wave Theory with Exact Results for the Interface Profile of the Two-Dimensional Ising Model?, Phys. Rev. Lett. 48, 368 (1982) 3 N. C. Bartelt, J.L. Goldberg, T. L. Einstein, E.D. Williams, J. C. Heyraud and J.J. Metois, ?Brownian Motion of Steps on Si(111)?, Phys. Rev. B 48, 15453 (1993) 4 C. Dasgupta, M. Constantin, S. Das Sarma, and S. N. Majumdar, ?Survival in Equilibrium Step Fluctuations?, Phys. Rev. E 69, 022101 (2004) 5 S. Kodiyalam, K.E. Khor, N.C. Bartelt, E.D. Williams and S. Das Sarma, ?Energetics of Vicinal Si(111) Steps Using Empirical Potentials?, Phys. Rev. B 51, 5200 (1995) 13 6 O. Bondarchuk, D. B. Dougherty, M. Degawa, E.D. Williams, M. Constantin, C. Dasgupta, and S. Das Sarma, ?Correlation Time for Step Structural Fluctuations?? Phys. Rev. B 71, 045426 (2005) 7 M. Wortis, in Chemistry and Physics of Solid Surfaces VII, eds. R. Vanselow and R.F. Howe, 367-428 (Springer-Verlag, Berlin, 1988) 8 B. A. Gregg, S-G. Chen, and R. A. Cormier, ?Coulomb Forces and Doping in Organic Semiconductors?, Chem. Mater. 16 (23), 4586-4599 (2004) 9 T. Taima, J. Sakai, T. Yamanari, and K. Saito, ?Doping effects for organic photovoltaic cells based on small-molecular-weight semiconductors?, Sol. Energy Mater. Sol. Cells, 93, 6-7, 742 (2009) 10 D. Bonifazi; A. Kiebebe, M. Stohr, F.Y. Cheng, T. Jung, F. Diederich, and H. Spillmann, ?Supramolecular Nanostructuring of Silver Surfaces via Self-Assembly of Fullerene and Porphyrin Modules?, Advanced Functional Materials, 17(7), 1051-1062 (2007) 11 M. A. Baldo and S. R. Forrest, ?Interface-limited Injection in Amorphous Organic Semiconductors?, Phys. Rev. B 64, 085201 (2001) 14 Chapter 2 ACA/Ag(111) System 2.1 Introduction Low-dimensional boundaries between phases and domains in organic thin films are important in charge transport and recombination. Experimentalists have studied the fluctuations of interfacial boundaries in an organic thin film, acridine-9-carboxylic acid (ACA) on Ag(111), using Scanning Tunneling Microscopy. The boundaries fluctuate via attaching and detaching of molecules at the step edge. Although ACA has highly anisotropic intermolecular interactions, but experiments have shown that it forms nearly circular islands that have essentially identical thermodynamic and kinetic properties for different boundaries, in contrast to the anisotropic shape expected as depicted in Fig. 1.2. We argue that the physical basis of the modified symmetry is because of the significantly different substrate interactions induced by alternating configurations of adjacent molecules in the solid phase. All these different interactions can be modeled using a multi-component lattice gas model, and can straightforwardly reproduce the experimentally observed isotropic behavior. The general multi-component description allows the domain shapes and boundary fluctuations to be tuned from isotropic to highly anisotropic in terms of the balance between intermolecular interactions and molecule-substrate interactions. This chapter is heavily based on work that has been published: C. Tao, Q. Liu, B. S. Riddick, W. G. Cullen, J. Reutte-Robey, J. D. Weeks and E. D. Williams, ?Dynamic interfaces in an organic thin film?, Proc. Natl. Acad. 15 Sci. USA 105, 16418-16425 (2008). 2.2 Experimental Results and Analysis The ACA molecule, shown in Fig 2.1a, has a nominally rectangular structure, with the potential for strong hydrogen-bond (H-bond) intermolecular interactions along the short axis, and weaker van der Waals interactions along the long axis. The absorbed ACA molecules do not form ordered solid structures at low coverage (< 0.3 monolayer (ML)) [1]. As ACA coverage is increased beyond 0.3 ML, ordered areas (islands) begin to form, covering the substrate partially. The area around the islands shows no ordered overlayer at room temperature, but noise in the tunneling current suggests the presence of a mobile gas phase of ACA molecules. This is confirmed by measuring these areas at low temperatures, when the motion is quenched (14). Thus above 0.3 ML, ACA exists in two-phase equilibrium between a dense ordered phase and a disordered phase. High resolution imaging [2,3], see Fig. 2.1b, shows that alternating ACA molecules are inequivalent. There are two differently oriented ACA molecules per unit cell. Along the rectangular short edge direction, there is hydrogen-bonding interaction between ACA molecules, with the ring nitrogen acting as the H-bond acceptor and the carboxyl proton acting as the H-bond donor [2], shown in Figs. 2.1c. At condensed solid phase areas, ACA molecules arrange head-to-tail into chains connected by H-bonds between the carboxyl group on one side of the ACA molecule, and the ring nitrogen on the other side. These chains run along the [110] direction of the substrate, as shown in Fig. 2.1b. The asymmetry of 16 the ACA ordered structure provides two distinct types of edge boundaries (parallel to the hydrogen bond chain direction, we call ?S?, and perpendicular to the chains we call ?CE?). (a) (b) (c) Fig. 2.1 Structural models and images of ACA on Ag(111) (by C.-G. Tao). (a) Molecular structure of ACA molecule. (b) measured STM image of ACA on Ag(111). (c) molecular configuration models. Alternating ACA molecules along the head-tail hydrogen bonded chains are tilted with respect to the substrate. The molecule dimensions are ??~2? in the [110] direction and ???~2?3? in the [112] direction where a=0.289nm is the nearest neighbor distance on Ag(111) 17 Because of the simple rectangular ordered structure of the ACA solid phase, the simplest theoretical model we could imagine would be a simple lattice gas (Ising) model, to study the ordered overlayer. The anisotropic ACA interactions would represented by anisotropic nearest neighbor bonds in different directions, with the substrate defining a sub lattice and available three directions for the island growth. In this picture, the CE boundaries would have larger edge energies due to the broken H-bonds (magnitude typically ~0.3 eV per bond) and easy excitation of edge roughness due to the weak vdW interactions (magnitude typically 4 to 6 times smaller than the H-bonds for ACA molecules) along the edge (4,5). The energy/roughness property would be the opposite for the S boundaries. However, the experiment results showed some unusual properties of the island shape as well as the boundaries behaviors. We will explain these by introducing strong non-uniformity of the molecule-substrate interactions for the two molecules in an H-bonded pair due to their different orientations, which might be a general case for structures of molecular overlayers in organic monolayer films. Images of ordered regions formed at an average ACA coverage of ~ 0.6 ML are shown in Fig 2.2. The isolated ordered ACA islands were formed on wide (111) terraces and do not contact silver steps. Such isolated islands are typically compact, with an aspect ratio of the island width along the [110] direction of the substrate to the width along the [112] direction near to one. The aspect ratio ?? ???? of the island, shown in Fig. 2.2, is directly related to the edge formation free-energies ??? and ?? for the two types of steps shown in Fig. 2.1, as ?? ???? = ?? ???? [6,7]. In 18 the simple lattice-model discussed above, the T = 0 edge energy is determined by the interaction energies ??, ??? perpendicular to the boundary edge and the molecular lengths aS,CE along the edge as ??,?????,??/??,?? . Thus the observation of ?? ???? ~1 yields ???/??~?3, in distinct contrast to the expected value of 4-6 for the ratio of the H-bonding to vdW interaction strength. Fig. 2.2 STM image of ordered islands and surrounding disordered phase, the distance from edges to center are ??? and ?? , with aspect ratio close to 1. The edges of ACA ordered islands, i.e. the boundaries between the ordered phase and the disordered phase, as shown in Fig 2.3a, appear frizzy, clearly indicating thermal fluctuations at room temperature [8]. The boundary shown is CE, perpendicular to the [110] direction of the substrate. Based on the molecular arrangement in the ordered phase (see Fig. 2.1), this boundary is formed by ends of the ACA chains in the ordered phase, corresponding to the CE boundary of Fig. 2.1. 19 The inset shows a boundary roughly perpendicular to the [ 112 ] direction, corresponding to the S boundary of Fig. 2.1. For temporal imaging, the STM scan direction is chosen perpendicular to the boundary orientation, as indicated by white arrows. Examples of the temporal pseudo-images formed as the STM tip repeatedly scans over a point at the boundary are shown in Fig. 2.3b. The vertical axis is time t and the horizontal axis is the position of the point x(t). (a) 102s (2000 lines) 50nm (b) 20 Fig. 2.3 Measurements of boundary fluctuations (by C.G. Tao). The white arrows indicate the scanning direction. (a)STM images of the boundary of ordered and disordered phases. The boundary is perpendicular to the [110] direction and the scan direction is along [110] direction (CE boundary); (inset) boundary perpendicular to [112] direction with the scan direction is along the [112] direction (S boundary). (b) Pseudo-images of boundaries fluctuations, with line scan size 50 nm (horizontal axis), line scan time 51.2 ms, and total measurement time 102.4 s (vertical axis) for 2000 lines. (left) boundary perpendicular to the [110] direction with the scan direction along the [110] direction (CE boundary); (right) boundary perpendicular to the [112] direction with the scan direction along the [112] direction (S boundary). The distribution of displacements has a Gaussian distribution, as required for fluctuations in thermal equilibrium. There is no deviation from the Gaussian shape for scans measured in either direction across the boundary, providing strong confirmation that under the chosen tunneling conditions the STM scans are not perturbing the thermal distribution at the boundaries. The mean squared boundary width ?2 = ?(? ???)2? is obtained by fitting the Gaussian, yielding width w = 2.04 ? 0.04 nm for CE boundaries and 2.01 ? 0.08 nm for S boundaries (Fig. 2.4). 21 Fig. 2.4 Gaussian distribution of the step position x(t) (by C.G. Tao). The red solid line is the Gaussian fit for the experimental data, y0 +Ae?( ???0 ? ) 2 , with best fitting parameters, y0 = -0.25 ? 0.38, A = 67 ? 1, x0 = ? 0.06 ? 0.01 nm and w = 1.72 ? 0.02 nm. The free energies and time constants governing the behavior of the boundaries can be evaluated from the correlation functions of the boundary fluctuations. Given the experimentally measured boundary displacements, x(t), it is straightforward to determine the time correlation function G(t) and autocorrelation function C(t) respectively [9]: ?(?) = ?(?(? +?0)??(?0))2? (2.1) ?(?) = ?(?(? +?0)???)(?(?0)???)? (2.2) Each individual x(t) data set is used to calculate individual correlation functions. Fig. 2.5 shows typical results for G(t) and C(t), each the average of the individual 22 correlation functions for more than 10 x(t) sets. Fig. 2.5 (by C.G. Tao) Typical correlation function G(t) (black solid circles, right axis, CE boundary) and autocorrelation function C(t) (black solid triangles, left axis, S boundary) data. The red curve is the fit for G(t), and extracts the best fit values, G0 = 3.86 ? 0.28 nm2, 1/n = 0.52 ? 0.07. The blue curve is the fit for C(t) using Eq. 2.4, and extracts the best fit values, C(0) = 1.81 ? 0.36 nm2 and ?c = 1.77 ? 0.49 s. The continuum step model is used to relate the correlation functions to the physical properties of the boundary. As shown in Fig. 2.5, the measured time correlation function G(t) is well fit using [10]: ?(?) = ?0?1? 23 ?0 = (2?(1?1/n)? )(????? ) ??1 ? ??1/? (2.3) where ?a is the mobility of the boundary, ?(1-1/n) the gamma function and ?? the step stiffness. The fit to the time correlation function, shown as the solid red curve in Fig. 2.5, yields the exponent 1/n = 0.52 ? 0.07. The average value 1/n = 0.52 ? 0.09 for the CE boundary and 0.51 ? 0.05 for the S boundary. These values show that the boundary fluctuations are dominated by uncorrelated events for both boundaries (n = 2) [11], as opposed to correlated events such as diffusion along the boundaries, which would yield n = 4 [12]. The fits also yield the magnitudes of the time correlation functions, which are G0 = 3.76 nm2 for CE boundaries and 3.70 nm2 for S boundaries. Using Eq. 2.3, we can find the ratios of the physical constants governing the edge fluctuations: ??????? ? 11.1??4??1 for CE boundaries and 10.8 nm4s?1 for S boundaries. Since the mobility of a step edge is generally lower when the stiffness is larger, it is surprising to find these ratios similar for the two types of step edges. The significance of these values will be discussed below. Given the result that the boundary fluctuations are predominantly due to molecular exchange (n = 2), we can further analyze the autocorrelation function as shown in Fig. 2.5. For n = 2, the autocorrelation functions have the time dependence: ?(?) = ?(0)[?? ? ?? ??(1 2, ? ??)( ? ??) 1 2] (2.4) ?? = ( ?2?)2???? ??? (2.5) 24 where ?? is the correlation time, ?(12, ?? ? ) is the incomplete gamma function, L is the correlation length (or effective system size), and ?? the mobility of the boundary. The fitting curve, shown as the solid red curve in Fig. 2.5 yields the values for C(0) and ?c. The average value of C(0) is 2.74 nm2 for CE boundaries and 2.42 nm2 for S boundaries. The average correlation times are found to be ?? = 3.59 s for CE boundaries and 2.73 s for S boundaries, substantially smaller than the measurement time of 102.4 s. The edge length of the island boundaries, typically from 50 nm to 70 nm, was used to estimate the range of values for the system size, L, of Eq. 2.5. The lower limit of L (= edge length) is obtained if we assume that the island edge fluctuations are not constrained by the corners of the island. If we assume that the corners completely pin the fluctuations, then we obtain the upper limit of L (= 2 x edge length) [13]. Using the measured time constants and L = 60 nm (120 nm) in Eq. 2.5 gives the values ???? ??? = 0.039 nm-2s (0.16 nm-2s) for CE boundaries and 0.030 nm-2s (0.12 nm-2s) for S boundaries. Combining these values with the values of ????? ?? obtained from the values of ?0 above, we obtain, for CE boundaries, edge stiffness ?? = 39.1 meV/nm (78.2 meV/nm) and edge mobility ?? = 16.8 nm3/s (33.6 nm3/s); and for S boundaries, ?? = 45.5 meV/nm (91.0 meV/nm) and ?? = 19.0 nm3/s (38.0 nm3/s). Within the experimental uncertainties shown in Table 2.1, the values for the two different boundary types are the same: ??? ???? = 1.2?0.3? and ??,? ??,?? = 1.1?0.1? . The measured stiffness values can be used in a simple lattice approximation to estimate the lateral interaction energy (kink energy, ?) at the step edges. The 25 relationship ?? ? ???2? ? ??? ( ?? ?? )(1????(??/???))2 yields estimations ?CE ~ 32 meV (43 meV), and ?S ~ 43 meV (56 meV). The estimated magnitude of ?S is much smaller than the expected value of a hydrogen bond, ~0.3-0.4 eV. In addition, the ratio of ?S/?CE is also much smaller than the expected value of 4-6 for the ratio of hydrogen bonds to vdW interactions. Thus the measured values of the step stiffness are clearly inconsistent with the anisotropy that would be expected for a model based solely on intermolecular interactions between ACA molecules. Also, given the values for the edge mobilities, ??, the average time between molecular attachments/detachment events ?a can be expressed as [9]: ?? = ?? 2?? ?? (2.6) where ?? is the molecular dimension along the boundary direction and ?? the molecular dimension perpendicular to the boundary direction. This yields values of ?a ? 34 ms (17 ms) for CE boundaries and 18 ms (9 ms) for S boundaries. 2.3 Theoretical Model of Molecular Interactions The asymmetry of the chain bonding in the ordered phase, shown in Fig. 2.1, should impose a substantial difference in the behavior of the S and CE edge boundaries. Thus the similarity of the measured thermodynamic and kinetic values, specifically the step stiffness and mobility listed in Table 2.1, for these two boundaries is initially quite surprising, and has very interesting physical consequences. Here we will derive a specific model for the interactions between 26 ACA molecules, modeled as pairwise interactions within a lattice gas framework, and test the model?s ability to reproduce the experimentally measured edge free energies. CE-boundary (parallel to [112?]) S-boundary (parallel to [11?0]) 1/n 0.516 ? 0.087 0.512 ? 0.046 w (nm) 2.04 ? 0.04 2.01 ? 0.08 G0 (nm2s-1/2) 3.76 ? 0.25 3.70 ? 0.72 ?c (s) 3.59 ? 0.30 2.73 ? 0.43 L (limits) (nm) 60 ? 120 60 ?120 a// (nm) 1.001 0.578 a? (nm) 0.578 1.001 ?? (meV/nm) lower limit upper limit 39.1 ? 4.2 78.2?8.4 45.5 ? 12.4 91.0?24.8 ?a (nm3/s) lower limit upper limit 16.8 ? 0.4 33.6?0.8 19.0 ? 2.2 38.0?4.4 ?a (s) upper limit lower limit 0.034 ? 0.001 0.017 ? 0.0004 0.018 ? 0.002 0.009 ? 0.001 Table 2.1: Synopsis of values (1/n, w, G0 and ?c) directly determined from the distribution and correlation functions, and the derived values (??, ?a, ?a) obtained using the limiting values of L = 60 nm and L = 120 nm. All values correspond to T = 300 K. In the simple lattice model discussed in the Introduction, each site on a square lattice is either vacant or occupied by a single ACA molecule. Chain-like structures arise from the anisotropic nearest neighbor interactions in which adjacent molecules 27 in one direction (see Fig. 2.1) have a large favorable energy ??CE due to the H-bonding, while neighbors in the opposite (side) direction have a higher energy ???, with ??? ? ?? due to the weaker vdW interaction. ACA-substrate interactions are the same for all molecules in this model, which maps onto the usual spin 1 2? Ising model with anisotropic ferromagnetic coupling between nearest neighbor spins. As shown above, the experimental observations simply do not follow the intuitive expectations of highly asymmetric island shape and edge stiffnesses that follow from this model. A final, and important, problem with this simple lattice gas model is that it predicts a low density in the disordered phase except very near the critical temperature, in contrast to the large density (~ 0.3 ML) observed experimentally [2,3]. We resolve the apparent discrepancies, while still staying within a lattice gas framework, by taking into account the different molecular orientations. As noted in Fig. 2.1, the strong H-bonds along the chain require alternating tilt-orientations of adjacent ACA molecules. As illustrated in Fig. 2.6, we can treat these distinct orientations as different species of a multi-component lattice gas, where sites are now either vacant or singly occupied by tilted molecular species U (?up? shown as red/dark gray in Fig 2.6) or D (?down? shown as blue/medium gray). This leads to a three-component (U, D, and vacancy) lattice gas that maps onto an anisotropic version of the Blume, Emery, Griffiths (BEG) or spin 1 Ising model [14]. There have been many studies of different versions of the BEG model in other contexts, though not with parameters expressing the competition between anisotropic bonding 28 (intermolecular) and substrate (single particle field) terms needed here. Fig. 2.6 Schematic illustration of molecular system with variable substrate interactions as well as intermolecular interactions. U and D correspond to up/down orientations of the molecular tilt as shown in Fig. 2.1, while G corresponds to an orientation with most favorable interaction with the substrate, possibly perfectly horizontal. The parameters ???, ???, ???, ??? correspond to lateral interactions between molecules of orientation i = U, D or G, and the parameters ??, ??, ?? correspond to interaction of molecules of orientation i = U, D or G with the substrate. As per the physical model of Fig. 2.1, we assume that strong H-bonds with energy ????? form only if a U is next to a D in the chain direction. In the side direction, the vdW interaction ????? is much less favorable, as are energies ???? and ???? of (non-H-bonded) UU or DD pairs in the chain direction, which we take 29 equal to those in the side direction for simplicity (i.e., we have set ???? = ???? = ??? and ???? = ???? = ???). Experimentally (see Fig. 2.1b), the lateral ordering of the chains provides evidence that the symmetrical side interactions (???? and ????) are somewhat more favorable than the unsymmetrical ????? . It seems quite plausible that different conformers (tilt-orientations) will have different effective interaction energies with the substrate. We arbitrarily designate U as the tilt-orientation with a more favorable interaction ??? with the substrate. Then the difference in total energy between a H-bonded UD pair and a weakly bonded UU pair in the chain direction results from a competition between a favorable intermolecular energy ??? ????? and a less favorable substrate energy ?????. To allow independent control of the density in the disordered phase, we introduce another component G to describe the dominant orientation in the disordered phase, shown as green (light gray) in Fig. 2.6. An isolated G molecule is supposed to have an optimal orientation with respect to the substrate, thus an even more favorable substrate interaction energy ??? than ??? of the favorable U orientation. The substrate interaction energy, ??, physically could include some relaxation energy of the substrate around an isolated G molecule, some of which could be lost for two adjacent G molecules. This can be represented within a lattice gas framework by a weak repulsive effective intermolecular interaction between a GG pair, ??? < 0. This contrasts with the weak attractions between UU and DD pairs needed for lateral ordering of the chains. 30 The overall Hamiltonian for this 4-component lattice gas model in its most general form is: ?? = ? ???? ?=?,?,? +? ? ???????? ?=?,??,? (2.7) Here Ni is number of molecules of species i (= U,D,G) and ?i the corresponding substrate energy, ???? is the number of nearest neighbor pairs of species i and j in the chain (? = c) or side (? = s) directions, and ???? the associated pair interaction energies. The only attractive interaction for UD pairs occurs in the chain direction,???? , corresponding to the strong H-bonds. There are many parameters to be determined, as illustrated in Fig. 2.6, but we can use the experimental observations to help choose physically relevant values. An important relationship is obtained by requiring that the T = 0 boundary energies per unit length along the straight S and CE boundaries be essentially the same. Using Eq. 2.7, it is easy to calculate the total energy at T = 0 of a large rectangular island containing a fixed number of molecules ???? = ?? ????, made up of NCE perfect laterally ordered chains of length NS molecules, shown in Fig. 2.7 below: 31 Fig. 2.7 Illustration of a rectangular island of ACA?s, composed of ?? ???? molecules. (Assume ?? is odd, because U is favored at ends of chains) ????? = (?? ?1)?????? +12(?? +1)(??? ?1)??? +12(?? ?1)(??? ?1)??? +12(?? +1)????? +12(?? ?1)????? (2.8) Both CE boundaries will optimally contain U rather than D molecules because of the more favorable substrate interaction, so NS is odd. Minimizing the energy with respect to NS with NTOT held constant, and assuming for simplicity ?UU = ?DD, gives the relation ?? ??? = ???? ?(?? ???)/2 ??? (2.9) Thus, in this multicomponent model, the anisotropy in the T = 0 island shape arising from the strong H-bond interaction ???? relative to ??? is effectively reduced by the difference in the substrate interaction (?? ???)/2. To agree with the experimental observation, the minimum energy should occur for a square island with 32 length ???? = ??????, where ?? ???? = 1/?3, and thus ?? ???? = ?3. For the case here, with the bonding asymmetry set by the ratio of hydrogen-bond to vdW interaction strength, ??? = ??? ? ???? 5? , Eq. 2.9 yields the estimate ?? ??? ? 1.3???? . In other words, to produce a symmetric domain shape consistent with the experimental observation despite the very anisotropic pair interactions requires a substrate binding energy difference comparable to the strongest intermolecular interaction. Another constraint from experiment arises from the equilibrium between the disordered gas phase, with a rather high density of isolated G molecules, and the chain phase. At low temperature, addition or removal of an H-bonded UD pair at either the S or CE boundary of a chain island creates an effective ?kink? (repeatable excitation unit). Any significant channel for interface fluctuations must allow exchange between H-bonded UD pairs at such a kink and two distinct G molecules. Efficient exchange demands a small difference in the total energies of the two configurations, so that 2?? ? 2???? +??? +??? +?? +?? (2.10) Using the estimates above in Eq. 2.10 yields ?? ??? ? 0.55???? . The room temperature desorption of submonolayer ACA films from Ag(111) on a 10-hour timescale supports a value of ?? value around 0.65 eV [3]. 2.4 Monte Carlo Simulations 33 2.4.1 Introduction to Monte Carlo simulations One of the most common methods of studying interfaces by simulations is the Metropolis Monte Carlo (MC) technique [15,16]. The stochastic MC technique generates microscopic configurations of the system, and gets the structural and thermodynamic properties. Computer simulations of a bulk or surface system are usually performed by employing periodic boundary conditions, in order to reduce surface effect due to the finite size of the system. In our 2-D surface system, the simulation cell is a 139 nm x 150 nm (240aS x 150 aCE ) rectangular lattice system with periodic boundary conditions. The simulation is in a grand canonical ensemble with (T, V, ?) fixed. Simulations were initiated with two perfectly straight boundaries, and the simulation temperature, fixed at 371K, was chosen somewhat higher than experiment to permit more rapid equilibration. 2.4.2 Simulation details Using these estimates as a starting point, we performed Monte Carlo simulations of the lattice gas model to determine parameter values consistent with the experimental observations. The simulations used a grand ensemble where changes in the total number of molecules of any species i = (G, U, D) on the substrate are controlled by a common chemical potential and the change in local energy on addition or removal of the particle as determined by the Metropolis criterion. Alternatively, a molecule of one species can directly convert to another molecular 34 species with a probability based only on the local change in energy. By testing a range of relative parameter values, we arrived at the values given in Table 2.2. The substrate binding energy differences follow the trends expected from the T = 0 estimates discussed above, with ?? ??? = 1.29???? and ?? ??? = 0.46???? . 0.65 eV ?? 0.48 eV ?? 0.37 eV ???? 0.075 eV ???, ??? 0.0 eV ??, ???, ???, ???? -0.02 eV ??? Table 2.2: Parameters determined for MC simulations. A reasonable value of 0.37 eV was chosen for the H-bond energy ???? and the ratio of ???? to ??? = ??? was set to an assumed value of 5. The symmetry balance and step stiffness are mainly determined by this ratio and the difference ?? ??? as given in Eq. 2.10. We constrained our final parameter set to satisfy Eq. 2.10 for square T=0 islands, as suggested by experiment. The density in the gas phase increases with ??, which was chosen to give a high density gas phase at experimental temperatures. Given the values of ???? , this corresponds to orientationally induced changes in the substrate interaction energy from the strongest value of 0.65 eV to only 2.5 meV. This model suggests that half the ACA molecules (the unfavorable D molecules) in the ordered structure are primarily held in place by intermolecular interactions. The chemisorption of aromatic molecules such as ACA is dominated by pi-bonding to the 35 substrate. However, nitrogen-containing heterocycles are known to adopt tilted-configurations via additional nitrogen lone-pair interactions with the substrate [17,18,19,20]. Thus a reasonable hypothesis is to assign the most favorable G state to a planar configuration, the U state to the configuration with the N atom down, and the least favorable D state to the configuration with the carboxyl group down. Fig. 2.8 shows equilibrated Monte Carlo configurations of the lattice gas model with parameters from Table 2.2 for a state with fixed temperature of 300K, where finite strips of the solid, mainly the UD chain phase, are in equilibrium with a dense vapor of mostly isolated G molecules. (a) 36 (b) Fig. 2.8 Images of Monte Carlo simulations. A condensed strip with side (S) boundaries in equilibrium with the dense gas phase is shown in the upper panel (a), and one with end (CE) boundaries is shown in the lower (b). The areas enclosed in the boxes are shown in expanded view to the right. Both simulations use a (139 nm x 150 nm) (240aS x 150 aCE ) system with periodic boundary conditions. Simulations were initiated with two perfectly straight boundaries, and the images shown are after 300,000 sweeps. The simulation temperature 371K was chosen somewhat higher than experiment to permit more rapid equilibration between the phases and the simulations determined the coexistence chemical potential as ?0.675 eV. The similarity in fluctuations at the orthogonal boundaries is evident, with root mean square (rms) widths wS = 2.23 nm and wCE = 2.20 nm, consistent with the experimental ratio close to 1. We use the Monte Carlo configurations to directly determine the step stiffness ?? at each boundary from the equilibrium spatial 37 correlation function ?(?) = ?(?(? +?0)??(?0))2?. According to the capillary wave model, in a large system of size L with periodic boundary conditions this should equal ?(?) = ???? ?(1???) (2.11) for 0 ? ? ? ?. Fig. 2.9 shows the spacial correlation function for both the side and end boundaries with the fit to Eq. 2.11. This yields accurate absolute estimates for the stiffnesses, ????= 103 meV/nm and ???= 102 meV/nm for this parameter set. Their ratio ??? ???? = 1.01? is in reasonable agreement with the experimental ratio, ??? ???? = 1.2?0.3? . The absolute value of the MC stiffness is somewhat larger than the experimental upper limit (~ 80-90 meV/nm), suggesting that the true H-bond strength may be somewhat smaller than the value of 0.37 eV assumed in the simulation. 38 Fig. 2.9 Simulation of G(y) and fits for ?? using Eq. 2.11 for the boundaries shown in Fig. 2.1. For the CE boundary, L = 150 nm, and the best fit value ??= 103 meV/nm. For S boundary, L = 139 nm, and fit value ?? = 102 meV/nm. Fig. 2.10 and Fig. 2.11 are the time correlation function G(t) and auto correlation function C(t) for both the side and end boundaries with the fit to Eq. 1.23 and Eq. 1.22. These results also verified the similarity of the two boundaries in fluctuations. The simulated G0=0.032?0.002 nm2 here also agree with the experimental results (Table 2.1). 39 Fig. 2.10 Time correlation function G(t) for simulation. The red solid line is the fitting of the Simulation data. The one above is the S boundary and G(?) = 0.034?1? , 1? = 0.509?0.002. The lower one is CE boundary and G(?) = 0.031?1?, 1? = 0.528?0.002. 40 Fig. 2.11 Auto correlation function for simulation. The red solid line is the fitting of simulation data. The one above is the S boundary and lower one is CE boundary. ?? = 59900(??????), ??? = 49500(??????). 2.5 Discussion and Conclusions Boundaries within ACA thin films fluctuate dynamically with rms displacements of 2.0 nm and underlying molecular-exchange time constants of 10-30 ms at room 41 temperature (see Table 2.1). The fluctuations are well described by the equilibrium properties of line boundaries with non-conserved thermal excitations balanced by a restoring line tension. We expect that such fluctuations will influence electron transport properties and impact electronic applications of organic thin films. For instance, reduced carrier transport across a domain boundary or the formation or recombination of electron-hole pairs at a boundary will generate a time-dependent signature due to the fluctuations in length of the boundaries, which evolve with characteristic time constants ?? ? ??????2? ? , where q is the wavenumber corresponding to the mode of fluctuation [21]. The boundaries in the ACA films displayed surprising structural characteristics ? they are isotropic in island shape, line tension and fluctuation time constant despite the strongly anisotropic intermolecular interactions, with a ratio of orthogonal H-bond to vdW interactions ???? ???~5? . We address this by adapting the BEG model to incorporate an additional energetic effect for the ACA molecule, where configurations with alternate tilt angles are observed experimentally in the condensed phase. The energy costs of the tilted configurations are introduced as single-particle field effects, and are modeled as discrete, but interchangeable, molecular components, parameterized as ?? and?? . Physical analysis of the low-temperature limit of the model shows that the effective anisotropy of the overlayer can be changed dramatically, as represented by the shape ratio: ??? ?? = ???? ?(?????)/2? ?? . Monte Carlo simulation using the model yields quantitative agreement with the experimental observations, and suggests that molecular tilt can reduce the molecular-substrate interaction strength from its 42 optimum value near 0.65 eV to near zero. The effect of molecular tilt is frequently observed in N-heteroaromatics [16-19], and thus the symmetry modifications observed for ACA should also be common. By considering the shape ratio shown above, it is easy to predict the consequences of modifying the intermolecular interactions. In the ACA case, where the asymmetry of the molecular interactions ???? ???? is large, the tilt effect effectively symmetrized the boundaries, thus creating compact domains. Conversely, modifying the molecular structure to reduce the interaction asymmetry (e.g. reduce ???? ???? ) by increasing the molecular interactions orthogonal to the hydrogen-bond direction would exaggerate the asymmetric shape ratio for the same relative tilt-effect. Beyond this one example, many other combinations can be evaluated using this model, effectively yielding a molecular basis for predicting thin film morphology. The ability to tune domain shape will allow the systematic design of organic thin film systems that will self-assemble boundary configurations favorable for charge transport and recombination [22,23]. References: 1 B. Xu, C.G. Tao, W.G. Cullen, J.E. Reutt-Robey and E.D. Williams, ?Chiral Symmetry Breaking in Two-dimensional C60-ACA Intermixed Systems?, Nano Lett. 5, 2207-2211 (2005) 2 B. Xu, C.G. Tao, E.D. Williams and J.E. Reutt-Robey ?Coverage Dependent Supramolecular Structures: C60:ACA monolayers on Ag(111)?, J. Am. Chem. Soc. 128, 8493-8499 (2006) 43 3 B. Xu, B. Varughese, D. Evans and J. Reutt-Robey ?Morphology Selected Molecular Architecture: Acridine Carboxylic Acid Monolayers on Ag(111)?. J. Phys. Chem. B 110, 1271-1276 (2006) 4 G.A. Jeffrey, An Introduction to Hydrogen Bonding (Oxford University Press, 1997) 5 S. Tsuzuki, K. Honda, T. Uchimaru, M. Mikami and K. Tanabe, ?Origin of Attraction and Directionality of the Pi/Pi Interaction: Model Chemistry Calculations of Benzene Dimer Interaction?, J. Am. Chem. Soc. 124, 104, (2002) 6 M. Wortis, in Chemistry and Physics of Solid Surfaces VII, eds. R. Vanselow and R. F. Howe, 367-428 (Springer Verlag, Berlin, 1988) 7 E.D. Williams, R.J. Phaneuf, J. Wei, N.C. Bartelt and T.L. Einstein, ?Thermodynamics and Statistical Mechanics of the Faceting of Stepped Si(111)?, Surface Science 294, 219-242 (1993) 8 M. Giesen-Seibert and H. Ibach, ?On the Time structure of Tunneling Images of Steps?, Surface Science 316, 205-222 (1994) 9 H.-C. Jeong and E.D. Williams, ? Steps on Surfaces: Experiment and Theory?, Surface Science Reports 34, 171-294 (1999) 10 E. Le Goff, L. Barbier, and B. Salanon, ?Time?space Height Correlations of Thermally Fluctuating 2-d Systems: Application to Vicinal Surfaces and Analysis of STM Images of Cu(115)?, Surface Science 531, 337-358 (2003) 11 N.C. Bartelt, J. L. Goldberg, T.L. Einstein, E.D. Williams, J. C. Heyraud and J.J. Metois, ?Brownian Motion of Steps on Si(111)?, Phys. Rev. B 48 15453- 15456 (1993) 12 O. Bondarchuk, D.B. Dougherty, T.L. Einstein and E.D. Williams, ?Dynamics of Step Fluctuation on a Chemically Heterogeneous Surface of Al/Si(111)?, Phys. Rev. B 66, 085327 (2002) 44 13 A. Pimpinelli, J. Villain, D.E. Wolf, J.J. M?tois, J.C. Heyraud, I. Elkinani and G. Uimin, ?Equilibrium Step Dynamics on vicinal Surfaces?, Surface Science 295, 143-153 (1993) 14 M. Blume and V. J. Emery? ? Ising Model for the ? Transition and Phase Separation in He3-He4 Mixtures?? Phys. Rev. A 4, 1071 (1971) 15 N. Metropolis and S. Ulam, ?The Monte Carlo Method?, Journal of the American Statistical Association, 44, 335 (1949) 16 N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A.H. Teller and E. Teller, ?Equation of State Calculations by Fast Computing Machines?, J. Chem. Phys. 21, 1087 (1953) 17 M. Bader, J. Haase, K.-H. Frank and A. Puschmann, ?Orientational Phase Transition in the System Pyridine/Ag(111): A Near-edge X-ray-Absorption Fine-Structure Study?, Phys. Rev. Lett. 56, 1921 (1986) 18 A.L. Johnson and E.L. Muetterties ?Chemisorption Geometry of Pyridine on Pt(111) by NEXAFS?, J. of Phys. Chem. 89, 4071-4075 (1985) 19 J. Bonello and R.M. Lambert ?The Structure and Reactivity of Quinoline Overlayers and the Adsorption Geometry of Lepidine on Pt(111): Model Molecules for Chiral Modifiers in Enantioselective Hydrogenation?. Surface Science 498, 212-228 (2002) 20 J. Kubota and F. Zaera ?Adsorption Geometry of Modifiers as Key in Imparting Chirality to Platinum Catalysts?, J. Am. Chem. Soc. 123, 11115-11116 (2001) 21 E.D. Williams, O. Bondarchuk, C.G. Tao, W. Yan, W.G. Cullen, P.J. Rous and T. Bole, ?Temporal Step Fluctuations on A Conductor Surface: Electromigration Force, Surface Resistivity and Low-frequency Noise?, New Journal of Physics 9, 387 (2007) 45 22 B.A. Gregg, ?The Photoconversion Mechanism of Excitonic Solar Cells?, Materials Research Society Bulletin 30, 20-22 (2005) 23 R. Otero, D. Ecija, G. Fernandez, J. M. Gallego, L. Sanchez, N. Martin, and R. Miranda, ?An Organic Donor/Acceptor Lateral Superlattice at the Nanoscale?, Nano Lett. 7, 2602-2607 (2007) 46 Chapter 3 Pair Approximation: Theory for Competition Between Substrate and intermolecular Interactions 3.1 Introduction As illustrated by the system discussed in last chapter, many self-assembled monolayers (SAMs) on surface systems have been studied in experiments, and there are a variety of applications [1,2,3]. The structures adopted by these systems could depend on a number of things: the geometric shape of unit molecules, details of the intermolecular and substrate interactions (mostly weak vdW interactions, H-bonds and strong electronic interactions making up chemical bonds) and temperature. Compared to single atom units, the complexities of the geometry of molecules and interactions may result in the complicated pattern structure including anisotropic islands and sub-lattices formation, which give us difficulties in developing theoretical models. On the other hand, the large sizes of the molecule also give us the possibility to use only nearest neighbor interactions to build the model so that we may simplify our theory, because the vdW interaction is so short ranged compared to the molecule size. Therefore, the nearest neighbor lattice gas model could be a good choice to characterize many of these systems. Here we will focus on the competition between the intermolecular and substrate interactions as illustrated in the ACA/Ag(111) system from last chapter [4]. There we used a four-component anisotropic lattice-gas model, with many parameters to try to 47 make direct contact with experiment. But the basic idea of a competition between nearest-neighbor intermolecular interactions and substrate interactions seem very general, and likely to be seen in other systems. Thus, rather than focusing on details of a particular system, in this chapter we will study even simpler lattice gas models that can exhibit competition between intermolecular pair interactions and substrate interactions for each member of the pair. To simplify the theoretical calculations still further, we can also use mean field theory or the pair approximation to study general qualitative features of the model. The pair approximation is a generalization of mean field theory which still maintains much of the simplicity of a mean field treatment while often giving more accurate answers for many (non-critical) properties. We will say much more about the pair approximation later in this chapter. We use the ACA/Ag(111) system to motivate our simplified lattice gas model. The ACA molecule, with a roughly rectangular structure, has anisotropic interactions along the short and long axis. But anisotropy need not be an essential feature of the basic competition between substrate and intermolecular interactions we wish to focus on, so we will consider both anisotropic and isotropic versions of our simplified model. The STM study also revealed an ?up? ?down? sub-lattice structure along the short axis. Large scale experiments showed us a high density gas phase and a similar E-boundary and S-Boundary fluctuations [2]. In analogy to the ?up? and ?down? sub-lattice structures of the ACA monolayer, we imagine two different components in the lattice gas, here corresponding to two 48 different orientations of a single molecular species on the surface. It is convenient to think of a vacant lattice site as another component, so we arrive at a three component lattice gas system, where each site on a square lattice is either a ?vacancy?(V) or occupied by an ?up?(U) or ?down?(D) molecule. Generally, we think that U and D have different substrate interactions ??? and ???. The U molecule has a more favorable interaction ??? with the substrate compared to the D molecule which has substrate interaction ???. The substrate interaction for a vacancy is 0 because there is no physical interactions between V and substrate. In this chapter, to simplify our model, we will not include the G molecules so we have only a three component lattice gas. The intermolecular interactions in the ACA experiments, are anisotropic: along the short axis, chain direction, there can be strong U-D hydrogen bond with ???? = ????? . But along the long axis, side alignment, two adjacent U or D molecules have weak side-side (U-U or D-D) vdW interaction ??? = ???,?? = ???,??. The U-D pairs in side alignment and U-U, D-D pairs in chain directions are very rare to see. If the pair interaction ????? is strong enough, we could expect the lowest energy pattern on the surface should be chains ?U -D-U-D?. But forming a U-D pair causes an unfavorable D substrate interaction as well. If ????? is not that strong, relative to the difference of the substrate interaction ??? and ???, the substrate could dominate and the whole surface could be covered with substrate favored ?U? 49 molecules. So, the competition of the substrate interaction and intermolecular interactions will determine the monolayer patterns that form. As shown in the last chapter, this competition can change the ACA island shape and control boundary fluctuations [2]. The overall Hamiltonian for our simplified anisotropic 3-component lattice gas is: ?? = ? ???? ?=?,?,? + ? ? ???? ???? ?=?,??,?=?,? (3.1) Here N? is number of molecules of species ?(= U,D) and ?? the corresponding substrate energy, ???? is the number of nearest neighbor pairs of species ? and ? in the chain (i=c) or side (i=s) directions, and ???? the associated pair interaction energies. The parameter values can be chosen relevant to the experiments. In this paper, we will focus on the global phase diagrams corresponding to a wide range of parameter values. Finally, we will consider an even simpler lattice gas model that exhibits competition between intermolecular and substrate interactions in its simplest possible form. This isotropic model has the same U-D interactions ????? in both directions and no U-U and D-D interactions. It will be interesting to see how the phase diagram will change in going from the anisotropic to the isotropic model. The Hamiltonian for the isotropic lattice gas model can be written as: 50 ?? = ? ???? ?=?,?,? + ? ?????? ?,?=?,?,? (3.2) The isotropic three-component (U, D, and V) lattice gas maps onto a version of the Blume, Emery, Griffiths model or spin 1 Ising model [5]. If ?? = ??, it will change into the regular antiferromagnetic model. If ?? ? ??, it will change into the single component, U, regular lattice gas or spin 1 2? Ising model. 3.2 Pair approximation One general method to treat the Ising model is mean field theory, which assumes that there is no short ranged order apart from that which follows from long-range order, and it ignores the possibility of local correlation between species. Here instead of considering only one molecule, we will consider the distribution function of two adjacent molecules and the local pair interactions will come in. This method is called pair approximation [6,7,8] (also called Bethe-Peierls-Weiss or the quasi-chemical approximation), which is qualitatively more accurate than the mean field theory and corrects an inadequacy in the mean field theory in the anisotropic limit [9]. The standard derivation of the pair approximation uses xi, the probability of a lattice point having the respective species, and yi, the probability of appearance of a pair having one of the configurations, as the basic parameters. The total energy E as well as the entropy S of the system can be expressed in terms of xi and yi using approximate and complicated combinatorial arguments. By minimizing the grand potential ?? = ? ??? ???, we can get a set of equations with unkowns xi and yi. 51 Solving these equations will finally give us information about the system [10,11,12,13]. Here, we use a simpler and more intuitive method to derive the pair approximation equations [9]. Instead of enumerating the occupation possibilities, we will define effective interactions for each species, and use these effective interactions as the basic parameters to express the singlet and pair distribution functions, the energy, entropy and grand potential. This method will establish close connections between mean field theory and the pair approximation and show why the latter can be more accurate in many applications. Compared to other methods (Kikuchi method and Bethe-Peierls method [8]), this derivation is physically more direct, clearer and easier to understand. We even do not need the entropy to get the final equations about the neighbor effective interactions, and they are more convenient to solve by iteration methods. 3.2.1 Pair approximation for two species First let?s think about the simpler two species situation: ?A? and ?V?, where ?A? means atom and ?V? means vacancy. In spin language this corresponds to the usual spin 1 2? Ising model discussed in many textbooks. Here we use ?? and ?? to denote the effective interactions around A and V, which means, the ratio of the probabilities to find the specie A and V at a particular site is (Fig. 3.1): 52 A A ? A A A 1 1 ? V 1 1 Fig. 3.1 Effective interactions for two species 2-D square lattice gas. The circled species in the middle are input species we are thinking about, the surrounding spaces are the effective interactions. ?(?) ?(?) = ?4???+? ?4??? (3.3) where ? is the chemical potential of component A. We assume ?? = 0, because there is no interaction between V and other species, and also we only need one parameter ?? to satisfy above equation. This can be seen as the definition of ??, and we can apply this representation to various calculation methods, including mean field theory (MFT) and pair approximations. For example, the MFT assumes that the role of the neighboring particles is to form an average molecular field which acts on the center particle. Combining this assumption with our definition of effective interaction, we have the relation ?? = ? ??, where n is the average density ? ? P(A), and ? is the AA pair interaction energy. By taking this relation into Eq. 3.3, and taking count of P(A)+P(V)=1 will give us the fulfilled equation of MFT. For pair approximation, instead of only one particular site, we focus on two adjacent sites and make another assumption that, even if we put two adjacent atoms 53 AA or AV or VV together, their distribution possibilities can still be expressed with their surrounding effective interactions (six in case here) the same as the situation for only one atom at the center (Fig. 3.2). ?? ?? ?? ? A ? A ?? ?? ?? ?? 1 ?? ? A ? V 1 ?? 1 Fig. 3.2 Pair approximation for two species 2-D square lattice gas. The circled species in the middle are input species we are thinking about, the surrounding spaces are the effective interactions. And the assumed relation for approximation: ?(??) ?(??) = ?6???+2? ?3???+? (3.4) For the two species situation, we can easily map it to the spin representation, where "?" ? "+", "?" ? "?" (Fig. 3.3): ?+??? ?+??? ? + ?+??? ?+??? ????? ????? ? - ????? ????? Fig. 3.3 Spin representation of two species 2-D square lattice gas. The circled spin in the center are spins we try to put inside, and the surroundings are effective environment spins interact with center spins. 54 ?(+) ?(?) = ?4???+??????? ?4????????+?? = ? 4??(2??)?2?? (3.5) where ?? = (?+???+?????)/2. What seems tricky is that we can?t solve for ?+??? or ????? separately, but we can solve for ??, both in MFT and pair approximations [9]. And the mapping relations between the atom representation and spin representation are: ? = ?2H?8J ? = 4J ?? = 2???+2? (3.6) The Hamiltonian can be written as: ? = ?4???+???+? (for +) ? = 4????????? (for ?) (3.7) Or ? = 4????+,??????+?? (3.8) 3.2.2 Pair approximation for three species Now we come to our three species model. For simplicity, we first consider the isotropic situation. Assume we have U and D, two kinds of molecules, and they prefer to stay next to each other, as in an antiferromagnetic model, they may form an interlaced pattern in solid phase on the surface. The possibilities for different species on two adjacent sites might be different. So there might be a sub-lattice composed of ?white? and ?shaded? lattices (Fig. 3.4): 55 U D U D U D U D U D U D U D U D U D U D Fig. 3.4 Sub-lattices formed by the U and D molecules. When the U and D formed an interlaced structure, unlike the mean field model, the two neighbor site would have different properties and probabilities for species put on it. First let us think about a ?white? site. On this site, the probabilities of putting ?(U, D or V) depends on the property of ? itself (substrate interaction ?? and chemical potential ?) and the surroundings (most possible ?D? species), which we assume give an effective interaction ?? ?(?) = ? 4???+??+??? ? ?(?) = ? 4???+??+??? ? ?(?) = 1? (3.9) More generally, let P(?) denote the probability to find an ?(U, D or V) specie at the white site. Each of the ? (here ?=4) neighbors of this site contributes an effective factor exp (???) to the probability P(?). More explicitly, ?? is a number that makes the probability of having an ? at a white site, ?(?) = exp(4??? +???? +??)/?, where ?=U, D, and V. Z is the normalization factor ? = ? exp(4??? +????)? , ?? = 1 for ?=U, D, and ?? = 0 for ? =V. This number ?? always exists for a given uniform phase of a system. 56 Similarly, let ?(??) denote the probability to find an ?(U, D or V) at the shaded site and each of the ? (here ?=4) neighbors contribute the factor exp (????) to the probability of putting ? at this site. Notice that ??? may equal to ?? (in which case we do not need of sub-lattice at all). Since the vacancy has no interactions with other species, the contribution factors exp(???) = exp(????) = exp(? ?0) = 1. The possibilities to find U, D and V at a particular site, expressed with all those effective interactions, are shown in Table 3.1 below: White Site ? = ?4???+??+??? +?4???+??+??? +1 ?? ?? ? U ?? ?? ?(?) = ? 4???+??+??? ? ?? ?? ? D ?? ?? ?(?) = ? 4???+??+??? ? 1 1 ? V 1 1 ?(?) = 1? 57 Shaded Site ?? = ?4????+??+??? +?4????+??+??? +1 ??? ??? ? U ??? ??? ?(??) = ? 4????+??+??? ?? ??? ??? ? D ??? ??? ?(??) = ? 4????+??+??? ?? 1 1 ? V 1 1 ?(??) = 1?? Table 3.1 Schematic representation of the effective contribution factors for three species isotropic 2-D square lattice gas. The species circled in the center are species we try to put in, while the surroundings are effective interactions. White sites and shaded sites may have different preference. Z is the normalization factor. Now we make the pair approximations: For two species ? and ? on a pair of adjacent sites, the other 6 neighbors still contribute the same factor ??? (???) and ??? (????) to the possibilities of ? and ?. The possibilities of appearing two specific species configurations are shown below in Table 3.2. 58 W = ??(3??+3???+???+2??+2?) +??(3??+3???+???+??+??+2?) +??(3??+??+?) +??(3??+3???+???+??+??+2?) +??(3??+3???+???+2??+2?) +??(3??+??+?) +??(3???+??+?) +??(3???+??+?) +1 ?? ??? ?? ? U ? U ??? ?? ??? ?(???) = ? ?(3??+3???+???+2??+2?) ? ?? ??? ?? ? U ? D ??? ?? ??? ?(???) = ? ?(3??+3???+???+??+??+2?) ? ?? 1 ?? ? U ? V 1 ?? 1 ?(???) = ? ?(3??+??+?) ? ?? ??? ?? ? D ? U ??? ?? ??? ?(???) = ? ?(3??+3???+???+??+??+2?) ? ?? ??? ?? ? D ? D ??? ?? ??? ?(???) = ? ?(3??+3???+???+2??+2?) ? ?? 1 ?? ? D ? V 1 ?? 1 ?(???) = ? ?(3??+??+?) ? 59 1 ??? 1 ? V ? U ??? 1 ??? ?(???) = ? ?(3???+??+?) ? 1 ?? 1 ? V ? U ?? 1 ?? ?(???) = ? ?(3???+??+?) ? 1 1 1 ? V ? V 1 1 1 ?(???) = 1? Table 3.2 Schematic representation of the pair approximation method for a three species isotropic 2D square lattice gas. W is the normalization factor. Now we can determine the parameters ??, ??, ???, and ??? from the consistency conditions below: ?(?) ?(?) = ?(???)+?(???)+?(???) ?(???)+?(???)+?(???) ?(?) ?(?) = ?(???)+?(???)+?(???) ?(???)+?(???)+?(???) ?(??) ?(??) = ?(???)+?(???)+?(???) ?(???)+?(???)+?(???) ?(??) ?(??) = ?(???)+?(???)+?(???) ?(???)+?(???)+?(???) (3.10) Taking into the probabilities in Table 3.1 and Table 3.2, the equations should be: ?4???+?? 1 = ?3???+3????+????+2??+??? +?3???+3????+????+2??+??? +?3???+?? ?3????+??+??? +?3????+??+??? +1 60 ?4???+?? 1 = ?3???+3????+????+2??+??? +?3???+3????+????+2??+??? +?3???+?? ?3????+??+??? +?3????+??+??? +1 ?4????+?? 1 = ?3????+3???+????+2??+??? +?3????+3???+????+2??+??? +?3????+?? ?3???+??+??? +?3???+??+??? +1 ?4????+?? 1 = ?3????+3???+????+2??+??? +?3????+3???+????+2??+??? +?3????+?? ?3???+??+??? +?3???+??+??? +1 3.11 Let ? = ????,? = ????,?? = ?????,?? = ?????, the equations can be simplified as: ? = ? ?3??(???+??+?) +??3??(???+??+?) +1 ??3??(??+?) +??3??(??+?) +1 ? = ? ?3??(???+??+?) +??3??(???+??+?) +1 ??3??(??+?) +??3??(??+?) +1 ?? = ? 3??(???+??+?) +?3??(???+??+?) +1 ?3??(??+?) +?3??(??+?) +1 ?? = ? 3??(???+??+?) +?3??(???+??+?) +1 ?3??(??+?) +?3??(??+?) +1 (3.12) We can just simply numerically solve these equations by iterative technique. The final convergence results depend on the appropriate initial values of u, d, u? and d?. And at some particular temperature T and chemical potential ?, we can get two set of solutions based on different initial values. Then the real physical solution should be the one with lower grand potential: ?? = ? ??? ??? (3.13) To get the grand potential ?G, we need to use the Kikuchi method to calculate the system entropy S. 61 ? = ??? [ 3 2 ? ?(?)ln?(?)?=?,?,?, ??,??,?? ?2 ? ?(?,?)ln?(?,?) ?=?,?,? ?=??,??,?? ] (3.14) The two equations Eq. 3.13 and Eq. 3.14 above are also the starting point of the standard derivation of the pair approximation method. You need to differentiate the grand potential with respect to the independent variables P(?,?), and put the derivatives to zero to get the fundamental equations [11,12,13]. These equations will be equivalent to Eq. 3.12 we get earlier. In the method we use in this thesis, we calculate the entropy and grand potential after we derive and solve the equations (Eq. 3.12), to get the phase diagrams. 3.3 Calculation Results After we solved these equations, we can compute the values of P(?) and P(?,?) and then the energy, the grand potential and the coverage density. By fixing the temperature, we can get the ?G-? graph as shown in Fig. 3.5 below. We can see that, below the critical temperature T0, there is a region where there are two calculated ?G?s correspond to one ?, and the lower ?G is the physical one, the other is unstable. The crossing point at ?=?0 is the phase coexistence point. Below this point is gas phase; above this point is solid phase. 62 Fig. 3.5 Grand potentials ?G vary with chemical potential ?. ?UD=0.18eV, EU-Ed=0.35eV, T0~1465K If we always fix the ? at the phase coexistence ?0(T) at different temperatures T, and calculate the phase coexistent densities, we can get the phase diagram ??-T. In Fig. 3.6 we can see that, comparing with the regular Ising model, adding a third species can distinctly change the shape of the phase diagram. And you don?t need to be very near to the critical point to reach a higher gas density. 63 Fig. 3.6 Phase diagrams with different parameters. Ed=0.0eV If we divided the temperature by the critical temperature T0, then we get the scaled phase diagram Fig. 3.7. We can see more clearly of this change. We also noticed that this change depends on the ratio of the two competition parameters, (?U ??D) ?UD? . To reach the same gas density ?, the larger of the ratio (?U ??D) ?UD? , the lower the needed normalized temperature T T0? . When (?U ??D) ?UD? get large enough, say ?U ? ?D, then there is very low possibility of D molecules and the system goes to the limit of regular Ising model of U molecules with interactions ?UU. 64 Fig. 3.7 Phase diagram scaled with critical temperature T0 3.4 Conclusions and discussions We have presented a simple method for the calculation of molecule SAMs. This method uses three species lattice gas model and pair approximation. It?s easy to get the equations and solve numerically. There are some other systems could also be modeled in a similar way. For example: the up/down square patterns formed in the naphthalene-tetracarboxylicacid-dianhydride (NTCDA)/Ag(111) system [14], the binary phthalocyanine superstructure system [15] 65 3.4.1 Anisotropic situations Compared with the system we studied in Chapter two, we still need to consider the anisotropic situation. And we only need to make few modifications to apply it to this more complicated case. In anisotropic situation, ???? ? ???? ,???? ? ???? , which happened in the experiment example in Chapter 2, is shown in Fig 3.8 below: U U U U U D D D D D U U U D U D D D D D Fig. 3.8 Anisotropic 2D square lattice gas Instead of ??, ??, ???, ???, four parameters, we need ???, ???, ???? , ???? and ???, ???, ????, ???? eight parameters to characterize all possible interactions. The probability for a given site is in Table 3.3 below. White Site ? = ?2????+2????+??+??? +?2????+2????+??+??? +1 ??? ??? ? U ??? ??? ?(?) = ? 2????+2????+??+??? ? 66 ??? ??? ? D ??? ??? ?(?) = ? 2????+2????+??+??? ? 1 1 ? V 1 1 ?(?) = 1? Shaded Site ?? = ?2????? +2?????+??+??? +?2????? +2?????+??+??? +1 ???? ???? ? U ???? ???? ?(??) = ? 2????? +2?????+??+??? ?? ???? ???? ? D ???? ???? ?(??) = ? 2????? +2?????+??+??? ?? 1 1 ? V 1 1 ?(??) = 1?? Table 3.3 Schematic representation of the effective contribution factors for three species anisotropic 2-D square lattice gas. The species circled in the center are species we focus on, while the surroundings are effective interactions. Z and Z? are the normalization factor. 67 Correspondingly we can get the effective interaction table for pair approximations similar to Table 3.2, but more complicated since we have two adjacent sites both vertically and horizontally. Then we can get the self-consistency equations which composed of eight equations to solve these eight parameters. 3.4.2 Spin representation This 3 component lattice-gas model can also map on to a spin-1 Ising model: U molecule has spin +1, D molecule has spin -1, V molecule has spin 0. The Hamiltonian would be: ?? = ? ???? ?=?1,0,?1 + ? ??1?2??1?2 ?1,?2=+1,?1 ?? = ?? ???2 ? +?? +??2 ?2 (?+1 = ??,?0 = 0,??1 = ??) ??1?2 = ??? ????2 ?1 ??2 +??? +???2 (?12?22) (3.15) And the resulting Hamiltonian: ?? = ? (??? ????2 ?1?2 +??? +???2 ?12?22) ?1,?2 +?(?? +??2 ?2 +?? ???2 ? ? ) = ?(? ??1?2 +? ??12?22) ?1,?2 +?(????2 +? ?? ? ) (3.16) where ? = ??? ????2 ? = ??? +???2 68 ?= ??? +??2 ? = ?? ???2 (3.17) Compared to BEG model, there are two differences: one is that there is no constant term here because when all spins are zero, it means there are nothing on the surface, so Hamiltonian should be zero. Where in BEG model, S=0 means He-3, and there would be corresponding Hamiltonians for single He-3 and interactions between He?s. The second difference is that here we have a ? H?SS term corresponding to the different substrate interactions for U and D molecules. Where in BEG model there is not this term because there is no difference between S = 1 and S = -1 for a single He-4 (difference comes from the He-4 He-4 interactions). So basically speaking, our system is different from the BEG model. But when H = 0, which means U and D have same substrate interactions, our model would be one particular situation of BEG model, with ?3 = zK33, where ?3 is the chemical potential of He-3 and ?K33 is He-3 He-3 interaction energy. In this situation, if system is full with He-3, the total energy is zero, this is the same as our case: S=0 means vacancy. Under the H = 0 condition, we define ? ? ?? = ???+???? ?????? . R is negative if U-D pair is favorable than U-U pair, R is positive otherwise. The R positive situation has been thoroughly discussed in BEG?s paper. And the K negative situation has been discussed by A. Z. Akheyan and N. S. Ananikian [16]. As for J negative situation, normally we can just use this sub-lattice method and flip one sub-lattice and change 69 the sign of J. But here in this paper the two sub lattices have different substrate interactions so we have a ? H?SS term, so we can not simply flip the spin. This model is a small region in the phase diagram of BEG model that has not been discussed before. And we get the phase diagram by using pair approximations The case, H ? 0, has been examined by R.R. Levitskii, etc [17,18], using two particle cluster approximation, which will yield the same result with pair approximations. But only for situations when ?? > 0 and ?? > 0, which includes ferromagnetic, paramagnetic and quadrupolar phases. 3.4.3 Another method with fixed surroundings I can also use another method to describe the system. In this way, I can use a fixed ??, no matter if we put A or V at the center. ?? ?? ? A ?? ?? ?? ?? ? V ?? ?? Fig. 3.9 Fixed surrounding field illustration. The circled species in the middle are input species we focus on, the surrounding spaces are the effective field. But this ?? is a mixture of A and V, ?? = ?? +?? so that: ?(?) ?(?) = ?4??(?,??) ?4??(?,??) = ?4???? 1 (3.18) This method can also use the MFT or pair approximations. And it can also map 70 to the spin representation. But instead of a ?mixture?, I can just use a fixed number ??, to represent it. ?? ?? ? + ?? ?? ?? ?? ? - ?? ?? Fig. 3.10 Fixed surrounding spins. The circled species in the middle are input spins we focus on, the surrounding spaces are the average spin. The advantage to do this is that I can write the Hamiltonian simply as ? = ?4?????? ??? (3.19) No matter S is ?+? or ?-?. And this method can also be applied with MFT and pair approximations. Finally we found that the ?? is the same as the ?? in our previous method! Because the P(+)/P(-) has the exactly the same expressions with ?? both these two methods. But it?s not always like this. In the three species system, it will be totally different. The corresponding spin representation Hamiltonian with three species model, ? = ?J????? ?? ?????2??2 ?? ????? ? +????2 ? (3.20) in previous method we studied before, is to find ?+???, ?????, and ????? (= 0), so that 71 ?(?) = ? 4?(??+????+??+????2)+????? ? ?(?) = ? 4?(????????+???????2)?????? ? ?(?) = 1? (3.21) Now in this method, I need to find two particular numbers ?? and ?2??? (? (??)2) so that ?(?) = ? 4?(????+???2????)+???????2????? ? (3.22) No matter what S is. In summary, both these methods, they can be used in particle and spin representation; they both can be used with different approximations, MFT and pair approximation. The method we used before is very good at pair approximations (equations are simple and easy to be numerically solved by iteration), because pair approximation only needs interactions. No matter what effective parameters you use, finally you will have to calculate ??, the effective interactions. The new method now is often used in spin language and MFT like in the BEG paper, since the form of Hamiltonian will be simple to get, just plug in the ?? into the 72 ? = ?J????? ?? ????? ? (3.23) References: 1 J. Weckesser, J. V. Barth, C. Cai, B. Muller, and K. Kern, ?Binding and Ordering of Large Organic Molecules on and Anisotropic Metal Surfaces: PVBA on Pd (110)?, Surface Science 431, 168-173 (1999) 2 K. Braun and S. Hla, ?Probing the Conformation of Physisorbed Molecules at the Atomic Scale Using STM Manipulation?, Nano Lett. 5, 73-76 (2005) 3 R. Hiesgen, M. Rabisch, H. Bottcher, and D. Meissner, ?STM Investigation of the Growth Structure of Cu-phthalocyanine Films with Submolecular Resolution?, Solar Energy Materials & Solar Cells 61, 73-85 (2000) 4 C. Tao, Q. Liu, B. C. Riddick, W. G. Cullen, J. Reutt-Robey, J. D. Weeks, and E. D. Williams, ?Dynamic interfaces in an organic thin film?, Proc. Natl. Acad. Sci. USA, 105, 16418-16425 (2008) 5 M. Blume and V. J. Emery, ?Ising Model for the ? Transition and Phase Seperation in He3-He4 Mixtures?, Phys. Rev. A 4, 1071 (1971) 6 H. A. Bethe, ?Statistical Theory of Superlattices?, Proc. R. Soc. Lond. A 150, 552-575 (1935) 7 S. Katsura and m. Takizawa, ?Bethe Lattice and the Bethe Approximation?, Progress of Theoeretical physics, 51, No. 1 (1974) 8 K. Huang, Statistical Mechanics, 2nd edition, 342-367 (Wiley, 1987) 9 J. D. Weeks and G. H. Gilmer, ?Pair Approximation Equations for Interfaces and Free Surfaces in the Ising Model?, J. Chem. Phys. 63, 3136 (1975) 10 R. Salazar and L. D. Gelb, ?Application of the Bethe-Peierls Approximation to A Lattice-gas Model of Adsorption on Mesoporous Materials?, Phys. Rev. E 71, 041502 (2005) 73 11 R. Kikuchi, ?Theory of Domain Walls in Ordered Structures II?, J. Phys. Chem. Solids, 23, 137-151 (1962) 12 J.-C. Lee, Thermal Physics - Entropy and Free Energies, chapter 5, (New Jersey: World Scientific, 2002) 13 R. Kikuchi, ?A Theory of Cooperative Phenomena?, Phys. Rev. 81, 988-1003 (1951) 14 C. Stadler, S. Hansen, I. Kroger, C. Kumpf and E. umbach, ?Tuning Intermolecular Interaction in Long-range-ordered Submonolayser Organic Films?, Nature Physics, 5, 153 (2009) 15 A. Calzolari, W. Jin, J. E. Reutt-Robey, and M. B. Nardelli, ?Substrate-mediated Intermolecular Hybridization in Binary Phthalocyanine Superstructures?, J. Phys. Chem. C, 114, 1041 (2010) 16 A. Z. Akheyan and N. S. Ananikian, ?Global Bethe Lattice Consideration of the Spin-1 Ising Model?, J. Phys. A: Math. Gen. 29, 721-731 (1996) 17 R. R. Levitskii, O. R. Baran, and B. M. Lisnii, ?Phase Diagrams of Spin-1 Ising Model with Bilinear and Quadrupolar Interactions under Magnetic Field, Two-particle Cluster Approximation?, Eur. Phys. J. B 50, 439 (2006). 18 C. Temirci, A. Kokce, and M. Keskin, ?Equilibrium Properties of a spin-1 Ising System With Bilinear, Biquadratic and Odd Interactions?, Physica A 231, 673-686 (1996) 74 Chapter 4 C60/ZnPc 4.1 Introduction Zinc Phthalocyanine (ZnPc) and C60 are well known electron donors and accepters that can be used to construct organic photovoltaic (OPV) cells [1]. The particular arrangement of these molecules on organic and metal surfaces is clearly important in controlling charge transport and other properties of such cells. As discussed below, experiments show that the C60 molecules form meandering chains of one molecular width. This is very different from what happens on pure metal surfaces like Au, Cu or Ag [2,3,4,] where C60s form close packed islands. Here we visualize and simulate the patterns that C60 molecules form on zinc phthalocyanine (ZnPc) surfaces. We reproduce qualitative features of these patterns by modifying the standard Girifalco [5] potential between C60 molecules, which describes van der Waals interactions that are present even in the gas phase, to account for additional electrostatic components that may arise from charge transfer to the C60 molecules from the substrate. Finally we describe other aspects of the corrugated substrate on the C60 through a spatially-dependent static external field. C60 monolayer deposition on metal surfaces like Ag, Cu, and Au has been carefully studied both experimentally [2] and theoretically [3,4]. On metal surfaces, C60 molecules normally form close packed structures, and there is a charge transfer of 0.4-0.8 electron from the substrate to the C60 molecules. But the C60 monolayer on an organic surface is more complicated and not well understood. 75 There are potentially important applications of these organic bilayer systems. Organic materials like C60, ZnPc, and pentacene (Fig. 4.1) can be used to make photovoltaic devices [6]. The C60/ZnPc solution and deposition system has been studied by photoluminescence excitation and emission spectra methods [7] and co-deposition methods [8], and both indicate charge transfer from ZnPc to C60. Studying the details of the bi-layer interface can help us to explore the electronic properties of the materials and build more efficient devices. The motion of molecules on the surface is controlled by the substrate potentials and the intermolecular interactions. One example is the C60/ACA/Ag(111) system (Fig. 4.2), which has shown very interesting phenomena [9, 10]. C60 ZnPc Pentacene Fig. 4.1 C60, ZnPc, and Pentacene structures 76 Fig. 4.2 C60/ZnPc/Ag(111) System illustration. The C60 molecules are deposited on one monolayer of ZnPc molecules on Ag(1 1 1) surface. The ZnPc molecule, shown in Fig. 4.3, has a phase with D4h symmetry. Previous STM studies of ZnPc adsorption have shown that they can form a nearly square grid structure on Ag(111) surfaces. (a) (b) Ag Substrate ZnPc C60 293 01.037.1 02.044.1 ?? ?? ?? ? nmb nma 77 Fig. 4.3 (By W. Jin, et al.) ZnPc molecules on Ag (111) surface. (a): STM image of one monolayer of ZnPc molecules on Ag(111) surface. (b): the geometry structure of ZnPc monolayers. Upon deposition on the ZnPc layer, C60 aggregates into sub-monolayer islands, leaving large regions of exposed ZnPc. Island formation, and the lack of structural change in the bare ZnPc regions, indicate a relatively weak ZnPc-C60 interaction. As the C60 deposition continues, the C60 molecules finally cover the surface with a monolayer film, forming an abrupt C60-ZnPc interface. Molecularly resolved STM images reveal quasi-linear C60 chains of one molecule wide (Fig. 4.4 left). Chains exhibit one nm nearest-neighbor separation and branch lengths of ~5.12 nm. The C60 density, ~0.54 molecules/nm2, is less than half that of close-packed C60 films. After thermal annealing at 150oC for 10 min, some coarsening of the chain phase is observed: some close-packed C60 clusters are enfolded within wandering chains. These meandering chains are seen on several organic molecule surfaces like ZnPc, Pentacene and ?-Sexithiophene [9, 10] (Fig. 4.4). This suggests this is a more general phenomenon involving organic overlayers and we will try to provide some general theoretical insight into these patterns. 78 (a) C60/Pentacene/Ag(111) (b) C60/ZnPc/Ag(111) (c) C60/?-Sexithiophene/Ag(111) Fig. 4.4 C60 meandering chains formed on different organic surfaces. (a) C60/Pentacene/Ag(111). (b) C60/ZnPc/Ag(111). (c) C60/?-Sexithiophene/Ag(111) These patterns are distinctly different from the C60 deposition patterns on bare Ag, Cu, Au or Si surfaces, where C60 molecules normally form close packed structures [2] as shown in Fig. 4.5. The different behavior of C60 molecules on organic surfaces suggests more complicated effective interactions between the C60 molecules, affected by their presence on the bilayer surface. Here we will discuss some of the possible effective interactions between the C60 molecules, and try to reproduce the experimentally measured patterns. Fig. 4.5 C60 close packed structure on Ag(111) 79 4.2 Theoretical Model of Molecular Interactions The simplest interaction between two C60 molecules is the attractive van der Waals interaction between the carbon atoms on the different molecules. The integrated effect of all these interactions leads to a quite strong effective attractive interaction between two C60 molecules at a separation of about 10 ? with a strong repulsion at shorter separations from overlap of the repulsive carbon cores. This basic interaction is present even in the gas phase, and is well captured by the Girifalco potential, as described below [11]. This fundamental interaction will continue to play a key role even in condensed phases and indeed the spacing of adjacent molecules in the C60 chains is about 10?. Previous examples of the spontaneous formation of linear chain structures in other chemical systems have been explained by many different effects: polymerization of dipolar spheres [ 12 ], Friedel oscillations [ 13 ], anisotropic long-ranged interactions between atoms [14], and dipole-dipole repulsions [15,16]. Here we suggest that the most likely new feature in the effective interaction between adsorbed C60 molecules on organic bilayer substrates is electrostatic in origin, associated with the large charge transfer arising in such systems. This of course is why these systems are of interest as models of photovoltaic devices. As we argue below, electrostatic interactions seem to be the only interaction strong enough to compete with the basic Girifalco attraction and thus permit interesting pattern formation. 80 In general, at the interfaces of two materials, there can be charge transfer and it is possible to form large dipoles [17,18]. This dipole could be due to polarization effects [19] and charge transfer is certainly expected from donor to acceptors. Existing DFT calculations have shown that the dipole formed at electron donor and acceptor interfaces can be significant [20]. Another possible interaction affecting C60 surface patterns arises from the non-electrostatic corrugations in the substrate, here reflecting the size, shape and arrangements of the ZnPc molecules in the organic overlayer, as expressed mainly by short-ranged van der Waals interactions between C60 molecules and the ZnPc overlayer. Before discussing the (more important but more complex) electrostatic interactions, we will study the effects of the van der Waals interactions alone between the C60 molecules themselves and between the C60 and the ZnPc using a standard open-source molecular dynamics package, DL_POLY [21]. We will see that using the van der Waals interactions alone generates only close-packed 2D C60 islands, with no chain formation, as expected from the strong Girifalco intermolecular C60 interaction potential. Moreover the DL_POLY simulation based on the summation of all atom-atom interactions is very time consuming. So, later in this chapter, we will analyze some other possible effective interactions and develop an empirical model for more efficient simulations. 4.2.1 van der Waals Interactions Here, we use a standard software package DL_POLY [21] to perform a classical 81 molecular dynamics (MD) simulation of the C60/ZnPc system. DL_POLY was developed initially for simulation of molecular and biomolecular condensed phases using classical interatomic potentials between all atomic sites of the molecules. This calculation is similar to another C60 simulation on the GaAs(001) surfaces [22], which suggested charge transfer should also be considered in the calculation. We do the simulation within a monoclinic box, employing periodic boundary conditions in horizontal (x,y) directions. At both the top and bottom of the box, we put one monolayer of fixed ZnPc molecules (no silver substrate, since we consider only short-ranged van der Waals interactions here) based on the experimental data (Table 4.1) and use two ZnPc film as boundary walls on top and bottom of the simulation box (Fig. 4.6). The carbon atoms of C60 will interact with themselves and also with the atoms in ZnPc. Within one C60 molecule, all the carbon atoms interact with their neighbors with harmonic bonds: Vij = ?0(??? ??0)2 , where ?0=4.06x103kJ/mol when r0=1.46? for long bonds and ?0=5.94x103kJ/mol when r0=1.385? for short bonds. Between different molecules, we use the Lennard-Jones potential ??? = 4??? [(???? )12 ?(???? )6] to characterize the interactions between atoms, where i, and j correspond to two different atoms. There are standard parameters for each type of atom reflecting the van der Waals interaction as part of the simulation package. We use the potential parameters in Table 4.1 and the usual mixing rules given in Eq. (5.1) to describe our system [23,24]. The total system is described by a canonical ensemble (N,V,T) and we use Nos?-Hoover [25] thermostat to keep the system temperature constant. 82 ? (?) ? (10J/mol) C 3.471 39.81 H 2.846 6.363 N 3.262 32.23 Zn 4.045 23.02 Table 4.1 Parameters for Lennard-Jones interactions ??? = 12(?? +??) ??? = (?? ???)1 2? (4.1) The atoms in the ZnPc molecules are fixed at the experimental geometry depicted in Fig. 4.3, and only offer a rigid substrate potential to the mobile C60 molecules. The simulation begins with 18 C60 molecules put randomly in the box with random velocities in all directions (Fig. 4.6 left). After a while, because of the Nose Hoover dynamics, the system relaxes to its equilibrium configuration where the C60s deposit on both the top and bottom ZnPc surfaces, and aggregate into close-packed structures (Fig. 4.6 right). This indicates that the van der Waals interaction by itself will not generate meandering chains but rather form compact islands. Since ZnPc and C60 are considered to be electron donors and acceptors, respectively, we expect charge transfer, but these electrostatic interactions are completely ignored by the standard LJ potentials used. Also, this simulation with DL_POLY is very time consuming because it calculates all the individual atomic 83 forces, velocities, and coordinates in the system. Therefore in all the work discussed henceforth we will not use DL_POLY. Instead, we will create a simplified two dimensional model for the C60/organic molecules system that allows for more general electrostatic intermolecular interactions as well and we will carry out MD simulations of these 2D models using a special in-house MD program written by Qiang Liu. (a) (b) Fig. 4.6 DL_POLY simulation. (a) Initial configurations of the system in a (58?)x(58?)x(200?) box with ZnPc monolayer surface as top and bottom. (b) C60?s on bottom surface in final configurations. 4.2.2 Intermolecular interaction VG: Girifalco Potential First, to simplify our simulation, we are not going to use the individual atom-atom van der Waals interaction between two C60 molecules, instead, we can choose some effective potential between two whole C60 molecules. The Girifalco potential [5] is one of the simplest ways to describe C60-C60 interactions and 84 captures well most important features of interest here. If we only think about the van der Waals interactions, then the interaction between two C60s can be seen as the summation of the L-J interactions between two sets of 60 carbon atoms ??(?) = ? (? ?? ?? 6 + ? ???12) 60 ?,?=1 , where A and B are L-J parameters, i and j are two carbon atoms in two different C60 molecules. Following Girifalco, we assume all these atoms are uniformly distributed on the C60 sphere surface. Then the potential between C60 molecules 1 and 2 can be written as (Eq. 4.2): ??(?) = ?2 ???1??2(? ?? ?? 6 + ? ???12) (4.2) S1 and S2 are two C60 molecules? surfaces, r1,2 is the distance between two points on each of the C60 surface. ? is the Carbon density distribution: ? = 60 ?1,2? , and ?(?12) = ?2 (? ?? ?? 6 + ? ???12) is the Lennard-Jones potential between these two points with distance rij, and the integration is over two C60 sphere surfaces S1 and S2. After carrying out the integration, we get Eq. 4.3 ??(?) = ??[ 1?(??1)3 + 1?(?+1)3 ? 2?4]+?[ 1?(??1)9 + 1?(?+1)9 ? 2?10] ? = ?2? (2? = 7.1?) ? = 4.68?10?2??,? = 8.48?10?5?? (4.3) Here r is the distance of the centers of the two C60 molecules. Upon comparing the integrated Girifalco intermolecular potential to an individual atomic C-C LJ potential (Fig. 4.7), one can see that the Girifalco potential has a much lower energy minimum, almost 100 times lower than that of the atomic C-C interaction. This 85 strong interaction makes the C60 molecules stay together to form close packed structures on most metal surfaces. Thus, to destabilize the close packed structure, there must be some other strong interaction on organic surfaces to compete with the Girifalco interactions. Another important feature about the Girifalco potential is that it has a very narrow potential well compared to LJ potential with same potential well and ? (Fig. 4.8). Fig. 4.7 Girifalco Potential and the C-C LJ potential. The black line is Girifalco potential. The red one is the L-J potential of two Carbon atoms. -0.003eV 86 Fig. 4.8 Comparison of Girifalco potential (blue) and LJ potential (red) with same ? and ?. The Girifalco potential has a harshly repulsive core and is very short -ranged. 4.2.3 Intermolecular interaction Vd: dipole-dipole repulsion One possible way to destabilize the close packed structure is to add repulsive forces between the C60s. This is a possibility because there is evidence for charge transfer in both experiments [7,8] and theory [26]. S.J. Singer [15] showed that short ranged attraction and long ranged dipole repulsion could lead to the formation of striped domains qualitatively resembling the C60 patterns seen in experiments (Fig. 4.9), as was shown using a lattice gas model with the Hamiltonian below (Eq 4.4). 87 Fig. 4.9 (by S.J.Singer) Stripes formed by dipoles with vdW interactions. This Hamiltonian contains two terms, one is the nearest neighbor attraction between two occupied sites with parameter A and the other is a dipole-dipole repulsion from a vertical dipole with strength given by parameter B. nR and nR? represent the occupancy status at positions R and R?, and takes values 0 or 1. ? = ?? ? ????? +?2 ? ?????|? ???|3 (4.4) The stripes form because for appropriate values of B, when there are only few neighboring sites occupied, another adatom at the boundary nearby will feel the attraction from the A term much stronger than the longer-ranged repulsion, so the cluster can grow. But when the width of the cluster is large enough, the integrated repulsion on an additional atom will be large enough to push it away. Stripes form as a result of the competition between short-ranged attraction and longer-ranged dipole repulsions. But in our case with C60 molecules, the ?stripe? width is only one molecule. To prevent formation of close packed structures, and form single molecule chains in a 88 model like this, we show below that there must be a very large repulsion. But achieving such a large dipolar repulsion may not be as difficult as one might at first suppose when considering large molecules like C60. Suppose we have a dipole that comes from a one electron charge separation of 3.5 angstroms (a typical distance from the bottom of a C60 molecule to the ZnPc molecule is around 3 Angstroms): ? = 3.5(? ??) = 16.8(?). Then, two parallel dipoles will have a potential energy ?? = ? 2 4???3 = 127 ?3 (??) (4.4) where r is the distance between the two C60 molecule centers in Angstroms. When r=10? (equilibrium distance of neighbor C60 molecules both from experimental results and Girifalco potential calculations), ?? = 0.127(??), which is comparable to the strength of the Girifalco potential. By empirically adding a vertical dipole repulsion of varying strength to the Girifalco attraction, we can change the total effective potential between two C60 molecules Veff=VG+V d, which is shown in Fig. 4.10. In Fig. 4.10 we see that a dipole of 16.8D can significantly reduce the Girifalco potential depth and also give a longer ranged repulsive potential that could favor formation of chains. When this dipole is increased to 20.6D, there is still a negative local minimum energy at r=10?. But a dipole of 23.8D is so large that the local minimum energy at r=10? becomes positive, which could make all the C60 molecules prefer to stay far away from each other and so is not consistent with experiment. 89 Fig. 4.10 Effective C60-C60 interaction: Girifalco potential plus dipole-dipole repulsion. 4.2.4 Charge-Charge repulsion Proceeding strictly empirically for the moment, we could imagine another scenario where the repulsive force between C60 molecules arises from a negative charge associated with each C60 molecule generated by the negative charge transfer from ZnPc to C60 molecules. The neutralizing positive charges remaining in the ZnPC layer could be most simply modeled as a uniform 2D neutralizing background layer, which then has no direct effect on the effective force between C60 molecules. In effect then we replace the vertical dipole in the previous model with an effective negative point charge of appropriate magnitude. As we will see this model too can 90 lead to single molecule chains under some conditions. The Coulomb interaction is much more long ranged than the dipole one since the potential is proportional to 1 ?? , instead of the leading term in dipolar interactions 1 ?3? . But we are only interested in the short-ranged effects of these Coulomb repulsions on the nearest-neighbor length scale of the Girifalco attractions. Thus it seems reasonable to use a short-ranged truncation of the Coulomb interaction at some larger distance instead of formally dealing with the full Coulomb repulsion. A useful truncation method involves modifying the 1/r potential from a ?-function point charge using a neutralizing Gaussian charge distribution with appropriately chosen width ?. [27] Suppose we have charge density distribution function from an array of point charges ?(?) = ???(?) ? = ????(? ???) ? (4.5) The potential field generated by this charge distribution is the solution of Poisson?s equation ?2?(?) = ??(?)? 0 (4.6) and the solution is ?(?) = 14?? 0 ? ?(?)|? ???|?3?? (4.7) Now we split each delta function ?(? ???) of a single point charge ?? into two 91 terms by adding and subtracting a Gaussian distribution function ??(?) = 1?2??2 ?? ?2 2?2 (4.8) and the result is: ??(?) = ???(?)+???(?) ???(?) = ??[?(? ???)???(? ???)] ???(?) = ????(? ???) (4.9) The potential field generated by ?? can be split in the similar way ??(?) = ???(?)+???(?) ???(?) = ??4?? 0 ??(? ???)???(? ?? ?) |? ???| ? 3?? ???(?) = ??4?? 0 ???(? ?? ?) |? ???| ? 3?? (4.10) The potential field generated by a Gaussian charge distribution ??(?) also satisfies, Poisson equation ?2??(?) = ???(?)? 0 (4.11) and is given by ??(?) = 14?? 0? erf( ??2?) (4.12) where the error function erf(?) = 2??? ???2???0 . Therefore, 92 ???(?) = 14?? 0 ?? |? ???|erfc( |? ???| ?2? ) ???(?) = 14?? 0 ?? |? ???|erf( |? ???| ?2? ) (4.13) and erfc(?) = 1?erf (?) is complementary error function (Fig. 4.11). Fig. 4.11 Error function and complementary error function The initial 1? function, separated into short-ranged and long-ranged parts is shown below in Fig. 4.12: 93 Fig. 4.12 Long-ranged and short-ranged Coulomb interactions with ? = 5 We also note that there is not much difference between the repulsive force from the short-ranged part and total potential at distances less that ?? as shown in the Fig. 4.13 below. 94 Fig. 4.13 Forces of short-ranged and whole Coulomb interactions with ? = 5 4.2.5 Substrate potential The short-ranged C60-ZnPc substrate interaction could also affect properties of the meandering chains, as briefly mentioned before in Sec. 4.2. Rather than rely on the crude model considered there that considers only van der Waals interactions from each atom, we decided to create an effective empirical substrate potential based on the symmetry of the unit cell and experimental results, as shown in (Eq. 4.14): ????(?,?) = ?10 { ?[cos(4??14.4)+1]?[???(4??14.4)+1] ?0.2[???(4??14.4+?)+1]?[???(4??14.4+?)+1] ?1.5[???(2??14.4)+1]?[???(2??14.4)+1] } (4.14) 95 ????_??? ?????_??? ? ? (4.15) ? is a parameter with units of energy, roughly equals to the maximum potential difference on the surface. We can see the shape of the potential in Fig. 4.14 below with ?=6.0 eV. This potential has the periodicity similar to experiments, and its shape repels C60?s on Zn atom sites, which we assume are at the higher energy location based on experiment. Fig. 4.14 Part of the manipulated substrate potential ???? with ?=10.0 eV Now, the overall system Hamiltonian can be written as: ? = ? ????(??) ?=1,? +?[??(???)+??(???)] (4.16) ??? 96 where i and j are C60 molecules tags. 4.3 Simulation method: Langevin Dynamics Now we can do the molecular dynamics simulation with our simplified intermolecular potential model using a special MD program we wrote. Here we use Langevin dynamics [28], with governing Equation 4.17: ? ? 2 ??2 ? = ???? ??? ? ??? +? ?????+?? = 2??????(?) (4.17) x symbolically represents the position coordinate of one C60 molecule. R is a random force that represents thermal fluctuations, which inputs energy into the system constantly. ? is the friction constant, which takes kinetic energy away from the system. They work together to keep the system fluctuating around the temperature T. Ftot is the total force acting on the C60 molecule which comes from other molecules, including the Girifalco potential, dipole-dipole repulsion and the substrate potential. For efficiency, here we use a two dimensional canonical ensemble, in which temperature T, system area S and number of C60 molecules N are fixed. The ZnPc molecules are treated as a fixed substrate potential. The most important parameters to be determined are the vertical dipole d, charges on C60 molecules, and the substrate potential strength ?. 97 4.4 Simulation Results 4.4.1 The substrate effect First, we will study the situation where there is no dipole repulsion. We set the substrate potential strength ?=0.24eV and 0.48eV, so that ????_??? ?????_??? ? ? = 0.24?? and 0.48??, which is comparable to the Girifalco potential depth and van der Waals interaction strength. This is already an artificially strong substrate interaction, but the final configuration will still have close packed domains because there is no repulsion in regions away from the Zn atomic sites to compete with the strong Girifalco attraction (Fig. 4.15). (a) strength ?=0.24eV (b) strength ?=0.48eV Fig. 4.15 Simulation results with substrate potential. 685 C60 molecules on a (356?)x(356?) substrate, with dipole d=0. Picture taken after 1.6x104 ps of simulation. 4.4.2 The dipole-dipole repulsion effect Second, let?s assume there is no substrate interaction, but only the C60-C60 98 interactions including dipole-dipole repulsion. Fig. 4.10 shows that this dipole has to be at least 16.8D to increase the effective C60- C60 potential by ? 0.13eV when they are 10.0? from each other. We can see from Fig. 4.16 that, as the dipole d grows from 16.8D to 23.8D, the close packed structure disappears and the chains become more dominant. But the chain structures are not rigid, they change shape slowly all the time, unlike what is seen in the experiment. (a) d=16.8D (b) d=19.2D (c) d=20.6D (d) d=23.8D Fig 4.16 Final configurations with different dipole moments. No substrate potential, ?=0. Snapshots are taken 1.6x104 ps after simulation started. 99 4.4.3 Considering both dipole repulsions and C60-substrate interactions. The real situation should include both molecule-substrate interactions and dipole-dipole repulsions. The former can reduce fluctuations of chain structures if they form and the latter is the key interaction that destabilizes the close packed structures due to the strong Girifalco potential. The C60 molecular patterns formed on a periodic substrate potentials have greatly reduced fluctuations, similar to what is seen in the experiment. Some of the simulation results are shown here in Fig. 4.17. This shows that it is at least possible to reproduce many aspects of experimental patterns with a proper combination of dipolar and substrate interactions, though the magnitudes chosen have yet to be justified theoretically. 100 Fig. 4.17 Simulation results of 685 C60s on a 356?x356? surface at 300K, after 1500ps of simulation, with dipole d=19.6D, ?=0.16eV Similar patterns can also be seen in a C60 molecular system with properly chosen charges on C60?s, as shown below in Fig. 4.18. The Coulomb interactions are truncated with ?=50?, which is five times of the C60?s nearest distance. These results suggest that destabilization of the close packed structures by longer ranged repulsive forces may be a general and robust feature, independent of the precise details of the model. Fig. 4.18 Simulation results of 685 C60?s on a 356?x356? surface at 300K, after 1500ps of simulation, with 0.4e charge on C60 and the Gaussian truncation with ?=50?. Substrate interaction strength ?=0.16eV 101 4.5 Conclusions In this chapter, we analyzed the possible interactions between C60 molecule themselves and with the substrate. The basic van der Waals intermolecular interactions between two C60 molecules can be accurately represented as an effective Girifalco pair potential. The strong van der Waals attraction is the reason close packed clusters for C60?s usually form on most metals and other surfaces like Si. To form the unusual single chain structures on organic monolayers as ZnPc, Pentacene and ?-Sexithiophene, we need an additional longer-ranged, strong repulsive force between C60 molecules. This could be generated, by a large dipole of 20 Debye or a charged C60 model with a charge of 0.4e. These new contributions could arise because the C60 is good electron accepter while as the ZnPc is good electron donor. The C60 molecules on metal surfaces can also have large charge transfer from substrate, but the strong screening from metal probably greatly weakens the repulsion effect, leading to close packed cluster formation. We believe the charges in the organic ZnPc layer are less mobile and thus less able to screen the charge-induced effective C60 interactions. In principle microscopic quantum calculations could determine whether these speculations have merit. In the next chapter, we will use DFT calculation to try to estimate the detailed charge transfer and charge density distributions of the C60/ZnPc/Ag(111) and C60/Ag(111) system. References: 102 1 T. Taima, J. Sakai, T. Yamanari, and K. Saito, ?Doping effects for organic photovoltaic cells based on small-molecular-weight semiconductors?, Sol. Energy Mater. Sol. Cells, 93, 6-7, 742 (2009) 2 M.R.C.Hunt, S.Modesti, P.Rudolf, R.E.Palmer ?Charge Transfer and Structure in C60 Adsorption on Metal Surfaces?, Phys. Rev. B 51, 10039 (1995) 3 L.Wang, H.Cheng, ?Rotation, Translation, Charge Transfer, and Electronic Structure of C60 on Cu(111) Surface?, Phys. Rev. B 69, 045404 (2004) 4 L.Wang, H.Cheng, ?Density Functional Study of the Adsorption of A C60 Monolayer on Ag(111) and Au(111) Surfaces?, Phys. Rev. B 69, 165417 (2004) 5 L.A.Girifalco "Molecular Properties of Fullerene in the Gas and Solid Phases", J. Phys. Chem. 96, 858-861 (1992) 6 Z.R.Hong, B.Maennig, R.Lessmann, M.Pfeiffer, and K.Leo "Improved Efficiency of Zinc Phthalocyanine/C60 Based Photovoltaic Cells via Nanoscale Interface Modification", Appl. Phys. Lett. 90, 203505 (2007) 7 H.D.Burrows, A.A. Kharlamov, ?About Energy and Electron Transfer Processes in C60/Phthalocyanine Films? Progr Colloid Polym Sci 123, 52 (2004) 8 T.Toccoli, A.Boschetti, S.Lannotta, ?Molecular Materials for Optoelectronics by Supersonic Molecular Beam Growth: Co-deposition of C60 and ZnPc?, Synthetic Metals 122-1, 229 (2001) 9 W.Jin, Q.Liu, J.D. Weeks and J. Reutt-Robey, in preparation 10 H.L.Zhang et al., ?C60 Molecular Chains on ?-Sexithiophene Nanostripes? Small, 3, 2015-2018 (2007) 11 M. C. Abramo, C. Caccamo, D. Costa, G. Pellicane and R. Ruberto, ?Atomistic Versus Two-body Central Potential Models of C60: A Comparative Molecular Dynamics Study?, Phys. Rev. E 69, 031112 (2004) 103 12 J. Stambaugh, K. V. Workum, J. F. Douglas and W. Losert, ?Polymerization Transitions in Two-dimensional Systems of Dipolar Spheres?, Phys. Rev. E 72, 031301 (2005) 13 K. H. Lau and W. Kohn, ?Indirect Long-range Oscillatory Interaction between Adsorbed Atoms?, Surface Science 75, 69-85 (1978) 14 S. Koh and G. Ehrlich, ?Self-Assembly of One-Dimensional Surface Structures: Long-Range Interactions in the Growth of Ir and Pd on W(110)?, Phys. Rev. Lett. 87, 106103 (2001) 15 M.M. Hurley and S.J. Singer, ?Domain-array Melting in the Dipolar Lattice Gas? Phys. Rev. B 46 (9), 5783 (1992) 16 A.D. Stoycheva and S. J. Singer, ?Computer Simulations of A Two-dimensional System with Competing Interactions?, Phys. Rev. E 65, 036706 (2002) 17 C. Stadler, S. Hansen, I. Kroger, C. Kumpf and E. umbach, ?Tuning Intermolecular Interaction in Long-range-ordered Submonolayser Organic Films?, Nature Physics, 5, 153 (2009) 18 M. A. Baldo and S. R. Forrest, ?Interface-limited Injection in Amorphous Organic Semiconductors?, Phys. Rev. B 64, 085201 (2001) 19 M. Linares, D. Beljonne, J. Cornil, etc. ?On the Interface Dipole at the Pentacene-Fullerene Heterojunction: A Theoretical Study?, J. Phys. Chem. C 114, 3215-3224 (2010) 20 I. Avilov, V. Geskin and J. Cornil, ?Quantum-Chemical Characterizaiton of the Origin of Dipole Formation at Molecular Organic/Organic Interfaces?, Adv. Func. Mater. 19, 624-633 (2009) 21 W. Smith, C.W. Yong and P.M. Rodger, ?DL_POLY: Application to Molecular Simulation?, Molecular Simulation 28, 385 (2002) 104 22 H. Kamiyama, K. Ohno, and y. Kawazoe, ?Classical MD Simulation of C60 Adsorbed on GaAs(001) Surface?, Sci. Rep. RITU A41, 187-190 (1996) 23 S.Mayo, B.D.Olafson, W.A.Goddard, ?Dreiding: A Generic Force Field for Molecular Simulations?, J. Phys. Chem 94(26) 8897 (1990) 24 T. A. Halgren, ?Pepresentation of van der Waals (vdW) Interactions in Molecular Mechanics Force Fields: Potential Form, Combination Rules, and vdW Parameters?, J. Am. Chem. Soc. 114, 7827-7843 (1992) 25 S. Nose, ?A Unified Formulation of the Constant Temperature Molecular Dynamics Methods?, J. Chem. Phys. 81, 511 (1984) 26 H. Mizuseki, N. Igarashi, R.V. Belosludov, A.A. Farajian, Y. Kawazoe, ?An Organic Molecule Fullerene Mixture for A High Efficiency Photovoltaic Device: A Theoretical Study?, Nanotech, 2 (2003) 27 J.M. Rodgers and J.D. Weeks, ?Accurate Thermodynamics for Short-ranged Truncations of Coulomb Interactions in Site-site Molecular Models?, J. Chem. Phys. 131, 244108 (2009) 28 W.T. Coffey, Yu.P. Kalmykov, J.T. Waldron ?The Langevin Equation with Applications to Stochastic Problems in Physics, Chemistry and Electrical Engineering? World Scientific Series in Contemporary Chemical Physics, 14 (2004) 105 Chapter 5 Quantum Chemistry Calculation (DFT) In last chapter we used simulations with empirical potentials to study C60 meandering chains formed in the C60/ZnPc/Ag(111) system. We have shown that an effective longer ranged repulsive dipole or charge interaction added to the strong short ranged van der Waals (Girifalco) attraction between nearest neighbor C60 molecules can destabilize close packed clusters and under some conditions produce chainlike patterns similar to those seen in experiment. But there are still many open questions we might hope to answer. Most notably we argued that charge transfer could generate the strong electrostatic force needed to compete with the Girifalco attraction.. But can we provide some justification for these models based on more rigorous quantum mechanical calculations? Is there sufficient charge transfer to produce the large dipole or charge strength needed and how precisely could charge transfer affect effective C60- C60 intermolecular interactions on organic overlayers? If charge transfer is involved, why do C60 molecules form close-packed islands when deposited directly on Ag substrates, where arguably even more charge transfer could occur? Unfortunately, definitive answers to almost all of these questions for such a complicated system seems well beyond the capacity of even the best quantum treatments available today. And we are not expert in quantum calculations in any case. But we can use available standard state-of-the-art quantum density functional theory (DFT) packages to study some aspects of the charge transfer, and at least 106 check whether the order of magnitude effects we need seem reasonable. This is the subject of this chapter. This represents new and ongoing work and is a much less developed state than the rest of this thesis. 5.1 Comments on density functional theory DFT is the most widely used and powerful method currently available to study quantum effects in extended molecular and solid-state systems. Its power results from its ability to determine ground state electronic properties using in principle only the total electron density, a function of three spatial variables, rather than the total orbital-based N-electron wave function used in standard quantum treatments. That this simplification is theoretically possible was proved in work by Hohenberg and Kohn in 1964 [1]. In practice orbitals still are used in standard DFT treatments to describe the kinetic energy of a simpler model system of N non-interacting electrons with the same total density as that of the full interacting system. The model system is used to approximate the kinetic energy of the full system. Another large component of the total system?s energy is the classical (Hartree) electron-electron potential energy, which is readily expressed in terms of the total electron density. All corrections to these two large terms in the exact total energy are given by the small exchange-correlation density functional, whose exact form unfortunately is not known, but can often usefully be approximated. This key idea turned DFT into a practical computational tool and was proposed by Kohn and Sham in 1965[2]. For 107 this body of work Walter Kohn was awarded the Nobel Prize in Chemistry in 1998. Finding a generally accurate exchange-correlation functional is a focus of much current work in this field. Researchers have shown that with currently available functionals good results can often be found in a wide variety of condensed matter systems, particularly for problems involving charge transfer, as might be expected form the focus on the total electron density [3,4,5]. Here we will use a well-known and standard DFT program: the Quantum-ESPRESSO (open-Source Package for Research in Electronic Structure, Simulation, and Optimization) package, which is fully documented and available for free download at http://www.quantum-espresso.org, to analyze some aspects of the C60/ZnPc/Ag(111) system. Unfortunately, a major problem with almost all current implementations of DFT, including those in Quantum-ESPRESSO, is a poor treatment of van der Waals interactions. Thus we will not be able to answer one of the most basic and fundamental questions raised by our interpretation of the experimental results: how precisely does charge transfer from the organic layer to the C60 affect the intermolecular interactions, which certainly will still have a strong van der Waals attractive component? New functionals that give a better treatment of van der Waals interactions are currently being developed but we strongly doubt that this question can be definitively answered for many years to come. However the classical Girifalco treatment of the van der Waals interactions still seems very 108 reasonable, and we can at least use the quantum package where it can be more trusted to look at other aspects of charge transfer in the C60/ZnPc/Ag system. The Quantum-ESPRESSO package uses DFT with a plane wave basis set and pseudo potentials to calculate the ground state electronic structure. This method in principle permits charge transfer and can determine its magnitude. We will mostly use one part of the program, the Plane Wave Self-Consistent Field (PWscf) package. 5.2 C60/ZnPc system The pseudopotential files H.pbe-van_ak.UPF, C.pbe-van_ak.UPF, N.pbe-van_ak.UPF, Zn.pbe-van.UPF can be downloaded from the website: http://www.quantum-espresso.org. The plane-wave kinetic energy cutoff for the wave functions was 30 (150) Ry. A (2x2) k-point grid was used for summation over the surface Brillouin Zone. 5.2.1 C60/ZnPc Initially we ignored the silver substrate to reduce the computer time required, though in principle this can and should be included, as discussed below. The calculation automatically has periodic boundary conditions in x, y, and z direction. Based on experimental results, I chose a monoclinic box (?=90o, ?=90o, ?=95.0o, a=14.5?, b=13.8?, c=14.5?) as my calculation unit. The C60 molecule center was put initially 4.0? above the ZnPc molecule, at three different horizontal (x,y) positions ?O? (right above Zn atom), ?C? (at the corner of ZnPc), and ?S? (in the 109 middle between two adjacent ZnPc molecules) (Fig. 5.1). (a) ?O? Center of ZnPc (b) ?C? Corner of ZnPc (c) ?S? Side middle of ZnPc Fig. 5.1 Initial configurations for DFT calculation All the atoms in the calculation relax to their final positions O?, C?, and S? (Fig. 5.2), where the total system reaches the local lowest energy. The relative C60-ZnPc interaction energies, EC? - EO? = - 0.13eV, ES? - EO?= - 0.05eV, are consistent with the empirical corrugations we used in previous simulations. (a) ?O?? Center of ZnPc (b) ?C?? Corner of ZnPc (c) ?C?? Side middle of ZnPc Fig. 5.2 Final configurations of DFT calculation of C60/ZnPc We can also calculate the final charge distributions. The final configurations give us the electronic dipoles: dA? = 3.2D, dB?= - 0.8D, dC? = 0.8D, which are much smaller than what we need based on our empirical model discussed in the last chapter. 110 Moreover these results are in qualitative agreement with several experiments. T. Toccoli[6] reported that there is no charge transfer between C60 and ZnPc if they are deposited layer by layer, while there exists some charge transfer for the co-deposited films. H. D. Burrows[7] showed that vapor deposited C60 films on ZnPc on a glass substrate has no charge transfer from ZnPc, although the liquid deposited C60/ZnPc bilayers do exhibit electron transfer from ZnPc to C60. These results were initially rather discouraging. However, none of these experiments were done on a metal surface, and we left out the Ag substrate, in our initial DFT calculation as well. It would be expected that its inclusion should permit much larger charge transfer and hence the formation of larger dipoles. Thus we next considered simulations of the full C60/ZnPc/Ag(111) system. Of course this greatly increased computer time needed, and many results are only now becoming available. We expect much more rapid progress in the immediate future, since we can now run on the newly added Theoretical Chemistry nodes of the Deep Thought computer cluster at the University of Maryland, which became available for use only on Nov 16th, 2010. 5.2.2 C60/ZnPc/Ag(1 1 1) The 2-D islands formed on metal surface are due to the strong C60-C60 Girifalco interactions [8,9]. Although there are charge transfer from metal substrate to the adsorbed C60 molecules [3,4], any charge-charge or dipole-dipole interactions between C60 molecules likely get screened by the electron gas on metal surfaces [10]. 111 C60 molecules on organic surfaces with less mobile free electrons may not get screened by electron gas, but on the other hand, they may not get as large electron transfer from the organic substrate as the metal one. So one important question we ask is, how much charge transfer and how large the dipole would be, on the fullerene organic bilayer? There are some experiments results that try to answer this question: T. Toccoli[11] report that there is no charge transfer between C60 and ZnPc if they are deposited layer by layer, while there exist charge transfer for the co-deposited films. H. D. Burrows [7] showed that on a glass substrate, vapor deposited C60 films on ZnPc have no charge transfer from ZnPc while the liquid deposited C60/ZnPc bilayers do. But none of these experiments are done on a Ag surface, so we try to study the full system in what follows. Here, I am using a quantum chemistry software Quantum ESPRESSO package to do the first principles electronic structure calculation of the C60/ZnPc/Ag(111) system. The plane-wave kinetic-energy cutoff was chosen to be 30(150) Rydberg (Ry). A (2x2) k-point grid was used for summations over the surface Brillouin zone. The unit super cell contains one ZnPC, one C60 and three layers of Ag(111) with ( 5x3?3 ) periodicity. The box (14.4?x15.0?x26.0?) has periodic boundary conditions in all x, y, z direction (Fig. 5.3). In z direction, slab replicas were separated by 10?. In x, y direction, ZnPc forms a 2D film, while C60?s are isolated with a distance of 14.4 ?. 112 Fig. 5.3 Unit cell for calculation The ZnPc molecules are put 2.7A above the three layers of Ag(111) surface[12]. The C60 is put at O(center), S(side) and C(corner) three different horizontal positions with respect to ZnPc (Fig. 5.4), and 2.4 to 3.7 ? above ZnPc. 113 Fig. 5.4 Different positions of C60 on top of ZnPc For the ?C? position, we take the total charge distribution after placement of C60, and subtract the isolated C60 and ZnPc/Ag(111) charges distributions. The differential charge distribution is shown in Fig. 5.5. We can see clearly that there is a significant charge transfer from ZnPc/Ag(111) to C60. 114 Fig. 5.5 Differential charge density isosurfaces. Above: ?C60 ZnPc Ag?? ??C60 ??ZnPc ??Ag, below:?C60 ZnPc Ag?? ??C60 ??ZnPc Ag? . Red is negative charge, Blue is positive Charge. Isosurface = ? 0.2e/?3 115 If we integrate the charges along x and y direction, we can get the projected charge distribution along z direction as shown in Fig. 5.6. The integrated charge in the C60 part gives us a total charge transfer of 0.4e. The dipole formed in this system is 20 Debyes. Note that both these values are in very good agreement with the empirically determined values needed in the previous chapter in the dipolar or charged C60 models to produce chain-like patterns! Unfortunately, this calculation cannot resolve the final effect of charge transfer on the effective C60 interactions because of the poor treatment of van der Waals interactions in DFT. But it does suggest that the magnitudes are at least possible. Fig. 5.6 Differential charge distribution along Z directions for C60/ZnPc/Ag(111) system. C60 at ?C? position, 2.8? above ZnPc. 116 If we put C60 at different positions ?O?, ?C?, and ?S?, and from 2.4 to 3.7? above ZnPc, we still get the similar results as shown in Fig. 5.7: (a) (b) (c) (d) Fig. 5.7 Differential charge distribution isosurfaces: ??60 ???? ???? ???60 ? ????? ??? . Isosurface = ? 0.2e/?3. Red is negative charge, blue is positive charge. (a) 2.4? from bottom of C60 to ZnPc center. (b) 3.7? from bottom of C60 to center of ZnPc. (c) 3.0 ? from bottom of C60 to corner of ZnPc. (d) 3.0 ? from bottom of C60 to ?S? point of ZnPc 117 In addition to calculating the z direction differential charge distribution, we can apply another method widely used in quantum chemistry calculations, to calculate the transferred charge: the Bader analysis [13], Using this method we get the exact same result: the charge transfer to C60 is 0.4e. In summary, there is charge transfer to C60?s on ZnPc/Ag(111) and the formation of an interface dipole. As shown in Chapter 4, the dynamics of C60?s on ZnPc surface could be empirically modeled using both charge-charge and dipole-dipole repulsions. 5.2.3 C60/Ag(111) Since C60 on Ag should have a larger charge transfer of 0.6e as discussed in [3,4], why do they only form the close packed structures? Below in Fig. 5.8 is the differential charge distribution of C60/Ag(111) system, where we can see complicated isosurfaces. Although the Bader analysis tells us it does have a 0.6e charge transferred to the C60 molecule, the deformed Ag electron gas below C60 probably screens out any induced C60 charge-charge or dipole-dipole repulsions. We plan to investigate this important question further. 118 Fig. 5.8 C60/Ag(111) system differential charge distribution isosurfaces: ??60 ??? ? ??60 ????. Red is negative charge, Blue is positive Charge. Isosurface = ? 0.2?3. Charges for each part: Q(C60)=0.65e, Q(Ag)=-0.65e 5.2.4 C60/ZnPc/Au(111) Since the Au is not a good charge donor as the Ag, we would expect a smaller charge transfer from ZnPc/Au(111) to C60 than ZnPc/Ag(111). But experiments showed a similar pattern of C60 single chains formed on ZnPc/Au(111) surface. So we will calculate the charge distribution of this system in the future. 5.3 Discussion and Conclusions The C60 molecules on ZnPc surfaces showed meandering chain behavior, which is very different from the close packed structures seen on pure metal surfaces. These 119 chains indicate the presence of strong competing interactions other than the intermolecular C60 Girifalco potential. Two possible interactions, the substrate and dipole-dipole repulsions are investigated for this purpose. The substrate potential alone is not sufficient to prevent nucleation of close-packed structures but is very helpful in stabilizing the added molecules in relatively rigid chain-like structures. The dipole or charge repulsion directly competes with the Girifalco interactions and destabilizes the close-packed structures, but a large charge transfer is needed. This is consistent with results for the deposition of C60 molecules on hexagonal boron nitride (h-BN) monolayer on Ni(111). The h-BN layer a good insulator which prevents charge transfer to C60 molecules [14], and the C60?s indeed form only close packed structures as expected [15]. This study indicates that the experimental patterns can be reproduced by a strong interaction with the organic substrate and a large charge transfer leading to a intermolecular repulsion. We expect some more detailed experimental measurements or quantum chemistry calculations can unveil the underling physics that could lead to such effects. We also hope that by understanding the formation of these chains it can help us to tune the domain shapes and allow the systematic design of organic thin film systems. Appendix 1. C60/ZnPc/Ag(111) Quantum Espresso input file Here is on example of Quantum Espresso calculation input file for pw.x: 120 &CONTROL calculation = 'nscf', restart_mode = 'from_scratch' , outdir = './' , pseudo_dir = '$pseudo_potential_dir$' , prefix = 'c60znpcag' , / &SYSTEM ibrav=14 , celldm(1)=27.3 , celldm(2)=1.0392 , celldm(3)=1.8 , celldm(4)=0.0 , celldm(5)=0.0 , celldm(6)=-0.0 , nat=207 , ntyp=5 , ecutwfc=30 , ecutrho=150 , occupations='smearing', 121 smearing='marzari-vanderbilt', degauss=0.06 / &ELECTRONS diagonalization='cg', / &IONS pot_extrapolation= 'second_order', wfc_extrapolation= 'second_order', / ATOMIC_SPECIES N 14.01 N.pbe-van_ak.UPF H 1.01 H.pbe-van_ak.UPF Zn 65.38 Zn.pbe-van.UPF C 12.01 C.pbe-van_ak.UPF Ag 107.87 Ag.pbe-d-rrkjus.UPF ATOMIC_POSITIONS {angstrom} N 5.70248 6.25271 10.51454 1 1 1 N 6.91789 4.13978 10.47263 1 1 1 N 8.51016 5.99139 10.46524 1 1 1 122 N 3.80377 7.86605 10.55521 1 1 1 N 7.52861 10.87346 10.52738 1 1 1 N 5.93635 9.02184 10.53477 1 1 1 N 8.74403 8.76052 10.48547 1 1 1 N 10.64265 7.14723 10.44475 1 1 1 H 4.60886 2.21602 10.33565 1 1 1 H 2.31127 1.90989 10.29899 1 1 1 H 1.61228 12.45721 10.41641 1 1 1 H 12.89006 9.09379 10.25334 1 1 1 H 13.64804 11.37432 10.51702 1 1 1 H 11.03406 1.08093 10.61893 1 1 1 H 1.55645 5.91944 10.74667 1 1 1 H 12.42930 4.90450 10.41663 1 1 1 H 3.41245 13.93231 10.38107 1 1 1 H 9.83765 12.79722 10.66436 1 1 1 H 8.76535 1.88657 10.69821 1 1 1 H 2.01712 10.10877 10.58333 1 1 1 H 12.83423 2.55602 10.58360 1 1 1 H 12.13528 13.10334 10.70093 1 1 1 H 5.68115 13.12667 10.30180 1 1 1 123 H 0.79846 3.63891 10.48299 1 1 1 Zn 7.22327 7.50666 10.50001 1 1 1 C 3.57539 5.37276 10.43227 1 1 1 C 2.74437 10.71565 10.51089 1 1 1 C 12.32278 9.83919 10.41153 1 1 1 C 10.84723 2.01226 10.60808 1 1 1 C 11.81848 12.21057 10.62846 1 1 1 C 2.49832 12.11519 10.43634 1 1 1 C 9.99362 10.68047 10.61611 1 1 1 C 12.72309 11.15727 10.51771 1 1 1 C 10.40859 12.03825 10.64102 1 1 1 C 4.45289 4.33277 10.38390 1 1 1 C 10.38680 4.71133 10.52712 1 1 1 C 5.10520 11.18338 10.43380 1 1 1 C 3.59928 13.00098 10.39193 1 1 1 C 10.87112 9.64048 10.56774 1 1 1 C 9.88436 6.08316 10.54996 1 1 1 C 9.51945 2.46122 10.63591 1 1 1 C 11.70213 4.29759 10.48912 1 1 1 C 4.56206 8.93012 10.45001 1 1 1 124 C 4.37802 6.67386 10.47130 1 1 1 C 9.34130 3.82986 10.56621 1 1 1 C 10.06848 8.33938 10.52871 1 1 1 C 4.92705 12.55202 10.36410 1 1 1 C 2.12372 5.17405 10.58848 1 1 1 C 6.37581 10.38691 10.39155 1 1 1 C 8.71952 10.10645 10.41432 1 1 1 C 1.72341 3.85596 10.48230 1 1 1 C 11.94819 2.89805 10.56367 1 1 1 C 5.72699 4.90678 10.58569 1 1 1 C 4.05971 10.30191 10.47289 1 1 1 C 2.62802 2.80266 10.37155 1 1 1 C 4.03792 2.97499 10.35899 1 1 1 C 8.07070 4.62633 10.60846 1 1 1 C 13.71476 13.79124 20.33919 1 1 1 C 15.17866 13.79124 20.33919 1 1 1 C 15.87096 14.99044 20.33919 1 1 1 C 15.13896 16.25824 20.33919 1 1 1 C 13.75426 16.25814 20.33919 1 1 1 C 13.02236 14.99034 20.33919 1 1 1 125 C 13.26236 12.68484 19.49409 1 1 1 C 14.44666 12.00114 18.97169 1 1 1 C 15.63096 12.68494 19.49399 1 1 1 C 16.75116 12.83764 18.69449 1 1 1 C 17.05516 15.15184 19.49399 1 1 1 C 15.87086 17.20314 19.49389 1 1 1 C 15.17846 18.09694 18.69449 1 1 1 C 13.71456 18.09694 18.69459 1 1 1 C 13.02226 17.20304 19.49399 1 1 1 C 11.83796 16.51924 18.97169 1 1 1 C 11.83796 15.15174 19.49409 1 1 1 C 11.41006 14.10524 18.69469 1 1 1 C 12.14206 12.83754 18.69469 1 1 1 C 14.44656 11.50704 17.67819 1 1 1 C 13.26226 11.66834 16.83309 1 1 1 C 12.14206 12.31514 17.32719 1 1 1 C 11.41006 13.26004 16.48199 1 1 1 C 10.95766 14.36644 17.32719 1 1 1 C 10.95766 15.65994 16.83309 1 1 1 C 11.40996 16.76634 17.67819 1 1 1 126 C 12.14186 17.71124 16.83299 1 1 1 C 13.26206 18.35804 17.32699 1 1 1 C 17.05506 16.51944 18.97159 1 1 1 C 15.63066 17.34164 14.66589 1 1 1 C 16.75096 17.18894 15.46529 1 1 1 C 17.48296 15.92124 15.46529 1 1 1 C 17.05496 14.87474 14.66589 1 1 1 C 15.87066 15.03614 13.82079 1 1 1 C 13.71436 16.23524 13.82079 1 1 1 C 13.26206 17.34154 14.66599 1 1 1 C 14.44636 18.02534 15.18829 1 1 1 C 14.44636 18.51944 16.48179 1 1 1 C 15.63076 18.35804 17.32699 1 1 1 C 16.75096 17.71134 16.83289 1 1 1 C 17.93526 15.66004 16.83279 1 1 1 C 17.93536 14.36654 17.32689 1 1 1 C 17.48296 13.26014 16.48179 1 1 1 C 17.05506 13.50724 15.18829 1 1 1 C 15.87076 12.82344 14.66599 1 1 1 C 15.13876 13.76834 13.82079 1 1 1 127 C 13.75406 13.76834 13.82089 1 1 1 C 13.02206 15.03604 13.82089 1 1 1 C 12.14186 17.18884 15.46549 1 1 1 C 11.40996 15.92104 15.46549 1 1 1 C 11.83786 14.87464 14.66609 1 1 1 C 11.83786 13.50704 15.18849 1 1 1 C 13.02216 12.82334 14.66609 1 1 1 C 13.71456 11.92954 15.46549 1 1 1 C 15.17846 11.92954 15.46549 1 1 1 C 15.63086 11.66844 16.83289 1 1 1 C 16.75116 12.31524 17.32699 1 1 1 C 17.48296 16.76644 17.67799 1 1 1 C 17.48306 14.10544 18.69449 1 1 1 C 15.17826 16.23524 13.82079 1 1 1 Ag 0.00000 0.00000 7.80000 1 1 1 Ag 0.00000 5.00441 7.80000 1 1 1 Ag 0.00000 10.00883 7.80000 1 1 1 Ag 2.88930 0.00000 7.80000 1 1 1 Ag 2.88930 5.00441 7.80000 1 1 1 Ag 2.88930 10.00883 7.80000 1 1 1 128 Ag 5.77860 0.00000 7.80000 1 1 1 Ag 5.77860 5.00441 7.80000 1 1 1 Ag 5.77860 10.00883 7.80000 1 1 1 Ag 8.66790 0.00000 7.80000 1 1 1 Ag 8.66790 5.00441 7.80000 1 1 1 Ag 8.66790 10.00883 7.80000 1 1 1 Ag 11.55720 0.00000 7.80000 1 1 1 Ag 11.55720 5.00441 7.80000 1 1 1 Ag 11.55720 10.00883 7.80000 1 1 1 Ag 1.44465 2.50221 7.80000 1 1 1 Ag 1.44465 7.50662 7.80000 1 1 1 Ag 1.44465 12.51104 7.80000 1 1 1 Ag 4.33395 2.50221 7.80000 1 1 1 Ag 4.33395 7.50662 7.80000 1 1 1 Ag 4.33395 12.51104 7.80000 1 1 1 Ag 7.22325 2.50221 7.80000 1 1 1 Ag 7.22325 7.50662 7.80000 1 1 1 Ag 7.22325 12.51104 7.80000 1 1 1 Ag 10.11255 2.50221 7.80000 1 1 1 Ag 10.11255 7.50662 7.80000 1 1 1 129 Ag 10.11255 12.51104 7.80000 1 1 1 Ag 13.00185 2.50221 7.80000 1 1 1 Ag 13.00185 7.50662 7.80000 1 1 1 Ag 13.00185 12.51104 7.80000 1 1 1 Ag 1.44465 0.83407 5.44090 1 1 1 Ag 1.44465 5.83848 5.44090 1 1 1 Ag 1.44465 10.84290 5.44090 1 1 1 Ag 4.33395 0.83407 5.44090 1 1 1 Ag 4.33395 5.83848 5.44090 1 1 1 Ag 4.33395 10.84290 5.44090 1 1 1 Ag 7.22325 0.83407 5.44090 1 1 1 Ag 7.22325 5.83848 5.44090 1 1 1 Ag 7.22325 10.84290 5.44090 1 1 1 Ag 10.11255 0.83407 5.44090 1 1 1 Ag 10.11255 5.83848 5.44090 1 1 1 Ag 10.11255 10.84290 5.44090 1 1 1 Ag 13.00185 0.83407 5.44090 1 1 1 Ag 13.00185 5.83848 5.44090 1 1 1 Ag 13.00185 10.84290 5.44090 1 1 1 Ag 0.00000 3.33628 5.44090 1 1 1 130 Ag 0.00000 8.34069 5.44090 1 1 1 Ag 0.00000 13.34511 5.44090 1 1 1 Ag 2.88930 3.33628 5.44090 1 1 1 Ag 2.88930 8.34069 5.44090 1 1 1 Ag 2.88930 13.34511 5.44090 1 1 1 Ag 5.77860 3.33628 5.44090 1 1 1 Ag 5.77860 8.34069 5.44090 1 1 1 Ag 5.77860 13.34511 5.44090 1 1 1 Ag 8.66790 3.33628 5.44090 1 1 1 Ag 8.66790 8.34069 5.44090 1 1 1 Ag 8.66790 13.34511 5.44090 1 1 1 Ag 11.55720 3.33628 5.44090 1 1 1 Ag 11.55720 8.34069 5.44090 1 1 1 Ag 11.55720 13.34511 5.44090 1 1 1 Ag 1.44465 4.17035 3.08179 0 0 0 Ag 1.44465 9.17476 3.08179 0 0 0 Ag 1.44465 14.17917 3.08179 0 0 0 Ag 4.33395 4.17035 3.08179 0 0 0 Ag 4.33395 9.17476 3.08179 0 0 0 Ag 4.33395 14.17917 3.08179 0 0 0 131 Ag 7.22325 4.17035 3.08179 0 0 0 Ag 7.22325 9.17476 3.08179 0 0 0 Ag 7.22325 14.17917 3.08179 0 0 0 Ag 10.11255 4.17035 3.08179 0 0 0 Ag 10.11255 9.17476 3.08179 0 0 0 Ag 10.11255 14.17917 3.08179 0 0 0 Ag 13.00185 4.17035 3.08179 0 0 0 Ag 13.00185 9.17476 3.08179 0 0 0 Ag 13.00185 14.17917 3.08179 0 0 0 Ag 0.00000 1.66814 3.08179 0 0 0 Ag 0.00000 6.67255 3.08179 0 0 0 Ag 0.00000 11.67697 3.08179 0 0 0 Ag 2.88930 1.66814 3.08179 0 0 0 Ag 2.88930 6.67255 3.08179 0 0 0 Ag 2.88930 11.67697 3.08179 0 0 0 Ag 5.77860 1.66814 3.08179 0 0 0 Ag 5.77860 6.67255 3.08179 0 0 0 Ag 5.77860 11.67697 3.08179 0 0 0 Ag 8.66790 1.66814 3.08179 0 0 0 Ag 8.66790 6.67255 3.08179 0 0 0 132 Ag 8.66790 11.67697 3.08179 0 0 0 Ag 11.55720 1.66814 3.08179 0 0 0 Ag 11.55720 6.67255 3.08179 0 0 0 Ag 11.55720 11.67697 3.08179 0 0 0 K_POINTS {automatic} 2 2 2 1 1 1 References: 1 P. Hohenberg and W. Kohn, ?Inhomogeneous Electron Gas?, Phys. Rev. B 136, 864-871 (1964) 2 W. Kohn and L.J. Sham, ?Self-Consistent Equations including Exchange and Correlation Effects?, Phys. Rev. 140, A1133 (1965) 3 M.R.C. hunt, S. Modesti, P. Rudolf and R.E. Palmer, ?Charge Transfer and Structure in C60 Adsorption on Metal Surfaces?, Phys. Rev. B 51,10039 (1995) 4 L. Wang and H. Cheng, ?Density Functional Study of the Adsorption of A C60 Monolayer on Ag(111) and Au(111) Surfaces?, Phys. Rev. B 69, 165417 (2004) 5 H. Mizuseki, N. Igarashi, R.V. Belosludov, A.A. Farajian, and Y. Kawazoe, ?An Organic Molecule-Fullerence Mixture for A High Efficiency Photovoltaic Device: A Theoretical Study?, Nanotech, 2, 72-75 (2003) 6 T. Toccoli, A. Boschetti and S. Iannotta, ?Molecular Materials for Optoelectronics by Supersonic Molecular Beam Growth: Co-deposition of C60 and ZnPc?, Synthetic Metals, 122-1, 229 (2001) 7 H.D. Burrows and A.A. kharlamov, ?About Energy and Electron Transfer Processes in C60/Phthalocyanine Films?, Progr Colloid Polym Sci, 123, 52 (2004) 133 8 L. A. Girifalco, ?Molecular Properties of C60 in the Gas and Solid Phases?, J. Phys. Chem. 96, 858-861 (1992) 9 L.A. Girifalco, ?Interaction Potential for C60 Molecules?, J. Phys. Chem. 95, 5370-5371 (1991) 10 J.-N. Chazalviel, Coulomb Screening by Mobile Charges: Applications to Materials Science, Chemistry, and Biology, (Birkh?user, 1998) 11 T. Toccoli, A. Boschetti and S. Iannotta, ?Molecular Materials for Optoelectronics by Supersonic Molecular Beam Growth: Co-deposition of C60 and ZnPc?, Synthetic Metals 122-1, 229 (2001) 12 A. Calzolari, W. Jin, J. E. Reutt-Robey, and M. B. Nardelli, ?Substrate-Mediated Intermolecular Hybridization in Binary Phthalocyanine Superstructures?, J. Phys. Chem. C 114, 1041 (2010) 13 W. Tang, E. Sanville, and G. Henkelman, ?A Grid-based Bader Analysis Algorithm Without Lattice Bias?, J. Phys.: Condens. Matter 21, 084204 (2009) 14 J.G. Che and H.-P. Cheng, ?First-principles Investigation of A Monolayer of C60 on h-BN/Ni(111)?, Phys. Rev. B 72, 115436 (2005) 15 M. Muntwiler, W. Auwarter, A.P. Seitsonen, J. Osterwalder, and T. Greber, ?Rocking-motion-induced Charging of C60 on h-BN/Ni(111)?, Phys. Rev. B 71, 121402(R) (2005) 134 Chapter 6 Conclusions Organic molecules can form self-assembled monolayers on metal and organic surfaces. These self-assembly processes are governed by intermolecular and substrate interactions, including long ranged charge-charge Coulomb interactions and short ranged vdW, dipole-dipole interactions and hydrogen bond formation. Although vdW interactions, dipole-dipole interactions and hydrogen bonds are weaker than covalent chemical bonds, they still play a key role in the molecular assembly process. In some cases, there could be large charge transfer between molecules, and this charge transfer could significantly change the surface morphology. This in turn can affect the electronic properties, which may have important applications in electronic industry. Because of the complicity of the interactions, the organic molecular monolayer can form various patterns on the surface, both regular ordered structures and irregular patterns. These patterns are not only important for their intrinsic properties, e.g. electronic property, and mass transportation, but also give us a key to understand the underling physics and chemistry in the system. The study of the structures and dynamics of these organic monolayers has been the main focus of this thesis. We use different modeling methods to study different system. For the regular ACA/Ag(111) system, we use multi-component lattice gas models, which we analyzed using both computer simulations and analytic methods like mean field theory and the pair approximation. This lattice gas model is simple, but can explain 135 the basic competition between intermolecular and substrate interactions, which we believe is an essential part of pattern formation on organic monolayers. This model can be modified with different parameters and can be used to describe various monolayer systems: isotropic and anisotropic, with and without sub-lattices, etc. For the irregular C60 patterns on ZnPc/Ag(111), we use an effective Girifalco potential to model the C60-C60 vdW interaction and use effective substrate potential to describe the substrate interaction. By using 2D Langevin dynamics simulation, we show that to form the single C60?s chain structures, there must be strong longer ranged repulsions between C60 molecules. We suggest these could come from dipole-dipole repulsions, or Coulomb repulsions between charged C60?s. This is a reasonable assumption because C60/ZnPc is good electron acceptor/donor system. To further investigate charge transfer in the C60/ZnPc/Ag(111) system, we used DFT. The results are preliminary but seem promising, and consistent with our basic assumption of large charge transfer in the C60/ZnPc/Ag(111) system. It seems plausible that this could modify the effective interactions between C60 molecules, as suggested in our 2D empirical models and the magnitudes are comparable. We are expecting a definitive analysis of the effective C60 interactions on organic overlayers which might be done with the development of a new generation of exchange-correlation functionals in DFT which properly describe both van der Waals and charge-charge interactions. 136 List of References 1. H.-C. Jeong and Ellen D. Williams, ?Step on Surfaces: Experiment and Theory?, Surface Science Report, 34, 171-294 (1999) 2. M.P.A. Fisher, D.S. Fisher and J.D. Weeks, ?Agreement of Capillary-Wave Theory with Exact Results for the Interface Profile of the Two-Dimensional Ising Model?, Phys. Rev. Lett. 48, 368 (1982) 3. N. C. Bartelt, J.L. Goldberg, T. L. Einstein, E.D. Williams, J. C. Heyraud and J.J. Metois, ?Brownian Motion of Steps on Si(111)?, Phys. Rev. B 48, 15453 (1993) 4. C. Dasgupta, M. Constantin, S. Das Sarma, and S. N. Majumdar, ?Survival in Equilibrium Step Fluctuations?, Phys. Rev. E 69, 022101 (2004) 5. S. Kodiyalam, K.E. Khor, N.C. Bartelt, E.D. Williams and S. Das Sarma, ?Energetics of Vicinal Si(111) Steps Using Empirical Potentials?, Phys. Rev. B 51, 5200 (1995) 6. O. Bondarchuk, D. B. Dougherty, M. Degawa, E.D. Williams, M. Constantin, C. Dasgupta, and S. Das Sarma, ?Correlation Time for Step Structural Fluctuations?? Phys. Rev. B 71, 045426 (2005) 7. M. Wortis, in Chemistry and Physics of Solid Surfaces VII, eds. R. Vanselow and R.F. Howe, 367-428 (Springer-Verlag, Berlin, 1988) 8. B. A. Gregg, S-G. Chen, and R. A. Cormier, ?Coulomb Forces and Doping in Organic Semiconductors?, Chem. Mater. 16 (23), 4586-4599 (2004) 9. T. Taima, J. Sakai, T. Yamanari, and K. Saito, ?Doping effects for organic photovoltaic cells based on small-molecular-weight semiconductors?, Sol. Energy Mater. Sol. Cells, 93, 6-7, 742 (2009) 10. D. Bonifazi; A. Kiebebe, M. Stohr, F.Y. Cheng, T. Jung, F. Diederich, and H. Spillmann, ?Supramolecular Nanostructuring of Silver Surfaces via 137 Self-Assembly of Fullerene and Porphyrin Modules?, Advanced Functional Materials, 17(7), 1051-1062 (2007) 11. M. A. Baldo and S. R. Forrest, ?Interface-limited Injection in Amorphous Organic Semiconductors?, Phys. Rev. B 64, 085201 (2001) 12. B. Xu, C.G. Tao, W.G. Cullen, J.E. Reutt-Robey and E.D. Williams, ?Chiral Symmetry Breaking in Two-dimensional C60-ACA Intermixed Systems?, Nano Lett. 5, 2207-2211 (2005) 13. B. Xu, C.G. Tao, E.D. Williams and J.E. Reutt-Robey ?Coverage Dependent Supramolecular Structures: C60:ACA monolayers on Ag(111)?, J. Am. Chem. Soc. 128, 8493-8499 (2006) 14. B. Xu, B. Varughese, D. Evans and J. Reutt-Robey ?Morphology Selected Molecular Architecture: Acridine Carboxylic Acid Monolayers on Ag(111)?. J. Phys. Chem. B 110, 1271-1276 (2006) 15. G.A. Jeffrey, An Introduction to Hydrogen Bonding (Oxford University Press, 1997) 16. S. Tsuzuki, K. Honda, T. Uchimaru, M. Mikami and K. Tanabe, ?Origin of Attraction and Directionality of the Pi/Pi Interaction: Model Chemistry Calculations of Benzene Dimer Interaction?, J. Am. Chem. Soc. 124, 104, (2002) 17. E.D. Williams, R.J. Phaneuf, J. Wei, N.C. Bartelt and T.L. Einstein, ?Thermodynamics and Statistical Mechanics of the Faceting of Stepped Si(111)?, Surface Science 294, 219-242 (1993) 18. M. Giesen-Seibert and H. Ibach, ?On the Time structure of Tunneling Images of Steps?, Surface Science 316, 205-222 (1994) 19. E. Le Goff, L. Barbier, and B. Salanon, ?Time?space Height Correlations of Thermally Fluctuating 2-d Systems: Application to Vicinal Surfaces and Analysis of STM Images of Cu(115)?, Surface Science 531, 337-358 (2003) 138 20. O. Bondarchuk, D.B. Dougherty, T.L. Einstein and E.D. Williams, ?Dynamics of Step Fluctuation on a Chemically Heterogeneous Surface of Al/Si(111)?, Phys. Rev. B 66, 085327 (2002) 21. A. Pimpinelli, J. Villain, D.E. Wolf, J.J. M?tois, J.C. Heyraud, I. Elkinani and G. Uimin, ?Equilibrium Step Dynamics on vicinal Surfaces?, Surface Science 295, 143-153 (1993) 22. M. Blume and V. J. Emery? ? Ising Model for the ? Transition and Phase Separation in He3-He4 Mixtures?? Phys. Rev. A 4, 1071 (1971) 23. N. Metropolis and S. Ulam, ?The Monte Carlo Method?, Journal of the American Statistical Association, 44, 335 (1949) 24. N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A.H. Teller and E. Teller, ?Equation of State Calculations by Fast Computing Machines?, J. Chem. Phys. 21, 1087 (1953) 25. M. Bader, J. Haase, K.-H. Frank and A. Puschmann, ?Orientational Phase Transition in the System Pyridine/Ag(111): A Near-edge X-ray-Absorption Fine-Structure Study?, Phys. Rev. Lett. 56, 1921 (1986) 26. A.L. Johnson and E.L. Muetterties ?Chemisorption Geometry of Pyridine on Pt(111) by NEXAFS?, J. of Phys. Chem. 89, 4071-4075 (1985) 27. J. Bonello and R.M. Lambert ?The Structure and Reactivity of Quinoline Overlayers and the Adsorption Geometry of Lepidine on Pt(111): Model Molecules for Chiral Modifiers in Enantioselective Hydrogenation?. Surface Science 498, 212-228 (2002) 28. J. Kubota and F. Zaera ?Adsorption Geometry of Modifiers as Key in Imparting Chirality to Platinum Catalysts?, J. Am. Chem. Soc. 123, 11115-11116 (2001) 29. E.D. Williams, O. Bondarchuk, C.G. Tao, W. Yan, W.G. Cullen, P.J. Rous and T. Bole, ?Temporal Step Fluctuations on A Conductor Surface: Electromigration Force, Surface Resistivity and Low-frequency Noise?, New Journal of Physics 9, 387 (2007) 139 30. B.A. Gregg, ?The Photoconversion Mechanism of Excitonic Solar Cells?, Materials Research Society Bulletin 30, 20-22 (2005) 31. R. Otero, D. Ecija, G. Fernandez, J. M. Gallego, L. Sanchez, N. Martin, and R. Miranda, ?An Organic Donor/Acceptor Lateral Superlattice at the Nanoscale?, Nano Lett. 7, 2602-2607 (2007) 32. J. Weckesser, J. V. Barth, C. Cai, B. Muller, and K. Kern, ?Binding and Ordering of Large Organic Molecules on and Anisotropic Metal Surfaces: PVBA on Pd (110)?, Surface Science 431, 168-173 (1999) 33. K. Braun and S. Hla, ?Probing the Conformation of Physisorbed Molecules at the Atomic Scale Using STM Manipulation?, Nano Lett. 5, 73-76 (2005) 34. R. Hiesgen, M. Rabisch, H. Bottcher, and D. Meissner, ?STM Investigation of the Growth Structure of Cu-phthalocyanine Films with Submolecular Resolution?, Solar Energy Materials & Solar Cells 61, 73-85 (2000) 35. C. Tao, Q. Liu, B. C. Riddick, W. G. Cullen, J. Reutt-Robey, J. D. Weeks, and E. D. Williams, ?Dynamic interfaces in an organic thin film?, Proc. Natl. Acad. Sci. USA, 105, 16418-16425 (2008) 36. H. A. Bethe, ?Statistical Theory of Superlattices?, Proc. R. Soc. Lond. A 150, 552-575 (1935) 37. S. Katsura and m. Takizawa, ?Bethe Lattice and the Bethe Approximation?, Progress of Theoeretical physics, 51, No. 1 (1974) 38. K. Huang, Statistical Mechanics, 2nd edition, 342-367 (Wiley, 1987) 39. J. D. Weeks and G. H. Gilmer, ?Pair Approximation Equations for Interfaces and Free Surfaces in the Ising Model?, J. Chem. Phys. 63, 3136 (1975) 40. R. Salazar and L. D. Gelb, ?Application of the Bethe-Peierls Approximation to A Lattice-gas Model of Adsorption on Mesoporous Materials?, Phys. Rev. E 71, 041502 (2005) 140 41. R. Kikuchi, ?Theory of Domain Walls in Ordered Structures II?, J. Phys. Chem. Solids, 23, 137-151 (1962) 42. J.-C. Lee, Thermal Physics - Entropy and Free Energies, chapter 5, (New Jersey: World Scientific, 2002) 43. R. Kikuchi, ?A Theory of Cooperative Phenomena?, Phys. Rev. 81, 988-1003 (1951) 44. C. Stadler, S. Hansen, I. Kroger, C. Kumpf and E. umbach, ?Tuning Intermolecular Interaction in Long-range-ordered Submonolayser Organic Films?, Nature Physics, 5, 153 (2009) 45. A. Calzolari, W. Jin, J. E. Reutt-Robey, and M. B. Nardelli, ?Substrate-mediated Intermolecular Hybridization in Binary Phthalocyanine Superstructures?, J. Phys. Chem. C, 114, 1041 (2010) 46. A. Z. Akheyan and N. S. Ananikian, ?Global Bethe Lattice Consideration of the Spin-1 Ising Model?, J. Phys. A: Math. Gen. 29, 721-731 (1996) 47. R. R. Levitskii, O. R. Baran, and B. M. Lisnii, ?Phase Diagrams of Spin-1 Ising Model with Bilinear and Quadrupolar Interactions under Magnetic Field, Two-particle Cluster Approximation?, Eur. Phys. J. B 50, 439-443 (2006). 48. C. Temirci, A. Kokce, and M. Keskin, ?Equilibrium Properties of a spin-1 Ising System With Bilinear, Biquadratic and Odd Interactions?, Physica A 231, 673-686 (1996) 49. M.R.C.Hunt, S.Modesti, P.Rudolf, R.E.Palmer ?Charge Transfer and Structure in C60 Adsorption on Metal Surfaces?, Phys. Rev. B 51, 10039 (1995) 50. L.Wang, H.Cheng, ?Rotation, Translation, Charge Transfer, and Electronic Structure of C60 on Cu(111) Surface?, Phys. Rev. B 69, 045404 (2004) 51. L.Wang, H.Cheng, ?Density Functional Study of the Adsorption of A C60 Monolayer on Ag(111) and Au(111) Surfaces?, Phys. Rev. B 69, 165417 (2004) 141 52. L.A.Girifalco "Molecular Properties of Fullerene in the Gas and Solid Phases", J. Phys. Chem. 96, 858-861 (1992) 53. Z.R.Hong, B.Maennig, R.Lessmann, M.Pfeiffer, and K.Leo "Improved Efficiency of Zinc Phthalocyanine/C60 Based Photovoltaic Cells via Nanoscale Interface Modification", Appl. Phys. Lett. 90, 203505 (2007) 54. H.D.Burrows, A.A. Kharlamov, ?About Energy and Electron Transfer Processes in C60/Phthalocyanine Films? Progr Colloid Polym Sci 123, 52 (2004) 55. W.Jin, Q.Liu, J.D. Weeks and J. Reutt-Robey, in preparation 56. H.L.Zhang et al., ?C60 Molecular Chains on ?-Sexithiophene Nanostripes? Small, 3, 2015-2018 (2007) 57. M. C. Abramo, C. Caccamo, D. Costa, G. Pellicane and R. Ruberto, ?Atomistic Versus Two-body Central Potential Models of C60: A Comparative Molecular Dynamics Study?, Phys. Rev. E 69, 031112 (2004) 58. J. Stambaugh, K. V. Workum, J. F. Douglas and W. Losert, ?Polymerization Transitions in Two-dimensional Systems of Dipolar Spheres?, Phys. Rev. E 72, 031301 (2005) 59. K. H. Lau and W. Kohn, ?Indirect Long-range Oscillatory Interaction between Adsorbed Atoms?, Surface Science 75, 69-85 (1978) 60. S. Koh and G. Ehrlich, ?Self-Assembly of One-Dimensional Surface Structures: Long-Range Interactions in the Growth of Ir and Pd on W(110)?, Phys. Rev. Lett. 87, 106103 (2001) 61. M.M. Hurley and S.J. Singer, ?Domain-array Melting in the Dipolar Lattice Gas? Phys. Rev. B 46 (9), 5783 (1992) 62. A.D. Stoycheva and S. J. Singer, ?Computer Simulations of A Two-dimensional System with Competing Interactions?, Phys. Rev. E 65, 036706 (2002) 142 63. M. Linares, D. Beljonne, J. Cornil, etc. ?On the Interface Dipole at the Pentacene-Fullerene Heterojunction: A Theoretical Study?, J. Phys. Chem. C 114, 3215-3224 (2010) 64. I. Avilov, V. Geskin and J. Cornil, ?Quantum-Chemical Characterizaiton of the Origin of Dipole Formation at Molecular Organic/Organic Interfaces?, Adv. Func. Mater. 19, 624-633 (2009) 65. W. Smith, C.W. Yong and P.M. Rodger, ?DL_POLY: Application to Molecular Simulation?, Molecular Simulation 28, 385 (2002) 66. H. Kamiyama, K. Ohno, and y. Kawazoe, ?Classical MD Simulation of C60 Adsorbed on GaAs(001) Surface?, Sci. Rep. RITU A41, 187-190 (1996) 67. S.Mayo, B.D.Olafson, W.A.Goddard, ?Dreiding: A Generic Force Field for Molecular Simulations?, J. Phys. Chem 94(26) 8897 (1990) 68. T. A. Halgren, ?Pepresentation of van der Waals (vdW) Interactions in Molecular Mechanics Force Fields: Potential Form, Combination Rules, and vdW Parameters?, J. Am. Chem. Soc. 114, 7827-7843 (1992) 69. S. Nose, ?A Unified Formulation of the Constant Temperature Molecular Dynamics Methods?, J. Chem. Phys. 81, 511 (1984) 70. H. Mizuseki, N. Igarashi, R.V. Belosludov, A.A. Farajian, Y. Kawazoe, ?An Organic Molecule Fullerene Mixture for A High Efficiency Photovoltaic Device: A Theoretical Study?, Nanotech, 2 (2003) 71. J.M. Rodgers and J.D. Weeks, ?Accurate Thermodynamics for Short-ranged Truncations of Coulomb Interactions in Site-site Molecular Models?, J. Chem. Phys. 131, 244108 (2009) 72. W.T. Coffey, Yu.P. Kalmykov, J.T. Waldron ?The Langevin Equation with Applications to Stochastic Problems in Physics, Chemistry and Electrical Engineering? World Scientific Series in Contemporary Chemical Physics, 14 (2004) 143 73. P. Hohenberg and W. Kohn, ?Inhomogeneous Electron Gas?, Phys. Rev. B 136, 864-871 (1964) 74. W. Kohn and L.J. Sham, ?Self-Consistent Equations including Exchange and Correlation Effects?, Phys. Rev. 140, A1133 (1965) 75. H. Mizuseki, N. Igarashi, R.V. Belosludov, A.A. Farajian, and Y. Kawazoe, ?An Organic Molecule-Fullerence Mixture for A High Efficiency Photovoltaic Device: A Theoretical Study?, Nanotech, 2, 72-75 (2003) 76. L. A. Girifalco, ?Molecular Properties of C60 in the Gas and Solid Phases?, J. Phys. Chem. 96, 858-861 (1992) 77. L.A. Girifalco, ?Interaction Potential for C60 Molecules?, J. Phys. Chem. 95, 5370-5371 (1991) 78. J.-N. Chazalviel, Coulomb Screening by Mobile Charges: Applications to Materials Science, Chemistry, and Biology, (Birkh?user, 1998) 79. W. Tang, E. Sanville, and G. Henkelman, ?A Grid-based Bader Analysis Algorithm Without Lattice Bias?, J. Phys.: Condens. Matter 21, 084204 (2009) 80. J.G. Che and H.-P. Cheng, ?First-principles Investigation of A Monolayer of C60 on h-BN/Ni(111)?, Phys. Rev. B 72, 115436 (2005) 81. M. Muntwiler, W. Auwarter, A.P. Seitsonen, J. Osterwalder, and T. Greber, ?Rocking-motion-induced Charging of C60 on h-BN/Ni(111)?, Phys. Rev. B 71, 121402(R) (2005)