Keywords: Monte Carlo simulations, Charge transport, Energetic fluctuations, Electro-Mechanical response, Acidity constant (pK a ), Polyanions Strain-Mobility trends in crystalline organic semiconductors

To a spirit, that drifts with the misty morning breeze, relishing the fragrance that prevails ... that can dance to the tunes of life, holding the beauty of the melodious unique self ... To a spirit, drenched by the mysteries of life, yet, standing tall, surmounting the challenges of time ... traversing through the mystic woods, yet, can embrace the unison with the nature that surrounds ... walking through the stormy winds, yet, retaining the self in the sanctum of will ... sailing through the oceans of unexplored self, yet, staying afloat against the tides of the perceptions that reigns ... To a spirit, that can evolve in time, challenging the experiential self, that can hold itself, through the dogmas of life, that can find its way, through the darkest of the thoughts, that can ignite the soul through the esteemed self, that can stand alone, in the seclusion of mind, that can relish the fire of life, in the solitude of a self, that is driven by the will to explore the sanities of life, TO THAT SPIRIT, within us, I dedicate this work ... v

R ÉSUM É

Les avancées technologiques et l'intégration massive de dispositifs électroniques nanométriques dans les objets de notre vie quotidienne ont généré une explosion des coûts de R & D, de conception et de production, ainsi que des inquiétudes sociétales quant à l'impact environnemental des déchets électroniques. En raison de procédés de production moins coûteux et à faible impact environnemental, de leur souplesse d'utilisation et de la possibilité de moduler leurs propriétés à l'infini, les molécules et polymères organiques constituent une classe de matériaux prometteuse pour la mise au point de nouveaux dispositifs électroniques. L'électronique organique couvre ainsi un vaste domaine d'applications, parmi lesquelles se trouvent les diodes électroluminescentes, les transistors à effet de champ ou les cellules photovoltaiques. Bien que certains de ces dispositifs soient déjà commercialisés, les processus gouvernant leur efficacité à l'échelle atomique sont loin d'être entièrement compris et maitrisés. C'est en particulier le cas des processus de transport de charge, qui interviennent dans tous ces dispositifs.

L'objectif de cette thèse est d'apporter une compréhension fondamentale des processus de transport de charge dans les semiconducteurs organiques, à partir d'approches théoriques combinant dynamique moléculaire, calculs quantiques et simulations Monte Carlo. Ce travail est développé suivant trois axes principaux:

[I] Étude des relations liant l'organisation structurale et les propriétés de transport de cristaux moléculaires, et du rôle des fluctuations énergétiques dans des matériaux polymères amorphes. Des simulations Monte Carlo Cinétique (KMC) couplés au formalisme de Marcus-Levich-Jortner pour le calcul des taux de transfert ont été effectués afin de déterminer les mobilités des électrons et des trous au sein de dix structures cristallines de dérivés phtalocyanines. Dans une deuxième étude, une approche similaire a été employée afin de décrire les propriétés de transport de charge au sein d'un copolymère amorphe de fluorène-triphénylamine, ainsi que l'impact des fluctuations énergétiques sur ces dernières. La méthodologie développée permet d'obtenir, pour un i faible coût calculatoire, une estimation semi-quantitative des mobilités des porteurs de charge dans ce type de système.

[II] Étude de l'impact de contraintes mécaniques sur les propriétés de transport de matériaux organiques cristallins. La réponse électronique et les propriétés de transport de matériaux organiques soumis à une contrainte mécanique ont été étudiés à l'aide de simulations de dynamique moléculaire et de calculs DFT. Le rubrène cristallin et ses polymorphes, ainsi que les dérivés du BTBT, ont été considérés pour cette étude, qui révèle un couplage électromécanique inhabituel entre les différents axes cristallographiques. Les résultats démontrent en particulier que l'anisotropie structurale des monocristaux organiques conduit à une anisotropie du couplage électromécanique.

[III] Étude du rôle du polyélectrolyte dans la conductivité des complexes conducteurs. Le polystyrène substitué par du bis(sulfonyl)imide est utilisé comme un contre-ion et un dopant dans les complexes conducteurs PEDOT-polyélectrolytes. En complément des analyses expérimentales, des simulations de dynamique moléculaire couplés à des calculs DFT ont été effectués dans ces systèmes afin d'analyser l'impact de la conformation et de l'état de protonation du polyélectrolyte sur la conductivité du complexe formé avec le PEDOT.

Les études décrites ci-dessus, réalisées sur différents types de matériaux en couplant différents types d'approches théoriques, ont permis d'apporter une compréhension fondamentale des propriétés de transport dans les semiconducteurs organiques. Elles mettent en particulier en évidence l'impact de l'organisation structurale, des interactions intermoléculaires et de l'application de contraintes mécaniques sur la mobilité des porteurs de charges dans ces matériaux.

Mots-clés : Simulations Monte Carlo, transport de charge, désordre énergétique, réponse électro-mécanique, constante d'acidité, polyanions.

ABSTRACT

With the advancement of technology, miniaturized electronic devices are progressively integrated into our everyday lives, generating concerns about cost, efficiency and environmental impact of electronic waste. Organic electronics offers a tangible solution paving the way for low-cost, flexible, transparent and environment friendly devices. However, improving the functionalities of organic (opto) electronic devices such as light emitting diodes and photovoltaics still poses technological challenges due to factors like low efficiencies, performance stability, flexibility etc. Although more and more organic materials are being developed to meet these challenges, one of the fundamental concerns still arises from the lack of established protocols that correlate the inherent properties of organic materials like the chemical structure, molecular conformation, supra-molecular arrangement to their resulting charge-transport characteristics. In this context, this thesis addresses the prediction of charge transport properties of organic semiconductors through theoretical and computational studies at the atomistic scale, developed along three main axes:

[I] Structure-charge transport relationships of crystalline organic materials and the role of energetic fluctuations in amorphous polymeric organic semiconductors. Kinetic Monte-Carlo (KMC) studies employing the Marcus-Levich-Jortner rate formalism are performed on ten crystalline Group IV phthalocyanine derivatives and trends linking the crystalline arrangement to the anisotropic mobility of electrons and holes are obtained. Subsequently, KMC simulations based on the simpler Marcus formalism are performed on an amorphous semiconducting fluorene-triphenylamine (TFB) copolymer, to highlight the impact of energetic fluctuations on charge transport characteristics. A methodology is proposed to include these fluctuations towards providing a semi-quantitative estimate of charge-carrier mobilities at reduced computational cost.

[II] Impact of a mechanical strain on the electronic and charge transport properties iii of crystalline organic materials. Crystalline rubrene and its polymorphs, as well as BTBT derivatives (well studied high mobility organic materials) are subjected to mechanical strain and their electronic response is analyzed. Employing tools like Molecular Dynamic (MD) simulations and plane wave DFT (PW-DFT) calculations, unusual electro-mechanical coupling between different crystallographic axes is demonstrated, highlighting the role of inherent anisotropy that is present in the organic single crystals which translates in an anisotropy of its electro-mechanical coupling.

[III] Protonation-dependent conformation of polyelectrolyte and its role in governing the conductivity of polymeric conducting complexes. Polymeric bis(sulfonyl)imide substituted polystyrenes are currently employed as counter-ions and dopants for conducting poly (3,4-ethylenedioxythiophene) (PEDOT), resulting in PEDOT:polyelectrolyte conducting complexes. Employing MD simulations and DFT calculations, inherent characteristics of the polyelectrolyte like its acid-base behavior, protonation state and conformation, are analyzed in conjunction with available experimental data and the role of these characteristics in modulating the conductivity of resulting PEDOT:polyelectrolyte conducting complexes is highlighted.

The above studies, performed on different organic electronic systems, emphasize the importance of inherent characteristics of organic materials in governing the charge transport behavior in these materials. By considering the inherent characteristics of organic electronic materials and systematically incorporating them into simulation models, accuracy of simulation predictions can be greatly improved, thereby serving not only as a tool to design new, stable and high performance organic materials but also for optimizing device performances.

• From the kritis (endarō mahānubhāvulu, Ghana Raga Pancharatnam) of Kakarla Tyagabrahmam # (Tyāgayya, Composer of Carnatic music, 1767-1847)

"If I have seen further it is by standing on the shoulders of Giants..."

• Issac Newton * (mathematician, astronomer and physicist, 1642-1726) Time, space and area of expertise would separate Kakarla Tyagabrahmam # (Tyāgayya, who is considered as one of the greatest composers of Carnatic music -Indian classical music) and Issac Newton * (one of the greatest scientist and theoretician of the modern era), yet, the two quotes are connected and holds truth in my case as this dissertation would have been a behemoth task to accomplish without the collaborative contributions from many people from different areas of expertise, be in computational chemistry or experimental. I consider it to be an honor to thank the individuals that, either directly or indirectly, provided insights, suggestions, inputs and/or collaboratively contributed to my work, which is presented in this dissertation.

First and foremost, I offer my sincere thanks to Dr Luca MUCCIOLI and Prof.

Frédéric CASTET, for their enormous support and encouragement in guiding me through this PhD work. As a graduate student, traversing through different fields from mechanical engineering to composite polymers and then on-to material physics and material chemistry, I was a novice to the field of computational chemistry with very little knowledge to understand the complex interplay in organic electronic systems. However, Dr. Luca MUCCIOLI, had been a complete mentor, nurturing me to develop, not only my scientific and computational temperament during the thesis vi progression, but also to build my scientific persona as a whole, and also in developing fool proof scientific arguments. He had also been a very good friend, philosopher and a councilor, offering me guidance both at scientific and personal fronts. I offer my sincere thanks and gratefulness to Dr. MUCCIOLI. My sincere thanks also to Prof. Frédéric CASTET, for his guidance and detail-focused suggestions during the numerous scientific discussions we had over time, that I was able to corroborate my results in a detailed and constructive manner. It is with the direction of both Dr. MUCCIOLI and Prof. CASTET, that I was able to explore the field of computational chemistry applied to organic electronics, and in the process understand the intricacies of associating the computational models to experimental observations, for which I am extremely grateful to them. Their mentor-ship helped me not only to gain insights from theoretical results but also developing an ability to communicate with experimentalists, as this area of computational chemistry could become more fruitful with a symbiotic relationship between experimental observations and computational results, so as to propel the field of organic electronics towards sustainable and efficient developments. Additionally I would like to also thank Dr. Lionel.A. TRUFLANDIER, Dr Gabriele D'AVINO and Dr. Yoann OLIVIER for many scientific discussions and their extensive guidance in developing my computational skills.

I would like to thank Dr. Micaela MATTA and Dr Claire TONNEL É, with whom I collaborated extensively during my thesis work, whose contributions were not only significant towards the progress of this thesis work, but also very fruitful for scientific accomplishments associated with this thesis. Within the same context, I thank also Dr. Guillaume WANTZ, Dr. Cedric AYELA, Mr. Marco PERIERA, (University of Bordeaux) Dr Anna.I.HOFMANN (former PhD student, University of Bordeaux) and Dr Benoît LESSARD (University of Ottawa) for their extensive experimental collaborations and insights that were essential both for the thesis work and also for scientific contributions. Additionally, during the development of this thesis work, I had the opportunity of visiting the group of Prof. Roberto LAZZORONI at the University of Mons-Laboratory for Chemistry of Novel Materials, and in the process gain some theoretical insights under the guidance of Dr. David BELJONNE and Dr. Jérôme CORNIL and start a collaborative work with Dr. Claudio QUARTI, Dr. Julien ID É and Dr. Vincent LEMAUR, which was also essential in developing this dissertation, for which I offer my sincere thanks to all of them. I would also like to thank Prof. Benoît CHAMPAGNE and his group at University of Namur, for their hospitality in accommodating me as a scientific and/or personal visitor.

I also thank Prof. Georges HADZIIOANNOU for introducing me to the field of organic electronics, and hosting me as a PhD student in his group for a considerable amount of time. I am also grateful to Ms. Aude MANSON for all the administrative assistance during the initial days of my PhD work in Bordeaux. I would also like to take this opportunity to thank all the scientific personnel of the group of Prof. HADZIIOANNOUI for fruitful collaborations and scientific discussions. Further, I also thank Dr. Pascal LARR ÉGARAY and all the personnel of Institut des Sciences Moléculaires (ISM)-THEO, for their hospitality and a friendly working environment. vii Additionally, I offer my sincere gratefulness to Mr. Philippe AUREL at ISM-THEO for his help in developing the essential scripts that were extremely helpful to accelerate my calculations. I am also grateful to LABEX-AMADEus for the financial support for my doctoral work, with special thanks to Prof. Etienne DUGUET, Prof. Corine MATHONIERE and Dr. Olivier FOUASSIER for their help and support during my scientific endeavors at University of Bordeaux.

I take this opportunity to thank all my friends in Bordeaux, Namur and Mons, names of whom I cannot take due to space constraints (for which I offer my apologies to them), and thanks to the multi-lingual and multi-cultural environment that was provided through their friendship, that I was able to gain a life-time experience with significant lessons, that would stay with me for the rest of my life. Also, I am grateful to all of them for making me feel at home and accepting me to be a part of their lives.

Of the friends I was able to make an acquaintance with, three people stand-out from the rest, and I take the privilege of expressing my gratefulness to my closest friends -Dr. Iionnis PETSAGKOURAKIS, Ms. Kornelia PIELAK, Mr. Bruno FLAMENT, whose encouragement and support served as a constant source of energy during my thesis work. I find myself lost in searching for words, to express my gratitude to the "trio" as they were my mentors during times of crisis, my flag-posts when I lost sight of my path and also my motivators when I lost the energy to fight. I thank the "trio" for being a part of my life and whom I consider a part of my family.

Finally, I take the opportunity to thank my family, my parents Madhu and Usha, my brother Venki and my uncle (maternal) Nagendra for their enormous support and encouragement during the whole of my academic endeavors. I further thank and express my sincere gratefulness to my parents with all my heart, whose mentor-ship and brought-up had cultivated in me a sense of independence, for teaching me the value of self-esteem and nurturing in me the ability to never-surrender and fight for my beliefs and for the cause. Above all, thanks to each and everyone who crossed my life, thanks to the time I spent with them and experiences I gathered, to those who stayed and also to those who left, that I learned valuable life lessons and gained sufficient strength in the process, to face the world that awaits, for there is nothing to loose or to gain, but only an opportunity to learn does the life presents.. and with it an option to grow and move forward....

Cheers and Thanks,

Sai-M, Gali.

viii CONTENTS

INTRODUCTION

Miniaturization of electronic devices has been at the forefront of many technological developments that were witnessed in the recent decades. Semiconductor technology, serving as a building block for many advancements in the field of electronic devices, has evolved (almost) in-line with Moore's predictions [1,2] that the number of transistors on integrated circuits doubles every two years. Although much of this progress is based on inorganic semiconductors based on Si, Ge, GaAs, etc., inorganic electronic devices present some challenges. Fabrication of inorganic semiconductor materials is often a multi-step process demanding huge amount ot energy and involving several toxic chemicals, posing several challenges in terms of human exposure to chemicals, disposal of chemical waste, waste management etc., restricting materials like Si or Ge to be far from being eco-friendly for long term mass production [3,4,5]. Organic semiconductors offers tangible solutions, constituting an alternative gateway towards miniaturized micro-electronics or large-area flexible electronics such as flexible conducting/semiconducting substrates. Low-cost, low-temperature manufacturing and often high volume and high-throughput production, offer, in principle, an uncanny edge to organic electronic devices in comparison to their semi-conducting counterparts [6].

Although preliminary investigations on the electrical/photo conductivity of organic materials dates back to the early 1900s [7,8], limitations like reproducibility of carrier mobility, poor control of material purity, structural imperfections, efficiency, solubility, processability etc., [9,10] restricted the potential use of organic electronic devices, until 1980s where the potential applicability of organic materials was demonstrated by successfully incorporating them into Organic Photo-Voltaic Cells (OPVCs) [11] and Organic Field Effect Transistors (OTFTs) [12,13], although the earliest reported OTFTs showed low efficiencies. Subsequently, Organic Light Emitting Diodes (OLEDs) [14,15] and OTFTs with improved efficiencies were reported 1 [16,17,18,19,20], propelling the field towards new scientific endeavors in terms of new materials for improved performance, increased efficiencies, stability over time and ease of synthesis.

Materials and Devices

Organic materials, crystalline or amorphous, are constituted of hydrocarbon-based molecules with the carbon atoms defining the molecular backbone. These materials can be classified as semiconductors, conductors and dielectrics [6], depending on the electronic/electric response and the associated conductivity. The conductivity depends on mobility of charge carriers (electrons and hole) in the material and density of charge carriers, which in-turn depend on availability of the loosely bound electrons that are delocalized along the molecule backbone and their degree of delocalization. Organic semiconductors, for example, exhibit sp 2 hybridization of atomic orbitals between the adjacent carbon atoms with alternating single and double bonds. This results in bonding and antibonding σ molecular orbitals with a strong overlap, that form the backbone of the molecule, whereas the π molecular orbitals overlap to a lesser degree resulting in spatial delocalization of electrons in the π and π * orbitals. These delocalized electrons are not bound to a specific carbon atom (unlike the electrons in σ orbitals) but belong to the whole π-segment of the molecule and are responsible for the macroscopic conductivity of the molecular assembly. The degree of delocalization depends on the specific molecule type, dimension, bond length alternation pattern and conformation of the molecule. The electric/electronic activity, arises from the conjugated network of the molecules that are either chemically or physically bound. It is widely accepted in the scientific community that, in chemically bound materials like polymers, the delocalization of the electrons is spread along the aromatic backbone whereas in physically bound materials like molecular crystals or amorphous polymers the delocalization depends on the spatial arrangement of the molecular subunits with respect to each other. Over the recent decades, several small molecules, oligomers and polymers were developed, with physical and chemical modifications, to meet specific design requirements and applicability [21,22,23,24,25], to reach efficiencies (sometimes) on par with the inorganic counterparts.

Further, in comparison to their conventional inorganic counterparts, organic electronic devices are comprised of several stacks of thin film layers (such as conducting, semiconducting and insulating parts), with each layer designed for a specific purpose, and the overall device performance depends on processing conditions of each layer as well as that of the device, and the interplay of each layer with the others [26]. For the sake of illustration, a schematic representation of OPV and OLED devices containing several functional layers, along with their basic working principle, is presented in Fig. 1.1.

Modeling: State of the art

Organic electronic materials are ofter described as soft materials, as molecules in crystalline or amorphous organic systems are bound by weak van der Waals interactions. Addition of electron or hole to any molecule in the system, during the charge transport phenomena or injection from electrodes, induces a significant amount of intra and inter molecular deformations, that are both physical and chemical in nature.

Charge carrier (electron or hole), considering a semi-classical picture, is then drifting amidst these deformations in the quantum space constituted by the interacting forces. This quantum space is in-turn is affected by the presence of charge carrier, movement of the charge and the associated deformations. This charge carrier, thus acts as a quasi-particle in the dynamic quantum space and is commonly referred to as polaron [29,30]. As the interactions involved are relatively weak, several factors such as energetic coupling of neighboring molecules that act as donors and acceptors of charge carriers, stacking, packing and relative orientation of these coupled neighbors, effect of molecular deformations and change in polarization of the surrounding medium due to presence of charge, vibrational modes etc., influence the charge transport with in the material. In addition to this, as elicited in the previous section and schematically represented in Fig. 1.1, almost all organic electronic devices are made of layers of organic materials, with the charge carrier exiting from one layer and entering into another layer through organic-organic or organic-inorganic interface. As a consequence, at the device scale, both the intra-material charge transport and the inter-layer charge transport very important to determine the overall efficiencies of organic electronic devices.

With the complex interplay of intricate effects (energetic or morphological), combined with the subtlety and sensitivity of these interactions, modeling charge transport in organic electronics, be it semiconducting, conducting or dielectric materials and even more the electronic devices, remains a challenging task, even after decades of significant research in the field. Several pioneering works combining modeling methodologies at different length scales helped to understand the charge transport mechanism, be it hopping-like charge transport or band like. The influence of different characteristic properties like energetic disorder [31,29,32,33,34,35], carrier density [36,37] inter and intra molecular factors like reorganization [38,39] and polarization energies [40,41,42,43,44], lattice vibrations [45,46], electron-phonon coupling and charge localization [47,48,49], molecular conformations, interactions and aggregations [50,51,52,53,54], fluctuations of energetic disorder [55,46], effect of percolation network [56,57,58,59] etc., within the semi-classical hopping formalisms [38,60,61,62] were discussed. Further, external factors like imperfections [63,64,65], interfacial effects [66,67,68,69,70,71], mechanical strain [START_REF] Wu | Strain effects on the work function of an organic semiconductor[END_REF][START_REF] Yao | Charge Transport in Organic and Polymeric Semiconductors for Flexible and Stretchable Devices[END_REF][START_REF] Zheng | Understanding lattice strain-controlled charge transport in organic semiconductors: A computational study[END_REF] etc., also influence the charge carrier mobilities. Several interconnected and intricate factors, therefore, govern the charge transport in organic materials, both at the material scale and at the device scale.

Challenges & Motivation

This dissertation is addressed at correlating the physical and chemical properties of organic materials like the chemical structure, molecular conformation, supra-molecular arrangement to their resulting charge-transport characteristics. This thesis can be split into four main parts wherein the prediction of charge transport properties of organic semiconductors through theoretical and computational studies at the atomistic scale is addressed.

Part-1: Theoretical methodology that is related to this dissertation is presented in Chapter 2. Kinetic Monte Carlo (KMC) simulations are performed based on a home built FORTRAN code to obtain the charge transport characteristics and some benchmark calculations validating the KMC code are discussed in Chapter 3.

Part-2: Structure-charge transport relationships of crystalline organic materials and the role of energetic fluctuations in amorphous polymeric organic semiconductors.

In Chapter-4, Kinetic Monte-Carlo (KMC) studies employing the Marcus-Levich-Jortner rate formalism are performed on ten crystalline Group IV phthalocyanine derivatives, which vary in terms of metal/metalloid core, type/position/nature of function groups on the phthalocyanine unit. Role of energetic disorder on charge carrier mobility is discussed and the link between the crystalline arrangement of the compounds and the mobility anisotropy is obtained. In Chapter-5, KMC simulations based on the simpler Marcus formalism are performed on an amorphous semiconducting fluorene-triphenylamine (TFB) copolymer, to highlight the impact of energetic fluctuations on charge transport characteristics. A methodology is proposed to include these fluctuations towards providing a semi-quantitative estimate of charge-carrier mobilities, with respect to available experimental observations, at reduced computational cost. These results are presented in Chapter-5.

Part-3: Impact of a mechanical strain on the electronic and charge transport properties of crystalline organic materials.

Crystalline rubrene and its polymorphs, as well as BTBT derivatives (well studied high mobility organic materials) are subjected to mechanical strain and their electronic response is analyzed. Employing tools like Molecular Dynamic simulations and plane wave DFT calculations, unusual electro-mechanical coupling between different crystallographic axes is demonstrated, highlighting the role of inherent anisotropy that is present in the organic single crystals which translates in an anisotropy of its electro-mechanical coupling. These results are discussed in Chapters 6 and 7.

Part-4: Protonation-dependent conformation of polyelectrolyte and its role in governing the conductivity of polymeric conducting complexes.

Polymeric bis(sulfonyl)imide polystyrenes are currently employed as counter-ions and dopants for conducting poly (3,4-ethylenedioxythiophene), PEDOT, resulting in PEDOTpolyelectrolyte conducting complexes. Employing MD simulations and DFT calculations, inherent characteristics of the polyelectrolyte like its acid-base behavior, protonation state and conformation, are analyzed in conjunction with available experimental data and the role of these characteristics in modulating the conductivity of PEDOT-polyelectrolyte complexes is highlighted. Results are discussed in Chapter-8.

CHAPTER 2

THEORETICAL BACKGROUND

Theoretical details that are essential and used in the dissertation are discussed in the subsequent sections. Details about DFT calculations at molecular & periodic levels, molecular mechanics (MM) & molecular dynamics (MD), charge transport models and Kinetic Monte-Carlo methods are presented.

Charge transport : Methodology & Models

Charge-transport in organic electronic systems can be described by two limiting cases, either by hopping or band transport models. It is generally agreed in the scientific community that crystalline materials at low temperature exhibit band transport regime while at high temperature they exhibit hopping transport. Amorphous materials in general exhibit hopping transport. Several theoretical models were proposed to describe the charge carrier mobility in organic semiconductors, either it is hopping transport or band transport. These models under appropriate conditions provide a semi-quantitative estimate and capture a realistic picture of organic semiconducting system of interest. The conditions for the applicability of the model and approximations used in these models are detailed in the subsequent sections.

Hopping picture of charge transport

In organic semiconducting system -crystalline (or amorphous), intermolecular coupling between valence bands (energy levels) and/or conduction bands is weak at room temperature [1,2]. Consequently charge carriers -holes or electrons are localized (in the adiabatic limit) and travel, in-between subsequent phonon scattering events, a distance that equals lattice constant. This picture is contrary to the mechanism observed in inorganic semiconductors where charge-carrier is de-localized and exhibit a band like transport. Therefore charge transport in organic semiconductors is described by transfer of charge carrier among adjacent molecular sites that are electronically coupled. These coupled sites are called hopping sites. The charge transport between these sites is associated by a hopping rate that is proportional to the strength of electronic interaction between the sites. Hopping rate not only depends on the electronic coupling between hopping sites, but also on temperature and external electric field, which acts as stimuli in modulating the activation energy associated with every hopping event. Series of hopping events lead to charge propagation under the influence of electric field and temperature. This charge propagation is associated with charge drifting or diffusion in the materials until the charges reach electrodes at device scale or reach equilibrium in simulation models, which is then quantified with mobility, obtained as per Eq. 2.1.

µ = eD k B T (2.1)
where, D is the charge diffusion coefficient. Mobile charge carriers -electrons or holes, in semi-conducting systems are responsible for observed conductivity. Mobility, an intensive property, can be related to conductivity (σ κ ), elementary charge (e) and charge density (n), through Eq. 2.2,

µ = neσ κ (2.2)
Considering the complete picture, the mobility, within a good approximation, can be expressed as originating from two main contributions [3], given by Eq. 2.3,

µ T otal = µ T un + µ Hop (2.3)
At low temperatures a coherent transport is prevalent, wherein charges can tunnel (µ T un ) over long distances as phonon scattering is low, while increase in temperature leads to a non coherent transport described by hopping (µ Hop ). The relative contribution of each term is associated to the microscopic parameters such as electron-phonon coupling, electronic and phonon band widths and phonon energy. Distinction can be made between these two mechanisms (µ T un & µ Hop ) based on the strength of interaction of local (intramolecular) electron-phonon coupling, quantified by coupling constant g 2 j which in turn is related to the Huang-Rhys factor for the specific vibration mode j, S j by g 2 j = S j [3]. When weak electron-phonon coupling (g 2 P << 1) prevails, tunneling dominates and charge carrier mobility is band-like. However, in organic semiconductors is it not often true as a strong-electron phonon coupling (g 2 P ≥ 1) is often observed. In such conditions, (i) at extremely low temperatures a band-like transport can still prevail, (ii) increase in temperature activates hopping-like transport where the charge carriers are localized and are diffusing in the semiconductor by a thermally/electrically activated process and (iii) further increase in temperature or electric field leads to an increase in phonon scattering which forces the charge carriers to drift leading to a decrease in mobility. Representation of different regimes is presented in Fig. 2.1, to give a bird's eye view of the description provided.

Hopping (VRH) ε i ε j E V E C E T E T E V E C E T E V E C

Hopping Mechanism & Models

The right panel of Fig. 2.1 represents two variants in the hopping picture. In the classical picture of hopping, called variable range hopping (VRH) the charge carrier is completely localized on donor or acceptor molecules and hops from time to time in the energetic sites that are spatially and energetically separated, while, in Poole-Frenkel picture, the charges are localized in specific states that act as traps and conduction happens through temperature and electric field assisted excitation from the trap states. Due to the mechanism of trapping and de-trapping, this model is called multiple-trap-and-release (MTR), and is applicable to several organic semiconductor systems as the mobility generally observed in organic semiconductors is thermally/electrically activated and exhibits a square root dependence to electric field (µ ∝ √ E, where E is the applied electric field) [4,6] and is proportional to the inverse of temperature (µ ∝ 1/T ).

Poole-Frenkel (MTR) model with its trapping of a charge carrier in energetic sites and temperature/electric field assisted de-trapping throws light on the importance of energetic disorder present in the material and its influence on charge-carrier mobility. Instigation of Gaussian Disorder Models (GDM) by Bässler [7], wherein energetic disorder, that generates energetic trap states in electronically coupled hopping sites helped to understand the impact of energetic disorder on hopping transport. In GDM models, energetic disorder is described by a Gaussian density-of-states (DOS) with a variance that quantifies the amount of disorder present in the system. Monte Carlo simulations were then performed for periodic systems as a function of temperature and electric field to extract the charge carrier mobility. Subsequently, a numerical approach was proposed by Pasveer et al [8] wherein energetic disorder, temperature, electric field and charge carrier concentration can be incorporated into a master equation to compute the charge carrier mobility. This numerical approach also emphasized the impact of charge carrier concentration on charge carrier mobility.

Although, these approaches were able to capture the charge transport mechanism providing a qualitative picture, they were not able to provide a complete fundamental understanding of several intricate parameters/phenomena that govern charge-transport mechanism in organic electronic systems. The need and importance of quantum descriptors that can influence the macroscopic charge-carrier mobility like polaron effects [9,10], local and non-local electron phonon coupling [11,12,13], reorganization energy [14] and polarization of the surrounding medium [4,14,15,16,17] were widely discussed in the organic electronic community. This led to the modifications in modeling mechanisms to include these parameters in charge transport studies through more sophisticated methodologies that include Marcus or Marcus-Levich-Jortner formalisms or Spectral-Overlap-Approach. Several studies were carried on different organic electronic systems that include these effects and a semi-quantitative description of charge transport as a function of the parameters mentioned above were provided [4,14,18,19,20,21,22]. At this juncture of availability of several model descriptors, the most important models that are relevant to the dissertation are presented in the following sections.

Miller-Abrahams model

Bässler [7,23] and co-workers employed Miller-Abrahams formalism in conjunction with Monte Carlo (MC) simulations, within the GDM scheme, to treat disordered systems where the hopping sites are distributed on a lattice of specific periodicity. The transfer/hopping rate in the Miller-Abrahams formalism reads (Eq. 3.1),

k ij = ν 0 exp(-2γR ij ) exp(- ∆G ij k B T ) if ∆G ij > 0 ν 0 exp(-2γR ij ) if ∆G ij < 0 (2.4)
where ν 0 is the frequency factor associated to rate of hopping and γ is the inverse localization radius associated to the electronic coupling matrix elements between adjacent neighbors. ∆G ij is the difference in energy between the initial and final states, given by Eq. 3.2.

∆G ij = j -i + e E • R ij (2.5)
where, e is the elementary charge, E is the electric field vector and R ij is the distance between any two hopping sites i and j, of energies i and j . Owing to the simplicity of the model and ease of application, several investigations were carried out employing Miller-Abrahams formalism [7,23] (Eq: 3.1), in particular for model systems where hopping sites are distributed in a regular lattice defined to have specific lattice constant a. However, Miller-Abrahams model can also be applied for systems where hopping sites are obtained using simulated morphologies (off-lattice models) [1,24].

Marcus Model

Formulated to explain the charge transport mechanism in chemical reactions wherein electron/charge is transferred from one species acting as electron donor to another acting as electron acceptor, Marcus Theory [25,26] has been widely used to compute charge transfer rates in systems where electronic coupling between the neighboring molecules is weak [4,14]. A priori employed in the model is that the charge (electron/hole) is fully localized on donor moiety before the charge transfer, and on acceptor moiety after the charge transfer. The transfer rate in the semi-classical Marcus formalism reads

k ij = 2π h J 2 ij √ 4πλk B T exp - (∆G ij + λ) 2 4λk B T (2.6)
where, ∆G ij is the difference in energy between the initial and final states, J ij is the electronic coupling (transfer integral), λ is the reorganization energy, k B is the Boltzmann constant and T is the temperature. The Marcus formalism describes the charge-transfer process as an Arrhenius activated process, with activation energy given by (∆G ij +λ) 2 /4λ. Reorganization energy comprises (i) an internal contribution (λ i ) associated to the change in geometry of the molecular site i upon charge transfer and (ii) an external contribution (λ e ) which reflects the structural changes in the surrounding medium. In the presence of an applied electric field ( E), ∆G ij includes the contribution from electric field (e E • R ij ). Applicability of Marcus formalism holds true for materials with weak electronic coupling (J < λ/4 ) and when thermal energy exceeds vibrational energy (k B T >> hω ef f , see next section). Marcus formalism is semi-classical in nature as the vibrations are treated classically wherein vibronic coupling is neglected. Vibrations are described by classical harmonic oscillators that are related to coordinates of donor and acceptor units involved in charge transfer and hence are included in the reorganization energy (λ).

Marcus-Levich-Jortner (MLJ)

In many real systems it is often the case that intramolecular vibrational energy exceeds thermal energy at ambient conditions (k B T << hω ef f ), entailing that quantum treatment of vibrational modes become important and Marcus formalism is replaced by MLJ formalism [27,28] wherein energy contribution from quantum mechanical vibration sublevels are incorporated to account for electron-phonon coupling (vibronic coupling) and tunneling mechanism. Inclusion of effective vibration modes allows the correct treatment of quantum level vibrational contribution at room temperature, thereby providing a more accurate description of charge transfer mechanism in comparison to Marcus formalism. Further, all the vibrational modes are merged into one effective vibrational mode with a specific energy and variance, which are described by the effective Huang-Rhys parameters, hω ef f and S ef f respectively. The transfer rate between electronically coupled sites in the MLJ formalism reads

k ij = 2π h J 2 ij √ 4πλ S k B T ∞ n=0 exp(-S ef f ) S 2 ef f n! × exp{- (∆G ij + λ S + nhω ef f ) 2 4λ S k B T } (2.7)
where ∆G ij is given by Eq. 3.2, J ij is the electronic coupling, λ S is the classical contribution to reorganization energy which is mostly due to surrounding medium, hω ef f and S ef f are effective Huang-Rhys parameters.

Reorganization Energy (λ)

Transfer of one charge either a charge donating group to an accepting group, is associated with a reorganization energy which takes into account the interaction of charge with surrounding nuclei and electrons. As electronically coupled molecules are bound by weak van der Waals interactions, reorganization energy can be divided into internal and external contributions. Internal contribution (λ i ) is the energy associated with geometry changes in the molecules during charge transfer process. External contribution (λ e ) is the energy associated by the change in energy of electron donor/acceptor sites during the charge transfer process due to the surrounding medium, caused by lattice distortions or polarization.

Internal contribution to reorganization energy (λ i )

The internal reorganization energy (λ i ) can be computed either by (i) a four-point adiabatic potential (AP) approach [29,30,31,32] or (ii) from normal mode displacements of vibration states of the molecule [14,15].

In the four-point adiabatic potential (AP) approach, the geometry of the isolated monomers is optimized for neutral and charged states. Then the energy of the monomer for both neutral and charged states in lowest energy geometries, E 0 N and E 0

C respectively, as well as the energy of neutral state in the geometry of the charged state, and the energy of charged state in the geometry of neutral state, E N C and E C N , are calculated. A schematic representation is presented in Fig. 2.2. Internal reorganization energy, λ i , is then computed using Eq. 2.8, Internal reorganization energy can also be obtained by expanding the potential energies of neutral and charged states in power series of normal mode coordinates

λ i = λ N + λ C = (E N C -E 0 C ) + (E C N -E 0 N ) (2.8) Charged (C) Neutral (N) E C 0 E C N E N 0 E N C λ C λ N Energy Geometry ΔQ
(Q N & Q C )
, taking into account the contribution of vibrational modes to λ i . In the harmonic approximation, the internal reorganization energy (λ i ) is given by, Eq. 2.9 and 2.10.

λ N i = j λ N j = j hω N j S N j = j k j 2 [∆Q N j ] 2
(2.9)

λ C i = j λ C j = j hω C j S C j = j k j 2 [∆Q C j ] 2 (2.10)
Under the assumption that Q C = Q N , ∆Q j is the displacement along normal mode (NM) j between the equilibrium positions of two electronic states (N and C), k j and ω j correspond to force constant and vibrational frequencies, S j denotes the Huang-Rhys factor of vibrational mode j. (Refer to section 2.3 below for details).

External contribution to reorganization energy (λ e )

The presence of a charge on any given hopping site not only influences the geometry of the host molecular site (which is considered in λ i ), but also the surrounding medium. Two types of effects on the surrounding medium as possible based on the response time, as inertial and inertial-less [33]. Inertial effect is the slow (nuclear) process, associated with the orientation and translation degrees of freedom of the surrounding medium due to presence of charge on a hopping site, whereas inertial-less is the fast (electronic) process, wherein the charge density of the surrounding medium is adjusted due to presence of charge on a hopping site [33].

Actually, the absolute value of external contribution to reorganization energy is debated in the scientific community, as some research groups proposed a value of λ e much smaller than λ i [34], some suggesting λ e λ i [4,14,35], and few others suggesting λ e much larger than λ i [36], entailing that computational methodology employed to compute external contribution and the approximations involved in doing so, would effect the absolute value of λ e . Within the Marcus-Levich-Jortner formalism the contribution of λ e becomes important, as it in included in the denominator (λ S ) [37,38,39]. Earlier studies comparing MLJ and Marcus formalisms demonstrated that mobilities obtained from MLJ formalism are significantly larger than those from Marcus formalism [15]. This difference was attributed to the type of reorganization energy considered in both models. MLJ formalism, in general, includes only reorganization energy from the classical vibronic modes, such that λ S contains only a fraction of the total reorganization energy considered in the original Marcus formalism (λ = λ i + λ e ) [15]. In order to provide a better estimate of calculated mobilities from MLJ formalism, efforts are in progress to include external contributions to reorganization energy such that λ S = λ cl + λ e , by correct treatment of external contributions to reorganization energy [35,34,37,36,40]. External contribution to reorganization energy can be obtained by employing polarizable force fields wherein distributed atomic multipole moments are considered to treat permanent electrostatic interactions [15], but this aspect is will not be discusses in this dissertation.

Electron-Phonon coupling & Huang-Rhys parameters

The presence of charge on a hopping site influences the intra-molecular vibrational normal modes, leading to a coupling between electron and phonon (vibrational) energies. If ω j is the vibrational frequency of any given vibration mode j, then the deformation energy can be written in terms of geometric displacements along the normal mode coordinates (see Eqs. 2.9 and 2.10). If only one vibrational ground state is considered (to be occupied), then the wave-function overlap integral of two harmonic oscillators between zero vibrational state and n th vibrational state of the neutral (N ) form of the molecule, χ N 0 and χ N n , respectively, turns into a Poisson distribution with an expectation value that corresponds to the intensity of the vibronic progression, given by Eq. 2.11, where S is the variance of the time distribution of overlap integral between χ N 0 and χ N n .

| χ N n |χ N 0 | 2 = S n n! exp -S (2.11)
If the normal mode displacement (∆Q) for any given vibrational mode j is zero, then S j → 0 and the transition probability of the electron or hole for this mode from donor state to acceptor state becomes one. However, when ∆Q j = 0, then S j = 0 and the electronic and vibronic excitations are coupled. So, the value of S j , for any given mode j can be considered as the strength of electron-vibration coupling. When all the vibrational modes are considered, electron-phonon coupling energy can be effectively described by the energy of single effective vibrational mode that is related to displacement of normal mode coordinates of neutral and charged states, hω ef f with an effective variance S ef f . S ef f is the Huang-Rhys parameter and hω ef f is the effective vibrational mode for a reorganization energy, such that λ i = hω ef f S ef f , together termed as effective Huang-Rhys parameters. A schematic representation of Huang-Rhys parameter (S ef f ) as a function of normal mode displacement vector (∆Q) and internal reorganization energy (λ i ) is given in Fig. 2.3.

For a molecule with m vibrational modes, the dimensional displacement vector B of the molecule between its neutral (N) and charged (C) state, assuming the harmonic approximation, is given by Eq. 2.12 and Eq. 2.13 and the Huang-Rhys parameter S of the molecule for the vibrational mode m, in its neutral and charged form is given by Eq. 2.14 and 2.15, respectively.

B N m = ω N m h {X N -X C }M 1/2 L m (N ) (2.

12)

Neutral (N)

E C 0 E C N E N 0 E N C λ C λ N Energy Geometry ΔQ (i) Charged (C) Neutral (N) E C 0 E C N E N 0 E N C λ C λ N Energy Geometry S i eff S ii eff S i eff > S ii eff Charged (C) ΔQ (ii)
B C m = ω C m h {X N -X C }M 1/2 L m (C) (2.13)
S N m = 1 2 (B N m ) 2 (2.14)
S C m = 1 2 (B C m ) 2 (2.15)
where X N ,X C correspond to 3N vectors of Cartesian coordinates of equilibrium geometries of neutral and charged species respectively, M is a 3N ×3N diagonal matrix containing atomic masses and L m (N ) and L m (C) are 3N × N matrices consisting the mass-weighted Cartesian displacements of normal coordinates of neutral and charged states, related to the vibration state m.

Using the above quantities, effective frequency for neutral (N) and charged (C) states are then computed as follows,

ω N ef f = Σ m ω N m S N m Σ n S n (2.16
)

ω C ef f = Σ m ω C m S C m Σ n S n (2.17)
Once ω ef f is estimated for neutral and charged states, S ef f can be evaluated using Eqs. 2.18 and 2.19, if the internal reorganization energies are known for neutral and charged states (λ N/C i from Eq. 2.8 or Eqs. 2.9 and 2.10).

S N ef f = λ N i hω N ef f
(2.18)

S C ef f = λ C i hω C ef f (2.19)
S ef f is computed using 2.20 and finally, hω ef f is computed using Eq. 2.21.

S ef f = S N ef f + S C ef f
(2.20)

ω ef f = λ i S ef f (2.21)

Transfer Integrals

As could be noticed from Marcus/MLJ formalisms charge-transfer rate strongly depends on electronic couplings or otherwise called transfer integrals (J). Electronic couplings between two interacting molecules (a "dimer") can be obtained either by (i) the energy splitting in the dimer method or (ii) the projection method.

Energy Splitting in the Dimer (ESD) method

The Energy Splitting in the Dimer (ESD) method is the simplest approach to compute the electronic couplings in organic semiconductors [4,13,14,37]. At a transition point where the excess charge between two interacting molecules (i and j considered as a dimer) is equally distributed (delocalized) over both sites (symmetric dimer) [3], the energy difference, E j -E i , between the unperturbed adiabatic states Ψ j and Ψ i is related to transfer integrals as, Eq. 2.22,

J ij = E j -E i 2 (2.

22)

. This method, rigorously, requires that the geometries of the molecules and that of the dimer at transition state (at the avoided crossing point) are available. However, in practice, the geometry of the neutral dimer is used. Further, ESD method is applicable only in case of symmetric dimers and relies on the Koopmans' theorem, as an additional approximation. The value of electronic coupling for charge transfer (electrons/holes) between two sites i and j is then approximately equal to,

J ij = L+1[H] -L[H+1] 2 (2.23)
.

Where L(L+1) and H(H-1) correspond to energies of LUMO (LUMO+1) and HOMO (HOMO-1) energies obtained from closed-shell configuration of the neutral state of the dimer. A schematic representation of energy level splitting is presented in Fig.

Projection Method

In the projection method [16], one-electron dimer states are defined in terms of localized monomer orbitals. Under the assumption that dimer HOMO and HOMO-1 energies result from the interaction of only monomer HOMOs, then the orbital energies of the dimer can then be defined following a secular equation,

HC -ESC = 0 (2.24)
Where H and S are the system Hamiltonian and overlap matrices in the basis of monomer HOMOs,

H = e 1 J 12 J 12 e 2 S = 1 S 12 S 12 1
The matrix elements e i and J ij that enter in matrix H and S are given by Eq. 2.25 and 2.26.

e i = Ψ i | Ĥ|Ψ i (2.25) J ij = Ψ i | Ĥ|Ψ j (2.26)
The monomer orbitals, Ψ i and Ψ j , considered to obtain e i and J ij are non-orthogonal. An orthonormal basis sets retaining the initial local character of the monomer orbital can be obtained from Ψ i , following Löwdin's symmetric orthogonalization scheme [41]. In a symmetrically orthonormalized basis Eq. 2.24 is given as,

H ef f =   e ef f 1 J ef f 12 J ef f 12 e ef f 2  
where, site energies e ef f and transfer integrals J ef f are given by, Eqs. 2.27 and 2.28.

e ef f 1(2) = 1 2 (e 1 + e 2 ) -2J 12 S 12 ± (e 1 -e 2 ) 1 -S 2 12 1 -S 2 12
(2.27)

J ef f 12 = J 12 -1 2 (e 1 + e 2 )S 12 1 -S 2 12
(2.28)

When S 12 = 0, then J ef f 12 = J 12 obtained from the ESD method (J ij in 2.23). Further, J ef f ij computed from projection method (Eq.2.28) can differ from J ij computed from ESD method (Eq. 2.23) based on the configuration of the dimers. To elicit this effect, earlier studies of anthracene dimers indicate that, when the long axis of molecules are aligned parallel one-to each other (face-on configuration) then both methods give similar magnitude of transfer integrals, whereas when the long axis of one molecule is rotated with respect to the other (face-to-edge configuration) the the transfer integrals obtained from both the methods differ, with ESD method drastically overestimating the value of transfer integral [42], entailing that, the projection method is a better estimator of the absolute magnitude of the transfer integrals.

Band Transport

When interacting molecules are arranged in ordered periodic crystals, the electronic picture is represented by combining the atomic orbitals to provide an equivalent of molecular orbital picture such that the atomic orbital overlap modulates but do not change the initial atomic levels. This is because, in a periodic crystal the interaction between atomic orbitals leads to a reorganization of energy levels that depends on the intermolecular/interatomic (nuclei) distance: the closer the distance between two interacting nuclei, the stronger is the perturbation of the atomic orbitals. The amount of perturbation is in turn a function of energy of the atomic orbital. This results in series of discrete energy bands, each energy band resulting from interacting atomic orbitals, with regions of forbidden energy called "energy gaps", thus providing a band picture, that replaces the atomic/molecular energy levels of the isolated units picture. Energy bands formed by such overlap are periodic in nature, entailing that electron density is identical at equivalent points in the crystal. It follows that the wavefunction at any point x in the crystal and at any equivalent point x + na (n being an integer), are identical, such that Eq, 2.29 follows,

|Ψ k (x)| 2 = |Ψ k (x + na)| 2 (2.29)
Here, a is the lattice constant of a periodic 1D lattice. Bloch's theorem can then be applied to obtain a phase relation of the wavefunction at any periodically connected points by, Eq. 2.30, where k is a wavevector in reciprocal space,

Ψ k (x + a) = 2π x exp(ika)Ψ k (x) (2.30)

Tight Binding Approximation (TBA)

When the wavefunctions that represent the energy bands in a periodic system are approximated from wavefunctions of free atoms, the resulting energy model is known as Tight Binding Approximation (TBA) or Linear Combination of Atomic Orbitals (LCAO) approximation. In TBA, the innermost energy bands, formed from inner electronic shells of atoms, are relatively localized while the energy bands from the outer electronic shells of atoms are diffuse/unconfined. The spread of the band represents the strength of energy perturbation. This entails that the width of a band typically as a function of energy of the band.

The energy of the system can be represented by Schrödinger equation, Eq. 2.31, where the wavefunction Ψ k takes the form, given by Eq. 2.32, where n correspond to the lattice vector ±na and the orbital located at site n is given by, φ α (x ± na), with φ α representing the atomic orbital and c α (k) the expansion coefficient.

HΨ(x) = EΨ(x) (2.31) Ψ k (x) = n α exp(ina)c α (k)φ α (x ± na) (2.32)
Indeed, the wave function from Eq. 2.32 is a Bloch function since it satisfies Eq.2.30, taking into account the Tight Binding Approximation (TBA).

Eigen values of Eq. 2.31, results in energy of the bands as a function of wavevector, with a bandwidth (W) that is equal to the energy difference between the highest and lowest energies of the given band. The variation of band energy in the k space is known as the band dispersion, that takes into account the orbital overlap in the given direction of k vector. The energy as a function of k is given by, Eq. 2.33, where na represent the periodicity of the crystal lattice,

E(k) = E 0 + n=±1,n=±2 E n exp (-ikna) (2.33)
Effective Mass (m * )

As energy in the crystal is quantified with discrete energy bands with forbidden/gap states, the electron or hole will move along the energy bands. Different forces, internal and external in nature, are in action for the electron movement in the bands and the total force acting on the electron can be represented as F = ma, m * is the electron/hole effective mass for a free electron. Kinetic energy of the electron represented as de Broglie particle, is given as E = p 2 /2m, where p = hk. As the electron moment is restricted in the band, the energy of the electron can be approximated to be equal to that of band energy as a function of k (Eq. 2.33)

From the above, it follows that, Energy, E(k) and wavevector k are related to effective mass, through, Eq. 2.34. Further, band dispersion and effective mass are inversely proportional indicating that higher the curvature of the energy band, higher is the band dispersion and lower is the effective mass. Energy of the electron/hole as a function of k vector, at the center of Brillouin zone (Γ point) can be given as Eq.2. 34, where E 0 correspond to the thermal energy of the electron when no external forces are acting and i and j correspond to the dimensionality of the periodic crystal.

E ± (k) = E 0 ± h2 k 2 2m * (2.34)
For a semiconductor, the electronic band at the center of Brillouin zone (Γ point) corresponds to extrema of conduction and valence bands, which are localized at Γ. The energy dispersion at Γ point are then related to the electron (e) and hole (h) mobility in semiconductors. This entails that, the energy of the electron at conduction/valance band edges in the TBA approach, can be approximated to that arising from orbital interaction of nearest neighbors. For a 1D system, the energy of the electronic band at the extrema of conduction or valance band, can then be given by, Eq. 2. 35, where E 0 =< φ 0 | Ĥ|φ 0 > is the onsite atomic energy and

J =< φ 0 | Ĥ|φ -1 >=< φ 0 | Ĥ|φ 1 >
is the electronic or transfer integral with the two nearest neighbors. Further, when band energy close to the bottom of the band is considered, where k → 0, Eq. 2.35, when no thermal energy is considered, can be approximated to 2.36,

E(k) = E 0 -2J cos(ka) (2.35)
on Taylor expansion, gives

E(k) = E 0 -2J -Jk 2 a 2 (2.36)
From Eq. 2.35, it follows that, for a one dimensional case, the width of electronic band W = 4J, and relating electronic band energy and electron energy in the band as a function of wave vector (k), from Eq. 2.36 and 2.34, effective mass (m * ) and band integral (J) can be related through, Eq. 2.37. A schematic of the energy dispersion for 1D tight binding approach is presented in Fig. 2.5. The band energy is presented as a function of k. Within the semi-classical Drude approximation mobility is given by µ = eτ /m * , e is the elementary charge and τ is the collision time of the particle, which corresponds to the life time of the charge-carrier between subsequent scattering events. Therefore, mobility of charge carriers (electrons/holes) can either be computed by considering hopping like charge transport model wherein (mobility) µ ∝ J 2 or by employing band transport model and computing the effective mass wherein µ ∝ 1/m * [43,44,45]. The applicability of hopping or band models (with limiting conditions) are already presented in the previous sections.

h2 k 2 2m * e,h ∼ 2J -Jk 2 a 2 =⇒ m * ∝ 1 J . ( 2

Kinetic Monte-Carlo

Kinetic Monte-Carlo simulations are performed to compute charge carrier mobility of semiconducting organic materials with a home built FORTRAN code. The algorithm employed for KMC simulations is presented Fig. 2.6.

Parameters

As shown in Fig. 2.6, different parameters obtained from quantum chemical calculations and described in this chapter, like reorganization energy (λ), Huang-Rhys parameters (S ef f & hω ef f ), energetic disorder introduced for a material under study are introduced in KMC simulation scheme by selecting the appropriate model to compute the transfer rate: Miller-Abrahams, Marcus, Marcus-Levich-Jortner formalisms. Diagonal energetic disorder and off-diagonal energetic disorder are introduced through disorder on site energies ( ) and disorder on transfer integrals (J), respectively. Different energetic disorder configurations are employed based on the simulation requirement, such as, using calculated site energies/transfer integrals or adding static or dynamic energetic disorder through a Gaussian distribution of specific standard deviation (see Chapters 4 & 5, for further details). KMC simulations are performed as a function of temperature and electric field.

Charge Propagation

For a given KMC run, a starting hopping site is selected randomly and three dimensional periodic boundary conditions are applied to follow the propagation of a single charge hopping along the sample as a function of electric field and temperature. Hopping rates (k ij ) from the hopping site i to all N available neighbor sites j are computed based on the rate equation selected. Instantaneous hopping times are generated by Eq: 2.38,

τ ij = 1 k ij ln 1 ξ (2.38)
where ξ is a random number drawn from an uniform distribution within the interval 0 and 1. Destination hopping site i is selected based on the "first reaction" [46,47] method given from the set of available neighbor sites j as the one having the smallest reaction time τ min ij . Simulation progresses by subsequent hopping steps and the charge is propagated. The simulated elapsed time is increased by τ ii , for each hopping event, and the distance is recorded. Simulation is continued until the distance traveled along the electric field direction (d) reaches a fixed end distance (d E ). Since only one charge is propagating at a time, the simulation targets very low charge densities, as charge carrier density can increase the calculated mobility [8]. For each simulation run (κ), the mobility is obtained as ,

µ κ = d |E| 1 τ κ (2.39)
where τ κ is the hopping times τ ii , ie., the time required by the particle to cover the distance d. The mobility is computed as an average of N KMC simulation runs. Average values are reported either as arithmetic mean, µ = µ A = ( µ κ )/N or as logarithmic mean, µ = µ L = exp ( κ ln [µ κ /N ]), as the latter is a better estimator of mobility in small samples [48] or in samples with high diagonal disorder. The effect of considering µ A or µ L is detailed in the next chapter, as a function of sample size. Compute hopping time :

μ N = (d) N E .( τ ) N μ Calc = ∑ μ κ N =exp( ∑ κ ln μ κ N )

YES

Choose hopping event, with smallest τ : τ ij

Record : d = d+r ij .E τ = τ+τ ij

Update system configuration

Move charge from " i " → " j "

Equilibrium : STOP "KMC N " → i i j=i "N" KMC runs:

t ij =- 1 k ij ln(random[0,1] ij ) d > d E ? KMC N+1 d < d E NO Figure 2.6: Algorithm employed in KMC simulations

Molecular Mechanics

In Molecular Mechanics (MM) methods, the energy of a system of interacting molecules under study is described classically. The energy of the system, given by Eq. 2.40, can be expressed as a sum of energetic terms involving bonded and non bonded interactions that depend only on the atomic coordinates and on a set of empirical parameters. Bonded interactions include bond length, bond stretching, bond/angle bending and torsional terms, while non bonded terms include van der Waals and electrostatic interactions. A schematic representation of the energy terms is presented in Fig. 2.7. Adapted from Ref. [49].

E M M = E Stretch + E Bend + E T orsion + E vdw + E Elec + E Cross (2.40)

Bonded Interactions

Bond stretching energy (E Stretch ) is the energy associated with interaction of two atoms A and B that are chemically bonded. The potential of bond stretching connecting two atoms A and B is given by a harmonic potential (Eq. 2.41) around an equilibrium bond energy E(R AB Eq ), R AB being the distance between the two atoms. However, the harmonic approximation breaks down at large distances, a limitation, that is corrected by Morse potential, given by Eq. 2.42, where D is the bond dissociation energy and κ is force/stretching constant of the bond.

E Stretch = κ(R AB -R AB Eq ) 2
(2.41)

E M orse = D[1 -exp (-κ/2D[R -R Eq ])] 2 (2.42)
Bond bending is the energy associated with bond angle (θ) between three connected atoms A, B and C, where atoms A, B and B, C are bonded. Bond bending energy is most often represented by a harmonic potential (Eq. 2.43), where κ ABC is the force constant associated to bending.

E Bend = κ ABC (θ ABC -θ ABC Eq ) 2 (2.43)
Torsional energy is the energy associated with the rotation around bonded atoms B -C where atoms A-B -C -D are connected sequentially. The torsional angle (τ ) between two planes formed by atoms A, B, C and B, C, D respectively, as presented in Fig. 2.7, can be expressed as Eq. 2.44 where γ correspond to a phase shift, n is an integer number and V n are empirical coefficients.

E ABCD T ors = n V ABCD n 1 + cos(nτ ABCD + γ) (2.44)
Bonded interactions are coupled to one another through internal coordinates as variation in one bonded term induces a variation on another, with the interaction between bond stretching and bond bending being the strongest. An explicit energy cross term (E Cross ) is sometimes introduced. Considering the interaction between bond stretching and bond bending as the most important contribution, E Cross can be represented using Eq. 2.45. Other interactions like stretch-stretch, bend-bend stretch-torsion, bend-torsion etc are generally included in Force Fields that vary on the usage of cross terms.

E Cross(Stretch/Bend) = k ABC (θ ABC -θ ABC Eq )[(R AB -R AB Eq ) -(R BC + R BC Eq ] (2.45)

Non-Bonded Interactions

Van der Waals interaction energy (E vdw ) represents the attractive, repulsive energy between neutral atoms that are not directly involved in bonding. At short distances repulsive forces between the atoms dominate while at intermediate distance the forces are attractive. Considering a quantum picture to represent E vdw , attraction part can be considered as an electron correlation effect, in which fluctuations of an atomic electron density produce a temporary dipole that in-turn induces a complementary dipole on surrounding atoms. Attractive forces are generally referred to as dispersion or London forces. As the distance between atoms become smaller, the electron-electron repulsion energy dominates due to Pauli repulsion and the energy of the system increase exponentially, leading to E vdw being positive, while at intermediate distance, the energy follows a 1/R 6 behavior. As the distance between atoms increases further E vdw 0. The Lennard-Jones potential (LJ) is employed to describe E vdw , given by Eq. 2.46, where σ represent the collision diameter (separation for which energy is zero) and which described the depth of the LJ potential. However, repulsion energy at short distances due to overlap of electronic wavefunctions (electronic structure theory), dies off approximately following an exponential behavior with increase of distance from nucleus. Therefore 1/R 12 from LJ potential can be replaced by the expression containing an exponential form given by Eq. 2.47, where A, B and C are adjustable parameters.

E LJ = 4 [( σ R ) 12 -( σ R ) 6 ]
(2.46)

E vdw = A • exp(-BR) - C R 6 (2.47)
As atomic elements are electro-negative or electro-positive in nature, they tend to attract/give away charge density, giving rise to non-uniform charge distribution in a molecule. A common practice is to represent the distribution of charges as an arrangement of fractional point charges throughout the molecule. These charges when restricted to occupy nuclear centers are referred to as partial atomic charges. Following Coulomb's law, electrostatic interaction between any two molecules is calculated as the sum of interactions between pairs of point charges, given by Eq. 2.48, where N A and N B represent the number of point charges in the two molecules.

E El (r ij ) = N A i=1 N B j=1 q A q B 4π 0 r AB (2.48)
The above mentioned energy terms (bonded and non bonded) are often quantified by quantum descriptors wherein DFT/Semi-empirical calculations are performed on the system of interest and bonded and non bonded parameters from quantum mechanical calculations are introduced in Molecular Mechanics calculations through parameterizing the force field. Molecular Mechanics can be employed to find the equilibrium geometry of the molecule, which is equivalent to optimization/energy minimization at quantum mechanical level (empirical/semi-empirical/DFT).

Molecular Dynamics (MD) Simulations

Molecular Dynamics (MD) simulations are employed to obtain thermodynamic stable configurations of the system of atoms/molecules interacting in a model inter-atomic potential by solving Newton's equations of motions. Starting from an initial system configuration with zero velocities, new atomic positions and velocities are generated within a specified time step. As the time increments (so do the atomic positions and velocities), series of trajectories are generated, each containing the ensemble of atoms/molecules, and the energy of the system evolves till a thermodynamic equilibrium is attained. These MD runs where the energy of the system reaches a thermodynamic equilibrium are typically called "equilibration". Once equilibration is attained, MD simulations are allowed to continue long enough, a step called "production", to analyze the system ensemble for properties of interest, ie., energy, volume, pressure, inter/intra atomic interaction etc, to name a few.

Newtonian equations of motion that are considered in MD simulations can be represented using, Eq. 2.49 .

δ 2 r i (t) δt 2 = F i m i where, F i = - δV (r i , r j , ..r n ) δr i (2.49)
where, force (F i ) and mass (m i ) correspond to atom i and time is denoted by t. V ( r) correspond to the molecular mechanics potential energy (Eq. 2.40) which is differentiable as a function of atomic coordinates r i (see Eq. 2.40). For each time step (∆t) forces are evaluated analytically for all the atoms/molecules and Newton's equations are integrated numerically. Verlet algorithm [50] is widely used to solve these Newton's equations of motion wherein atomic positions (r) and velocities (v) are calculated using Eqs.

Types of MD simulations

Algorithms presented in the previous section describe the time evolution of the system consisting of N particle in volume V , conserving the total energy, generating a microcanonical ensemble (NVE). However, many experimental conditions are performed at specified (constant) temperature and pressure and MD simulations in isobaricisothermal conditions are important to correlate with experiments. Accordingly, both constant temperature (N V T ) and constant pressure (N P T ) MD simulation conditions are presented.

In constant temperature MD simulations, temperature of the ensemble is related to time average of kinetic energy by, Eq. 2.56. By scaling the velocities, temperature is scaled and velocities are multiplied at each time step by a factor (T 0 /T (t)) 2 , where T 0 and T (t) are the desired temperature and temperature calculated from kinetic energy, respectively.

E kin (t) = 1 2 N i m i v 2 i = 3 2 N k B T (t), it follows that, T (t) = 1 3k B N N i=1 m i v 2 i (2.56)
Alternatively, temperature of the ensemble of particles is coupled to an external bath fixed at desired temperature. Velocities are scaled at each step such that the rate of variation of temperature is proportional to the difference of temperature between external bath and ensemble. The rate of temperature is given by Eq. 2.57, where τ is a parameter that determines the strength of coupling between system ensemble and external bath.

dT (t) dt = 1 τ (T 0 -T (t) ) (2.57)
At each step, velocities are rescaled by a factor, λ R , given by,

λ R = [1 + ∆t t ( T 0 T (t) -1)] 1 2
(2.58)

Similar to temperature control, pressure control in MD simulations is achieved either by scaling the ensemble volume such that the pressure is kept constant or by coupling the ensemble to an external pressure bath. The pressure in the system due to N interacting atoms , can be described by virial stress, given by σ V ij : Eq. 2.59, where i & j take values of x, y, z directions of the ensemble of volume

V , R A i & R B i represent the position of atom A & B in the direction i, F AB i
is the force acting between any two atoms A & B in the direction i, m represent the mass of the atom and v i & v j correspond to thermal excitation velocities of the atom along directions i & j respectively.

σ α.β = 1 V ( N i=1 m i v iα v iβ + N i=1 F iα r iβ ) (2.59)
where, α, β = x, y, z, v is the velocity, F is the force and r is the position. For an external pressure bath, the rate of exchange of pressure between the ensemble and pressure bath is given by Eq. 2.60, where τ P is a parameter that determines the strength of coupling between the system ensemble and external bath and P can be calculated employing Eq. 2.61.

dP (t) dt = 1 τ P (P -P (t)) (2.60) P = σ 2 xx + σ 2 yy + σ 2 zz 3 (2.61)

Softwares & Computational Packages

Electronic structure calculations

Gaussian09 [53] and ORCA [54] software packages were used to obtain molecular structures, vibrational frequencies and electronic structures of molecules and molecular clusters, giving access to information like partial charges, orbital energies, transfer integrals/electronic couplings. Level of theory employed for electronic structure calculations along with the specific package used, will be presented appropriately in the dissertation.

Periodic DFT

Periodic DFT calculations were performed with Quantum Espresso [55]. Level of theory employed for plane wave electronic structure calculations will be presented appropriately in the dissertation.

MD simulations

All molecular dynamics simulations in the dissertation were performed using NAMD [56] software package in conjunction with VMD [57] and MD-Analysis [58] suites, used appropriately based on the requirement.

KMC simulations

All Kinetic Monte-Carlo simulations presented in the dissertation were done by a home built FORTRAN-90 KMC simulation code, which has been tested and benchmarked with respect to previous studies using different simulation conditions and formalisms like Miller-Abrahams, Marcus, MLJ rate formalisms. Results of benchmark calculations validating the usability of the code is presented in the next chapter.

CHAPTER 3 CHARGE TRANSPORT SIMULATIONS: BENCHMARKS AND EXAMPLES

In the following sections, charge transport studies employing Kinetic Monte-Carlo (KMC) simulations performed with a custom built FORTRAN code are presented. The objective is to give an overview of the benchmark calculations validating the KMC code for Miller-Abrahams, Marcus and Marcus-Levich-Jortner (MLJ) formalisms, as they are used in the dissertation. The parameters entering in the kinetic models are detailed and discussed when required. For physical aspects and the terminology pertaining to different formalisms, refer to the previous chapter on Theoretical Methodology.

Miller-Abrahams formalism

Bässler [1,2] et al, instigated the Gaussian Disorder Model (GDM) to study the impact of energetic disorder (diagonal and off-diagonal contributions), on charge transport in amorphous and doped organic semiconductors. The morphology of the system consisted of hopping sites that occupy the positions in a cubic lattice with a specified periodicity. Diagonal and off-diagonal contributions to energetic disorder are described by Gaussian density of states (DOS). These KMC studies by Bässler et al were able to capture the main features of charge transport mechanism in organic semiconductors, even describing the Poole-Frenkel nature of hopping transport as a function of energetic disorder. In the dissertation, Kinetic Monte-Carlo (MC) simulations were performed on lattice model by constructing a cubic lattice of 8000 hopping sites with a lattice constant of a = 6 Å. Diagonal disorder, pertaining to site energies, is introduced by drawing a random number from a Gaussian distribution as per Ref. [1,2], with σ = 0 k B T, 2.0 k B T, 2.5 k B T, 3.0 k B T, 3.5 k B T . This is done to show the effect of diagonal contribution on mobility and the Poole-Frenkel nature of charge transport. Diagonal disorder considered is static/fixed in nature (refer to Sec. 3.3.1 for details). Distribution of diagonal disorder that correspond to σ = 2.0 k B T = 51.4 meV for 8000 lattice sites is presented Fig. 3.1.

The transfer rate between hopping sites, i and j, is computed using the Miller-Abrahams formalism (Eq. 3.1). Cut off distance on interaction energy between the hopping sites was set to 20 Å, 2γa was set to 10 and frequency factor ν 0 = 10 12 s -1 following reference [1].

k ij = ν 0 exp(-2γR ij ) exp(- ∆G ij k B T ) if ∆G ij > 0 ν 0 exp(-2γR ij ) if ∆G ij < 0 (3.1) with ∆G ij = j -i + e E • R ij (3.2)
where E is the electric field and R ij is the distance between any two hopping sites i &j and i and j are the energies of donor and acceptor sites, respectively. Kinetic Monte-Carlo simulations are performed as a function of applied electric field with temperature set to 300 K. Results obtained in this dissertation are presented in Fig. 3.2, shows a very good agreement with the reference data [1], validating the KMC simulation scheme for Miller-Abrahams formalism. A decrease in mobility associated with an increase in Poole-Frenkel nature of charge transport can be observed (Fig. 3.2) with the increase in σ . Increase in σ in-fact tends to increase the number of trap states and charge propagation happens through temperature and electric field activation, which leads to increase the Poole-Frenkel factor, β. Positive value of β (positive slope of mobility) indicates the Poole-Frenkel nature of charge transport.

Based on the value of β, two different regions of charge transport can be identified in Fig : 3.2 as a function of σ , diffusion driven (σ > 2k B T ) and drift driven (σ < 2k B T ) charge-carrier mobilities. From Fig. 3.2, it can be also observed that, below a specified limit of electric field, charge transport is governed by electrically and thermally assisted hopping wherein a positive slope of the curves with respect to β is observed. At high electric fields, the charges no longer see the disorder present (trap states) and are drifted along the direction of electric field leading to a decrease in mobility. This holds also true when the diagonal energetic disorder is zero, as no Poole-Frenkel nature of charge transport is observed.

Sample Size and Simulation equilibrium

In order to address the impact of the sample size on charge carrier mobilities, KMC simulations were performed on lattice models with σ (fixed disorder, static) and temperature set to 2.0 k B T and 300 K respectively, as well as the applied electric field. The size of cubic lattice boxes with a lattice constant of 6 Å was increased in terms of number of hopping sites from 8000 to 27000 and 64000. When simulations are stopped at d E = 4 µm no variation in mobility is noticed with respect to reference data [1], irrespective of the size of simulation box considered. However, when simulations are stopped with d E = d box , a decrease in mobility is noticed (Fig. 3.3) when the size of simulation box is increased, wherein the absolute mobility values tends to reach the values obtained when d E = 4 µm. This suggests that, equilibrium is not attained for d E = d box , even with an increased number of KMC simulation runs (10000 KMC runs in comparison to 100 for the earlier case when d E = 4 µm), and the choice of sample size of simulation box influences the computed mobility. It follows that the size of the simulation box should be very large if simulations were to be stopped at the end of the simulation box, so as to attain equilibrium, as the mere increase in number of KMC simulations runs does not influence the mobility for d E = d box . Attainment of equilibrium can be further verified by analyzing the mean energy extracted from the energy profile of the hopping sites during or at the end of the simulation run. If equilibrium is attained, the mean energy should be equal to -σ 2 /k B T , where σ is the initial standard deviation of the distribution of site energies [1]. In order to demonstrate this effect, KMC simulations are performed, with static diagonal disorder, on the lattice model with 8000 hopping sites, by progressively increasing the simulation end distance (d E ), and the energy profile of the hopping sites extracted from the KMC simulations as a function of d E is presented in Fig. When the simulation equilibrium is not attained, the final energy of the system is not saturated and the simulated mobility is overestimated, as could be observed from Fig. 3.3, compared to the equilibrated mobility. This effect is not specific to GDM lattice models but can be observed also in off-lattice models [3]. An alternative approach to check equilibrium, is to quantify the mobility through statistical means as a function of number of KMC simulation runs for the same end distance. This can be achieved by comparing arithmetic mean, µ A = ( µ κ )/N , where, N is the number of KMC simulation runs, with the logarithmic mean, µ L = exp( κ ln µ κ /N ). When equilibrium is attained with µ A = µ L , while far from equilibrium µ L is a better estimator of mobility [4]. The evolution of µ A and µ L as a function of electric field is shown in Fig. 3.5, wherein simulations are performed on lattice model with 8000 hopping sites by varying σ . KMC simulations are performed with 250 simulation runs (N ) and the simulation end distance is deliberately set to 2 µm for each KMC run, to verify the convergence of µ A and µ L . A deviation between µ A and µ L can be noticed with the increase in energetic disorder, indicating that equilibrium is not reached and longer simulations with larger simulation end distance are required. As reported by Bobbert et al, this effect is universal when energetic disorder is present in simulation samples, irrespective of lattice or off-lattice models [4]. All the KMC simulation results reported in the dissertation are run for longer simulation end distances, and it is always verified that µ A = µ L .

0.5 1.0 1.5 2.0

E 1/2 (V 1/2 /cm 1/2 )
10 - 910 -8 10 -7

10 -6 µ (a.

u)

σ 1. 0 σ 1. 5 σ 2. 0

Marcus-Levich-Jortner (MLJ) Formalism

KMC simulations on ideal stacks of parallel pentacene molecules, arranged in one dimensional (1D) π stacked configuration (with respect to Ref. [5], wherein the intermolecular distance is set to 4 Å) were performed using the MLJ formalism [6,7], wherein, the transfer rate reads,

k ij = 2π h J 2 ij √ 4πλ S k B T ∞ n=0 exp(-S ef f ) S 2 ef f n! × exp{- (∆G ij + λ S + nhω ef f ) 2 4λ S k B T } (3.3)
as already described in Chapter2. Hole and electron mobilities were computed with simulation conditions identical to Ref. [5], by setting the transfer integrals for hole and electron transport to 0.132 and 0.104 eV, respectively, effective vibrational frequency hω ef f to 0.2 eV and internal reorganization energies (λ i = λ S ) for holes and electrons to 0.095 eV and 0.133 eV respectively and S ef f = λ i /hω. Fig. 3.6 illustrates a good agreement between electron and hole mobilities obtained in this dissertation with respect to the results from reference [5], validating the implementation of MLJ formalism. It is to be noted that, the mobilities obtained here are 0.6 times smaller when compared to reference [5], due to the variation in KMC algorithm. "First reaction" method is employed in this dissertation where in the reaction rates and the subsequent hopping site are selected based on the fastest possible event where as in reference [5] both forward and backward reaction rates are employed to select the hopping event based on a probability factor applied on the event selection (see reference [5] for further details). However, apart the scaling factor, the electron and hole mobility profiles as a function of electric field, in this dissertation, exactly follow the profile reported in reference [5].

Å Å

.6: KMC simulations results obtained in this dissertation (empty dots) with respect to results from reference [5] (dashed lines) for an ideal 1D stack of parallel pentacene molecules, (dashed lines) as a function of electric field.

Marcus formalism & Effecct of energetic disorder

The Marcus formalism was validated by comparing 2D hole mobility of pentacene crystals with that reported in reference [8]. The transfer rate in the semi-classical Marcus formalism [9,10] reads, (see chapter 2 for details on terminology)

k ij = 2π h J 2 ij √ 4πλk B T exp(- (∆G ij + λ) 2 4λk B T ) (3.4)
Simulations were performed by replicating the pentacene unit cell to 10 × 10 × 10× along each crystallographic axis, using periodic boundary conditions (PBC). The total reorganization energy is set to 0.144 eV (λ = λ i (0.92 meV ) + λ e (0.52 meV ) = 0.144eV ). Transfer integrals for hole transport were taken from reference [11] and were assigned to the the replicated structures between the electronically coupled neighbors accounting for the PBC. In-plane hole mobility of the pentacene crystal with electric field (set to 1.0 × 10 5 V /cm) along the crystallographic ab plane is presented in Fig. 3.7, and is in good agreement with the results from reference [8].

-λ (i+e) --λ (i) • Exp.

Effect of energetic disorder

After testing the validity of the KMC algorithm for Marcus formalism, effect of energetic disorder on charge carrier mobility is studied by considering both diagonal and off-diagonal contributions to energetic disorder. In-plane hole mobility of the pentacene crystal is computed as a function of energetic disorder, employing Marcus formalism with λ and electric field set to 0.144 eV and 1.0 × 10 5 V /cm, respectively, similar to the previous section.

Diagonal disorder

The impact of diagonal disorder on charge carrier mobility and Poole-Frenkel nature of mobility was studied previously for both lattice and off-lattice models [1,2,12,13,14].

In the present dissertation, two different simulation conditions are considered for diagonal disorder: (i) Static disorder: site energies are extracted randomly from a Gaussian distribution with a specific standard deviation (σ ) and are kept fixed during the KMC simulation run, entailing that for a given KMC run the site energies are static in time, and (ii) Dynamic disorder: site energies are extracted from a Gaussian disorder with a standard deviation (σ ) and the site energy distribution is randomized at each KMC hopping event, hence the site energies are dynamic in time. Static and dynamic diagonal disorder are henceforth referred to as S G and G respectively.

Simulations are performed on the pentacene crystallographic structure with both static and dynamic diagonal disorder conditions with a standard deviation on site energies (σ ) set to 2 k B T . The in-plane hole mobilities obtained for both simulation conditions are reported in Fig. 3.8, along with the results when no energetic disorder in considered. Both static and dynamic diagonal disorders reduce the mobility with respect to simulations when no energetic disorder is present. This is due to the fact that diagonal disorder act as trap states and contributes to increase of ∆G in Miller-Abrahams, Marcus, Marcus-Levich-Jortner formalisms, which reduces the hopping rate and consequently the mobility. However, when dynamic diagonal disorder is employed, the absolute mobilities are larger than those obtained with the static case. This entails that, the life time of trap states is reduced by frequent randomization of site energies. When diagonal disorder is static in nature, a fixed percolation network for charge transport is established and the charge carrier follows specific hopping sites, whereas, when the diagonal disorder is dynamic, the percolation network is also dynamic and the charge transport is governed by hopping sites with minimal difference in site energies, i.e, the charge carrier is most likely to visit the hopping sites which occupy the Gaussian tail states. As Poole-Frenkel nature is associated to Multiple-Trapand-Release (MTR) mechanism, it subsequently follows that, Poole-Frenkel nature of the mobility would decrease when diagonal disorder is dynamic, in comparison to the static case, as the number of trap states or life time of trap states is reduced. This effect is illustrated in Chapters 5 and 6, where charge transport in crystalline phthalocyanines and amorphous semiconducting polymers is discussed, respectively.

Off-diagonal disorder

Inter and intra-molecular fluctuations such are variation in positional/conformation degrees of freedom, as a function of temperature, would also impact the transfer integrals. Indeed, Troisi and co-workers [15,16] performed Molecular Dynamic (MD) simulations on series of organic semiconductors and reported that computed transfer integrals fluctuate as a function of the temperature and these fluctuations produce a Gaussian distribution of transfer integrals. This entails that, the effect of thermal fluctuations which induce variations in molecular conformations/degrees of freedom can be included in the simulations as variations of the transfer integrals [12,17,18]. As reported by Cornil et al [17], incorporation of fluctuations on transfer integrals in simulation studies would increase the simulated mobility with respect to simulations where no disorder on transfer integrals is present. The increase in mobility is obtained by replacing the mean value of transfer integral (< J > in simulation study with the mean value of transfer integrals distribution summed up with the variance of transfer integral fluctuations (< J 2 >=< J > 2 +σ 2 ) [17,18].

Simulations are performed on pentacene structure by introducing dynamic off-diagonal disorder through disorder on transfer integrals, wherein each value of transfer integral is varied within 50% of its actual value, such that σ J = 0.5 × J. Two different sets of simulations are performed by replacing the absolute mean value of J with (i) < J 2 > +σ 2 J , similar to reference [17] (labelled as J σ G and (ii) J + x where x is a random number drawn from a Gaussian distribution centered on J with standard deviation σ J , (labelled as J G ). Results obtained for these simulation conditions along with no disorder on transfer integrals are reported in Fig. 3.10.

With respect to simulated mobility when no disorder (diagonal or off-diagonal) is considered (black-circles in Fig. 3.9), an increase in mobility is observed when dynamic off-diagonal disorder is introduced (J σ G , red-squares in Fig. 3.9), which is in line with the work of Cornil et al [17]. However, when the off-diagonal disorder is introduced as the actual dynamic fluctuations of transfer integrals using J G simulation conditions, a decrease in mobility is observed with respect to no disorder conditions (J G , blue-diamonds in 3.9). Further to the above, simulations are also performed on the same system with both diagonal ( S G and G simulation conditions) and off-diagonal (J σ G and J G simulation conditions) disorders with σ and σ J set to 2k B T and 0.5 × J, respectively, to elicit the impact of energetic disorder on simulated charge carrier mobilities. The results are reported in Fig. 3.10. When diagonal disorder is static ( S G ), introducing offdiagonal disorder through both J σ G and J G simulation conditions produces similar results (refer to blue-dashed and red-solid curves, respectively, in Fig. 3.10). However when diagonal disorder is dynamic ( G ), mobility obtained using off-diagonal disorder with J σ G is larger than that obtained using J G simulation condition (refer to bluehexagons and magenta-pentagons curves, respectively, in Fig. 3.10).

In summary, it can be stated that, the mobilities obtained for G simulation conditions are larger than those obtained with S G , whereas mobility obtained with J G simulation condition are lower than those obtained with J σ G simulation conditions. Further, introduction of both diagonal disorder (static or dynamic with σ = 2k B T ) and off-

0°4 5°9 0°1 35°1 80°2 25°2 70°3 15°5 10 15 a b J σ G J Q J G Figure 3
.9: Hole mobility of pentacene crystal as a function of electric field orientation in ab plane when no off-diagonal disorder is considered (black-circles) and when offdiagonal disorder is considered with J σ G simulation conditions (red-squares) and with J G simulation conditions (blue-diamonds). Radial grid correspond to mobility values in units of cm 2 /V s. diagonal disorders (J σ G or J G with σ J = 0.5 × J), simultaneously, reduce the mobility with respect to simulations when no energetic disorder is considered. The differences in simulated mobility for different simulation conditions of diagonal disorder and offdiagonal disorder are further discussed in the subsequent chapters.

CHAPTER 4 AMBIPOLAR CHARGE TRANSPORT IN CRYSTALLINE GROUP IV PHTHALOCYANINES

Introduction

Phthalocyanines (Pcs) comprising nitrogen-linked tetrameric diiminoisoindoline conjugated macrocycle often chelate metalloids/metals through covalent and/or coordination bonds, leading to stable Metalloid/Metal Phthalocyanines (MPcs). MPcs had been widely employed as dyes and pigments [1,2,3]. MPcs with divalent and trivalent metal atoms have also been successfully employed for organic electronic applications that include Organic Thin Film Transistors (OTFTs) [4,5], Organic Light Emitting Diodes (OLEDs) [6,7] and Organic Photovoltaics (OPVs) [8,9,10] due to their chemical stability, strong absorbance, ease of synthesis, opto-electronic properties with good charge-carrier mobility [4,5,6,7,9,10], wherein specific functionality such as pigment coloration, electronic and/or opto-electronic properties can be tuned by modifying the atoms bound to peripheral or bay positions [11].

Of-late tetravalent MPcs consisting of Si,Ge,Sn metalloid/metal cores [12] have emerged as promising candidates for applications in organic electronic devices with ambipolar charge transport characteristics. More specifically, Si and Ge based dichlorophthalocyanines, Cl 2 -SiP c and Cl 2 -GeP c (compound 3, refer to Fig. 4.1), were incorporated into planar heterojunction organic photovoltaics (PHJ-OPV) devices as both electron donating and electron accepting materials when paired with C 60 and pentacene respectively [13]. Cl 2 -SnP c (compound 2) crystals were demonstrated as high performance n-type materials in OTFTs [14]. In comparison to compounds Cl 2 -SiP c and Cl 2 -GeP c, bis(pentafluoro phenoxy) phthalocyanines (compound 4) based PHJ-OPVs showed better performance with improved device characteristics and a photo conversion efficiency (PCE) of 0.45% when paired with pentacene in comparison to 0.35% for Cl 2 -SiP c [13]. Subsequent to the improved performance of compound 4, a series of fluoro phenoxy phthalocyanines (F x OP h -SiP c derivatives, compounds 5 to 8) were synthesized by varying the number and position of fluorine atoms attached to axial phenoxy moieties [15]. These derivatives, which vary in crystalline packing and π -π stacking (direction, distance, orientation) as a function of number and position of fluorine atoms, were incorporated in PHJ-OPVs and experimentally characterized by Lessard and co workers [15]. It was found that compound 7 and compound 8, differing in the position of fluorine atoms on phenoxy group, are performing better in comparison to other F x OP h -SiP c derivatives. Specifically the performance of compound 7 is better when used with α -6T with a PCE of 1.0 % and that of compound 8 when used with pentacene and C 60 with a PCE of 0.32 % and 1.8 %, respectively [15]. Subsequently, further phenoxy phthalocyanines variants were developed by Lessard and co-workers thatinclude bis

(3-methylphenoxy)-Pc [(3M P ) 2 OP h -SiP c, compound 9] and bis(3-iodo-phenoxy)-Pc [(3IP ) 2 OP h -SiP c, compound 10
] with Iodine or CH 3 substitution on the phenoxy unit resulting in variation of π stacked arrangement of the molecules [16]. Further, it was reported by Lessard et al [15] that during the synthesis of

F x OP h -SiP c derivatives, F 2 -SiP c (compounds 1
) is obtained as a by product during the synthesis of other variants. However, F 2 -SiP c is stable and forms crystalline structure with well ordered π -π stacking.

Although tetravalent MPcs are shown as promising ambipolar materials [13,15], all the above mentioned tetravalent MPcs were incorporated into un-optimized PHJ-OPVs and comparison between different compounds was drawn based on the performance of devices and resulting external quantum efficiency. In light of non availability of full fledged experimental analysis to quantify the charge-transport characteristics of all the compounds, a detailed analysis of charge transport characteristics of tetravalent MPc materials is carried out, by performing Kinetic Monte-Carlo (KMC) simulations within the Marcus-Levich-Jortner (MLJ) formalism (KMC/MLJ). The motivation for the study is to provide insights on the ambipolar charge transport characteristics of these materials, as a function of type/number/position of functional groups attached at axial position or the type of metalloid core in tetravalent MPcs, therein attempting to identify the most promising compounds and the factors contributing to better/higher charge carrier mobility and multi-dimensional charge transport.

Illustrated in Fig. 4.1 are the molecules considered for this study, with the molecular reference numbers that will be henceforth used in the subsequent sections. These compounds are grouped into three categories, to establish relationships between molecular structures and charge carrier mobilities within the crystals, as (i) effect of metalloid core where compounds 1 -3 are considered (ii) effect of position and number of fluorine atoms: F x OP h -SiP c derivatives are considered, compounds 4 to 8 and (iii) effect of Iodine/CH 3 substitution on phenoxy groups (compounds 9 and 10). Crystallographic information for the molecules considered in this study is presented for reference in Table . 4.1.

N N N N N N N N X Y Y O F F F F F O F F F F O F F O F F F O F F F O O I 1: X = Si, Y = F 2: X = Sn, Y = Cl 3: X = Ge, Y = Cl 4: X = Si, Y =

Computational Details

Neutral and charged (positive/negative) geometries of the molecules were optimized at the DFT/B3LYP/6-31G(d) level of theory. Subsequently frequency calculations were performed ensuring proper optimization of the geometry. As reported in Chapter 2, Huang-Rhys parameters for neutral (N) and charged (C) states were computed based on the reference [17], and internal reorganization energies were computed using the four point adiabatic potential (AP) approach [18,19,20,21] at the B3LYP/6-31G(d) level. Transfer integrals between first neighbors were computed using ORCA software [22] by employing the projection method, at the same level of theory, by replicating the experimental crystal cell and using a spherical cut off of 20 Å. For molecules where two-fold (almost) degenerate LUMO were observed, the transfer integrals for electron transport (J e ) were computed using

J 2(e) = [J 2 (LU M O,LU M O) + J 2 (LU M O+1,LU M O) + J 2 (LU M O,LU M O+1) + J 2 (LU M O+1,LU M O+1) ]/2 [23].
Experimental crystal cells were replicated 10-15 times along the three crystallographic directions, such that the super cells obtained after replication are cubic with all sides larger than 100 Å and that the total number of hopping sites for each KMC simulation study varies from 2000 to 3350 depending on the number of formula units (Z in Table .  4.1) present in the experimental unit cell. Transfer integrals were assigned to the replicated structures between the electronically coupled neighbors accounting for the periodic boundary conditions (PBC). KMC simulations were performed on replicated geometries of the molecular structures employing PBC using Marcus-Levich-Jortner (MLJ) formalism. Simulations were performed at 300 K.

For all the KMC simulations reported in the subsequent sections, reorganization energy is fixed to 0.2 eV which is used as classical contribution (λ s ) in MLJ formalism, unless otherwise specified. This assumption is due to, (i) the objective of this work is to draw structure-charge transport property relation/trends for both hole and electron transport of different molecular crystals/structures under study and by fixing the reorganization energy the variations from one structure to another can be isolated to that arising from variation in transfer integrals, packing and to some extent on effective Huang-Rhys parameters and (ii) although several methods have been proposed to compute the classical contribution treating polarization from the surrounding medium explicitly or implicitly [17,24,25,26] the absolute values/magnitude of the classical contribution to reorganization energy are debated within the scientific community.

Further, as all the molecules considered here are similar in chemical and physical nature, it can be assumed that the external contribution to reorganization energy is also similar in magnitude for all the compounds. 4.2. Experimental IP of the compounds reported is extracted from the UPS spectra of the thin films, while the EA is extracted from the IP accounting for the optical band gap [13,15]. For all the molecules, the absolute magnitude of calculated IP are systematically larger than those reported experimentally, while that of EA are lower. This behavior is expected since the theoretical values are calculated for isolated molecules, while experimental values are obtained from thin films. In organic semiconductors, the local charges are stabilized by the polarizable environment, which induces a lowering of IP and an increase of EA compared to gas phase values. Further, absolute magnitude of HOMO energies are consistently smaller than that of IP whereas that of LUMO are consistently larger than EA. These differences can be attributed to the choice of DFT functional employed. Further, when molecular orbital picture is considered, all the compounds show similar distribution of molecular orbital coefficients for frontier orbitals (HOMO and LUMO), which are spread along the phthalocyanine backbone (refer Appendix 1 for details) Internal reorganization energies for hole transport are extremely small for all the molecules, when compared to other well studied crystalline organic semiconductors like rubrene (146 meV [27]) or pentacene (95 meV [27,28,29]), however similar to the values reported for Zn-Phthalocyanines and its derivatives [30]. Further, internal reorganization energies for hole transport (30-70 meV) are systematically smaller than those for electron transport (200-300 meV), despite small variations between the molecules (compounds 1 to 10). Consequently, the effective Huang-Rhys parameters for electron transport are larger than that for hole transport, hinting to electron mobilities being smaller than hole mobilities. Further, hω

Results & Discussions

(e/h) E

(eV ) and S

(e/h) E are relatively large for F x OP h -SiP c derivatives (compounds 4 to 8) when compared to Si, Sn and Ge variants (compounds 1 -3) or when compared to CH 3 /I variants (compounds 9 and 10), for both hole and electron transport. The relative difference between effective Huang-Rhys parameters when different categories of compounds are considered (ie., effect of core, number/position of fluorine atoms or CH 3 /I variants) or when small variations within the same category of compounds, would lead to relative differences in charge-transport characteristics. The impact of relative differences of effective Huang-Rhys parameters on charge carrier mobilities are commented in the subsequent sections with reference to Table 4.3 and Table 4.4 where and when required. Calculated transfer integrals (J Q , calculated using ORCA software [22] at B3LYP/631(d) level of theory) for hole and electron transport, between the reference molecule and its first neighbors (reported in crystal units) for all the compounds are reported in Table . 4.5. (1,0,0) 19 25 (0.5,0,1) 04 08 (0.25,1,0.5) 08 10 10 (0,0,1) 122 87 (1,0,0) 05 04 (0,0.5,0.5) 01 02

Effect of energetic disorder

Before venturing into the study of ambipolar charge transport characteristics of different compounds considered, emphasis is put on understanding the effect of energetic disorder on the charge transport behavior and magnitude of absolute mobility. KMC simulations are performed on compound 1 (Fig 4 .2), along the π stack direction as a function of electric field. Energetic disorder, acting through diagonal and off-diagonal contributions, play a vital role in estimating the charge transport characteristics obtained from KMC simulations. Emphasis was given to the role of static diagonal disorder within the Gaussian Disorder Models (GDMs) during the early simulation studies by Bässler and co workers [31,32]. In the present dissertation, dynamic diagonal disorder, wherein site energies are allowed to vary in time (within the initially considered Gaussian disorder) for every hopping event, is also included for the study. Moreover, Troisi and co-workers [33] have shown that transfer integrals also follow a Gaussian distribution and vary as a function of temperature. This entails that thermal fluctuations that induce variations in molecular positional and conformational space are bound to vary the transfer integrals as reported by Bredas et al [34]. Cornil et al [35], incorporated these fluctuations on transfer integrals by replacing the mean square of the transfer integral (< J > 2 ) by the mean value of the square of the transfer integral (< J 2 >) summed up with the variance of the transfer integral fluctuation (σ 2 J ) [35,34], resulting in the increases of the absolute value of simulated/calculated mobility. In this chapter, a different methodology is employed with respect to the work of Cornil et al [35] wherein each transfer integral associated to a specific hopping event (J Q ) is varied with 50% fluctuation of its current value. This is done by drawing a random number from a Gaussian distribution with a mean value of J Q and standard deviation (σ J ) of 0.5 × J Q and summing up with the value of J Q , which is then used in the KMC simulations, as J G , where J Q is the calculated transfer integral between sites i and j (J Q = J ij ). However, for the sake of completeness (in the present section) dynamic off-diagonal disorder is also introduced by actually replacing < J 2 Q > by

< J 2 Q + σ 2 J >
, similar to the work of Cornil et al [35]. Troisi et al [36] reported that fluctuations of transfer integrals (σ J ) in single crystals can be of order of the absolute value of transfer integrals, whereas Cornil el al [35] reported that the variations can also be lower. In both studies [36,35] transfer integrals were computed from dimers extracted from Molecular Dynamics simulations. In this dissertation, experimental crystal structures are taken from the group of Lessard et al (following the references [13,15,16]) and no Molecular Dynamics simulations were performed to account for the actual fluctuations of transfer integrals as a function of simulation temperature and the corresponding variance. Therefore, a fluctuation of each transfer integral to 50 % of its actual value is considered to be a sensible approximation, to represent the fluctuations on transfer integrals. These assumptions therein can act as scaling factors towards estimating simulated charge-carrier mobility close to experimental observations. However, this would only provide a semi-quantitative estimate as experimental observations/results and efficiencies are bound to several factors like purity of the material, defects/impurities, charge recombination and charge separation etc.

Diagonal disorder (static/dynamic) was incorporated through disorder on site energy ( ) and dynamic off-diagonal disorder was incorporated through disorder on transfer integrals, using the parameters given in Table-4 Gaussian distribution of site energies with σ = 51.4 meV , dynamic in time: disorder generated at every hopping event during the KMC simulation run J Q Calculated transfer integrals. No off-diagonal disorder is considered. 0.5 1.0 1.5 2.0 2.5 The effect of diagonal energetic disorder on simulated mobility, as presented in Fig. 4.3, provides two extreme cases with maximum and minimum mobility for 0 |J σ G (red curves) and S G |J G (magenta curves) simulation conditions respectively, at the same electric field, temperature, reorganization energy and Huang-Rhys parameters. It is of specific interest to note the variation of mobility when diagonal disorder is dynamic ( D G , cyan curve) with respect to static diagonal disorder ( S G , green curve), when no off-diagonal disorder is introduced (J Q ). An increase in mobility and decrease in Poole-Frenkel nature are observed when D G is used. Typical Poole-Frenkel nature of charge transport, wherein µ ∝ exp(β × √ E) with β being the Poole-Frenkel factor for the electric field (E), is observed for

J σ G J Q is replaced by J 2 Q + σ 2 J , introducing off- diagonal disorder during the KMC simulation run. σ J = 0.5 × J Q , is the standard deviation of Gaus- sian distribution centered at J Q . J G J Q is replaced by J Q + x,
E 1/2 ( × 10 3 V 1/2 /cm
0 |J σ G 0 |J Q 0 |J G D G |J σ G D G |J Q D G |J G S G |J σ G S G |J Q S G |J G 0.5 1.0 1.5 2.0 2.5 E 1/2 ( × 10 3 V 1/2 /cm 1/2 ) 0 1 2 3 4 µ (cm 2 /Vs) 0 |J σ G 0 |J Q 0 |J G D G |J σ G D G |J Q D G |J G S G |J σ G S G |J Q S G |J G (a) λ = 0.1 eV (b) λ = 0.2 eV
D G |J σ G , S G |J σ G , S G |J σ G
and S G |J G simulation conditions with the corresponding β values, in (cm/V ) 1/2 units, of 0.25, 1.42, 1.40 and 0.91, respectively (obtained using KMC/MLJ formalism for λ = 0.2eV , Fig.

4.3[b]: Right). The decrease in Poole-Frenkel factor for D

G simulation conditions with respect to S G , entails that when diagonal disorder is dynamic, the apparent diagonal disorder as seen by the charge carriers (electrons/holes) diffusing through the hopping sites is lower with respect to that when diagonal disorder is static. Subsequently, it follows that, the increase in mobility and decrease in Poole-Frenkel factor for charge carrier mobility when diagonal disorder is dynamic in nature can be attributed to the charge carrier (electrons/hole) visiting less often those hopping sites that act as trap states or the life time of trap states being reduced, eventually lead to an increase in mobility. More details on the effect of S G and D G simulation conditions are discussed in the next chapter (Chapter 5).

Introducing dynamic fluctuations of transfer integrals through a Gaussian distribution (J G simulation conditions), decreases the mobility with respect to the mobility when no off-diagonal disorder is considered, whereas introducing the same fluctuations (σ J ) through the variance (J σ G simulation conditions), increases the mobility, which is in line with the work of Cornil et al [35]: compare red (J σ G ) and crimson (J G ) curves with black curves (J Q ) in Fig. 4.3. This trend also holds true when diagonal disorder ( S G or G ) is introduced on top of the off-diagonal disorder.

In order to elicit the decrease in mobility when J G simulation conditions are used with respect to J σ G simulation condition, KMC simulations with MLJ formalism are performed by setting all the transfer integrals that are not along the π stack direction to zero and progressively increasing the percentage of disorder on transfer integrals along the π stack direction (J π , along a axis for compound 1). This is to isolate the contribution of nearest neighbors with non-zero transfer integrals to the mobility and to understand the impact of drawing random numbers for the Gaussian distribution and adding it to the value of J π .

Mobility computed with J G simulation conditions, with diagonal disorder set to zero ( = 0), by progressively increasing the percentage of disorder on transfer integral (i.e. increasing the value of x in σ J ) along the π stack direction, is presented as a function of η =< J > /σ J in Fig. 4.4-[a], with electric field and temperature set to 1.0 × 10 6 V /cm and 300 K respectively. Also, shown in Fig. 4.4-[b], is the mean of squares of transfer integrals summed up with the variance ( J 2 π + σ 2 ) as a function of percentage of disorder on J π , wherein J 2 π + σ 2 is extracted from KMC simulations as a function of percentage of fluctuations on the transfer integral, represented as J σ π . From Fig. 4.4-[b], it can be observed that J σ π , as a function of percentage disorder on J π , extracted from the simulations where dynamic off-diagonal disorder is introduced through J G simulation conditions, reproduces the value of √ J 2 + σ 2 exactly to that value if J σ G simulation condition is to be used. However, a decrease in mobility is observed until η = 1.0 (when σ G is increased), whereas a subsequent increase in σ J increases the mobility (see Fig. 4.4-[a]). Decrease in simulated mobility for J G simulation conditions with respect to J Q and J σ G simulation conditions (refer to Fig. 4.3) is expected as the hopping time hopping time, proportional to 1/ < J > 2 , to increase with respect to the two previous cases (J Q and

J σ G ), since ∞ -∞ 1 J 2 exp[(J-J ij ) 2 /(2σ 2 J )] √ 2πσ J dJ is always higher than 1/J 2 ij (J Q case) and 1/(J 2 ij + σ 2 J ) (J σ G case)
, and mobility to decrease accordingly. Further to the above, the variation in mobility as a function of η also depends (weakly) on reorganization energy (λ S ) and electric field considered. For fixed electric field (1.0 × 10 6 V /cm) and temperature set to 300 K, a decrease in mobility, as a function of η, is noticed when λ S is increased, which is also the case when electric field is increased at fixed reorganization energy (λ S = 0.2 eV ). The mobility plots as a function of η presented in Fig. 4.5 are normalized with respect to the mobility value when no dynamic disorder is considered, ie. with 0 |J Q simulation conditions (refer to Table However, all the variations considered above, ie., with variation in parameters -diagonal and off-diagonal disorder, give similar results, in the region of electric field considered, with variation in absolute value of mobility. Further, when diagonal disorder is introduced along with dynamic off-diagonal disorder on transfer integrals a decrease in mobility with respect to 0 |J σ G , 0 |J Q or 0 |J G simulation condition was observed. This is not surprising as presence of diagonal disorder is bound to reduce the absolute value of mobility [31,32].

Hole and electron mobilities of Compound 1 computed as a function of electric field as reported in Fig. 4.6, using MLJ formalism for different simulation conditions with temperature and reorganization energy (λ S ) set to 300 K and 0.2 eV respectively, shows a systematic variation of absolute values of mobilities, in line with that reported in Fig. 4.3. Further, in-plane mobility is also computed for compound 1 employing different simulation conditions, using MLJ formalism and setting the electric field and temperature to 1 × 10 6 V /cm and 300 K, respectively, as presented in Fig. 4.7, wherein diagonal and off-diagonal disorders are varied as per Table . 4.6. 0.5 1.0 1.5 2.0 2.5 between different factors (∆G ij , S e , hω e , J ij , λ S ) that enter into MLJ formalism. It is to be further noted that, (i) the magnitude of hole and electron mobilities obtained using S G |J σ G simulation conditions are always larger than those obtained using S G |J G and (ii) anisotropy in charge transport is always preserved for different simulation conditions, as reported also in Chapter 3.

E 1/2 (10 3 V 1/2 /cm 1/2 ) 10 -1 10 0 µ (cm 2 /Vs) 0 |J σ G 0 |J Q 0 |J G D |J σ G D |J Q D |J G S |J σ G S |J Q S |J G 0.5 1.0 1.5 2.0 2.5 E 1/2 (10 3 V 1/2 /cm 1/2 ) 10 -1 10 0 µ (cm 2 /Vs) 0 |J σ G 0 |J Q 0 |J G D |J σ G D |J Q D |J G S |J σ G S |J Q S |J G (a)

Effect of Huang-Rhys parameters & MLJ formalism:

Hole mobility obtained from MLJ formalism for compound 1 is consistently smaller than that obtained with Marcus formalism (see Fig. 4.3), for any disorder configuration considered. To further elicit the effect of Huang-Rhys parameters on charge transport, hole and electron mobilities along the π stacked direction (a axis), are computed for compound 4. For compound 4, transfer integral for electron transport along the π stack direction is larger compared to that for hole transport with values of 43 meV and 29 meV respectively (see Fig. 4.9 for details). So, one would expect that charge-carrier mobility for electron transport is larger than that of hole transport. KMC simulations are performed with λ S = λ M = 0.2 eV , where λ M is the reorganization energy used for Marcus formalism, temperature set to 300 K and with D G |J G simulation conditions. Electron and hole mobility for compound 4, along the π stack direction, as a function of electric field using both Marcus and MLJ formalisms are reported in Fig. 4.8.

Electron mobility obtained using Marcus formalism is larger than hole mobility. However, when MLJ formalism is used, electron mobility is smaller than hole mobility, as the value of S e E (2.319) is larger in comparison to S h E (0.644), leading to the reversal of absolute value of mobility of electrons and holes with respect to Marcus formalism, further, eliciting the impact of coupling the quantum treatment of effective vibrational mode with the charge transport characteristics and the importance of considering the quantum treatment of effective vibrational mode towards better estimation of charge carrier mobilities. 0.5 1.0 1.5 2.0 2.5 3.0

E 1/2 ( × 10 3 V 1/2 /cm

Charge transport properties of phthalocyanines

In-plane hole and electron mobilities are computed using MLJ rate equation and KMC simulation scheme, setting λ S to 0.2 eV, with Huang-Rhys parameters from Tables 4.3 and 4.4. All simulations are performed at temperature of 300 K and electric field E = 1 × 10 6 V /cm with D G |J G simulation conditions. (refer to Table 4.6 for details). This choice of electric field and simulation conditions is due to the fact that in the region where electric field is between 0.5×10 6 V /cm and 1.5×10 6 V /cm, the variation in mobility due to static and dynamic diagonal disorder is smaller for compound 1 (refer to Figs. 4.3 and 4.7) and this is expected to hold true for the other compounds, as all the compounds considered are similar in nature. Therefore, simulated mobility would depend less on the choice of nature of diagonal disorder considered. The maximum mobility values for all the compounds along the principle crystallographic axis (a, b and c) directions, obtained using D G |J G simulation conditions, are reported in Table-4.7. Also reported in Table-4.7 is the classification of the compounds with respect to their 1D, 2D and 3D charge transport characteristics for both electron and hole mobilities. The classification is done as, (i) 1D mobility: µ i > 2µ j and µ i > 2µ k , (ii) 2D mobility: µ i ≤ 2µ j and µ i > 2µ k and (iii) 3D mobility: µ i ≤ 2µ j and µ i ≤ 2µ k , where i, j and k correspond to the crystallographic directions of charge transport with µ i > µ j > µ k . Table 4.7: Maximum electron and hole mobility along the crystallographic axis directions for all the molecules considered for the study, to present a overview and quick reference. Mobility is reported in units of cm 2 /V s. Highest mobility values (electron and hole) for each compounds are represented in bold. Also represented in bold is the relative dimensionality of the compounds, wherein the corresponding dimensionality of the compound is higher in comparison to the others. In the subsequent sections, all the compounds (1 to 10, see Fig. 4.1) are categorized based on the type of metal/metalloid core, type/number/position of functional group attached to at the axial position. For all compounds, non-zero transfer integrals (in meV) and the corresponding in-plane mobility values along the represented crystallographic planes are reported. For the in-plane mobility plots reported in the subsequent sections, Red (dashed-line) represent hole mobility and Blue (solid-line) represent electron mobility, with grid representing the absolute mobility values in units of cm 2 /V s. Details of crystalline packing with relative distance between the molecules for ab, ac and bc crystallographic planes and the absolute values of non-zero transfer integrals along with PBC vector and scalar distances for all the compounds are provided in Appendix 1.

Molecule

Category-1: Effect of metalloid core

To illustrate the impact of the metal/metalloid core, electron and hole in-plane mobilities of compounds 1, 2 and 3 are reported in Figs. Transfer integrals for compound 2 along the π stacked direction for holes and electrons are 36 meV and 34 meV respectively whereas for compound 3 they are 33 and 30 respectively, which is due to the smaller π stacking distance of compound 2. Also, for compound 2 the in-plane mobility along ab, ac and bc planes are marginally higher than that of compound 3. For hole transport, as the effective Huang-Rhys parameter (S h E ) is only slightly larger than that of compound 3, the increased mobility of compound 2 can be attributed to the higher transfer integrals of compound 2. Whereas for electron transport, the increased mobility of compound 2 with respect to compound 3 can be attributed to both lower effective Huang-Rhys parameter (S e E ) and higher transfer integrals along a and b axis, and to lower S e E along c axis, as the transfer integral along c axis for both compounds is the same. In comparison to compounds 4 and 5, compound 6 does not show 2D charge transport characteristics, as the variation in relative orientation of molecules leads to smaller transfer integrals for compound 6, which in turn reduce the mobility (refer to Appendix 1 for details on inter-molecular distances and relative orientation of molecules).

Compounds 7 and 8 are compared in-terms of position of fluorine atoms and differences can be noticed between them in terms of π stacking and direction of π stacking, as reported above. Further compound 8 shows higher transfer integrals along the π stack direction. Again, the difference between effective Huang-Rhys parameters is small between these molecules to highlight the contribution originating from them. However, variation in position of fluorine atoms is found to strongly impact the charge transport characteristics. For uniaxial charge transport along the π stack direction, compound 7 (π stacking along a axis) and compound 8 (π stacking along c axis) exhibits similar electron mobilities, while the corresponding hole mobilities of compound 7 are lower. Also, electron mobilities of compound 7 along other two crystallographic axes are higher than that of compound 8. Electron mobilities for compounds 4, 5, 6 and compounds 7 and 8 are smaller than hole mobilities as the value of S e E are larger than S h E . Further, S e E and S h E for compounds 4, 5, 6 decrease progressively (refer to Tables 4.3 and 4.4), which is also the case for the transfer integrals corresponding to electron transport along the π stack direction. However, the corresponding transfer integrals for hole transport of compound 5 are larger than compounds 4 and 6, leading to higher hole mobility of compound 5. Therefore, a progressive decrease of relative electron mobilities for compounds 4, 5, 6 can thus be observed, with respect to the corresponding hole mobilities. This entails that, both the relative difference between transfer integrals (corresponding to electron or hole transport) and the relative difference between S e E and S h E of the compounds 4, 5 and 6 influence the magnitude of the electron and/or hole mobilities. This effect is more significant for compound 6, wherein the relative difference between transfer integrals corresponding to electron and hole transport is small while the relative difference between S e E and S h E is significantly large, eventually leading to a significant decrease in electron mobilities when compared to hole mobilities. Further to the above, comparing compounds from category-2 with those from category-1 (compounds 1, 2, 3), 1D hole transport in category-1 compounds is larger with respect to absolute mobility values along the π stack direction, where as 1D electron transport of compounds from category-01 is almost similar to that of compounds 7 and 8.

Category-3: Effect of CH 3 /I of the phenoxy groups

The effect of methyl or iodine functional group in the para of phenoxy unit on the charge transport characteristics of phthalocyanine derivatives is discussed in the following, comparing compounds 9 and 10. Both compounds show similar crystallographic symmetry but with variation in packing and π stacking. Compound 9 exhibits a slip stack packing along the a axis whereas compounds 10 exhibits an tilted columnar π stacking along the c axis. As a result, the transfer integrals are higher for both electron and hole transport in compound 10, while very small change in effective Huang-Rhys parameters is noticed, with respect to compound 9. In-plane mobilities of compounds 9 and 10 are shown in Fig. 4.17 and 4.18 respectively. Compound 9 shows (almost) shows lower 1D electron and hole mobilities in comparison to compound 10. Further, compound 10 shows the highest mobility of all the molecules that are considered in this study. Effect of highest transfer integrals resulting from efficient packing and reasonably smaller effective Huang-Rhys parameters for compound 10, when compared to the remaining compounds, can be considered as contributing factors for high mobility. Further, relative difference between electron and holes mobilities for compound 9 is the smallest with respect to other compounds (compounds 1 to 10), as compound 9 shows almost similar transfer integrals for electron and hole transport, in terms of magnitude of J and the relative difference between S e E and S h E for compound 9 is smaller than the remaining compounds.

Summary

KMC/MLJ based semi-quantitative theoretical assessment of tetravalent MPcs to establish the structure-charge carrier property trends for variations in metalloid core and type,number/position of atoms on the axial phenoly group is reported, providing a detailed assessment on ambipolar charge transport characteristics.

It generally emerges that, it is difficult to estimate the changes in mobility and to guess the 2D/3D mobility just by considering the crystalline structure. KMC simulations with different energetic disorder contributions brings semi-quantitative prediction of carrier mobility along with providing charge transport anisotropy measurements, which is often difficult to measure experimentally. In line with the same, KMC simulations indicated that, Tetravalent MPcs with variation in metalloid core and variation in position, number, functional units on the axial position exhibit very subtle variations in molecular packing and relative orientation that in-turn lead to large variations in charge transport characteristics. Further, with small variation on the axial functional groups of tetravalent MPcs large variations in charge transport can be obtained, with respect to 1D/2D/3D ambipolar charge transport. These results with a clear relation between structure and charge-carrier mobility, therefore are of importance in paving way towards optimizing the design criteria of the compounds towards developing efficient three dimensional ambipolar organic semiconducting materials. However, experimental characterization of these materials from the group of Lessard and co workers as reported in references [13,15,16] cannot be directly compared with the results from KMC simulations. Differences with respect to experimental observations would arise from different experimental factors like charge generation, separation, recombination, dissociation and also from defects, work-function etc which can affect the efficiency of these devices.

CHAPTER 5

ENERGETIC FLUCTUATIONS IN AMORPHOUS SEMICONDUCTING POLYMERS : IMPACT ON CHARGE CARRIER MOBILITY

Introduction

Amorphous semiconducting polymers, due to their proven versatility, offer substantial advantages in terms of ease of processing and low production costs, with promising applications in the field of organic electronics, such as solar cells, light emitting diodes and field effect transistors. Recent investigations suggest that charge transport in conjugated polymers is not simply governed by the degree of crystallinity, but in general by the presence of aggregates with short-range order, that promote inter-chain charge transport. [1] Although counter-intuitive, experimental evidence that aggregation of polymer chains with sufficient short-range order and interconnected domains, in a globally amorphous structure, is sufficient to promote charge transport [2,3] points towards the refinement of models conventionally employed for describing charge transport in crystalline organic semiconducting systems.

As polymer chains are bound by weak interactions with multiple degrees of freedom, parameters influencing charge carrier mobility that should be included in the modeling of charge (electron/hole) transport include (i) the length of polymer chains, (ii) the conformational freedom along the chain subunits, (iii) the conformation space of the polymer, (iv) the number of possible neighboring units, and (v) the orientation of each chain with respect to the direction of external electric field. [4,5,6,7,8] The theoretical frameworks developed to address the role of morphology on electronic transport in semiconducting polymers demonstrate that intra-chain transport domi-nates at shorter length scales with high mobilities, 10-100 cm 2 /V s, while inter-chain transport prevails at larger displacements, with mobilities a few orders of magnitude lower. The slowest inter-chain transport defines the limiting charge carrier mobility attainable at mesoscopic device scale. [9,5] Such theoretical approaches, in conjunction with Molecular Dynamics (MD) simulations, where reasonably accurate polymer morphologies can be obtained, [10,11] can provide substantial information about the effect of polymer morphology and local structure on charge transport properties of the polymeric system under study. [12].

The main complexity originates from the need of incorporating the whole morphological variability of the macromolecular system into the charge transport modelling. Indeed a theoretical study at the atomistic level, which considers both intra-and inter-chain transport and explicitly takes into account all possible fluctuations over many timescales and morphologies, is a behemoth task, both in terms of human effort as well as computational resources. [13] Such frameworks are particularly unpractical for preliminary material screening, whereas a simplified methodology, able to provide reasonably accurate results with respect to experimental observations at low computational cost, would be more appropriate.

In what follows, a step forward towards realization of this objective is described, in which, with a combination of MD and Kinetic Monte Carlo (KMC) simulations, computational predictions are systematically compared with experimental temperaturedependent mobilities. As a case study the alternating [(9,9-dioctylfluorenyl-2,7-diyl) (4,4 -(N-(4-sec-butylphenyl))] diphenylamine copolymer (TFB) (Fig. 5.1), is considered. TFB is an amorphous polymer [14] with relatively high hole mobility [15,16] that has been the subject of several fundamental studies, [17,18,19,20,21] as well as employed in practical applications such as light emitting devices [22] and solar cells. [23,24,25].

The computational approach is based on two main and interdependent assumptions: (i) that the holes are localized along short segments of the polymer chain, (ii) and that hole transport can be described as a hopping process between those oligomeric units. Evidence supporting these hypotheses is provided, and their consequences on calculated charge carrier mobilities are discussed in detail.

Computational details

Molecular Dynamics simulations

For simulations, focus is given to TFB oligomers composed by four 9,9-dioctylfluorene (FLU) units, alternated with five butyl triphenylamine (TPA) co-monomer units (see Fig. 5.1). This specific oligomer length was chosen as an upper limit to the con-jugation length of TFB, as estimated on the basis of quantum chemical calculations reported by Sancho-García et al. [26] This choice of the rather popular oligomer approach [10,27,28,29] in modeling charge transport in TFB, was also motivated by practical reasons, such as the possibility of dealing with relatively small box sizes and short equilibration times. A united-atom [30] force field is employed, parametrized with quantum chemical calculations for atomic charges and torsional potentials, as detailed in Appendix 2. All molecular dynamics simulations were performed with the NAMD software [31]. In order to simulate amorphous TFB samples, the starting geometries for the simulation were built introducing a certain degree of randomness, according to the following procedure: (i) the starting geometry was generated by placing the oligomers with random orientations on the nodes of a cubic mesh of 5×5×5 sites, with a lattice constant of 80 Å; (ii) the low-density sample was compressed at 1000 K and 100 atm up to a rough volume stabilization, that occurs in about 1 ns;

(iii) high-temperature annealing was performed for 10 ns of simulation at 1000 K and 1 atm; (iv) equilibration was carried out for 20 ns of simulation at T M D =300 K and 1 atm; (v) this procedure was applied to simulate four different and uncorrelated samples composed of 125 oligomers (33380 united atom centers). A snapshot of one of the four samples after equilibration at 300 K is shown in the right panel of Fig. 5.1.

The final density of the four samples is 0.987 ± 0.006 g/cm 3 . It is worth noting that since the equilibration temperature is well below the experimental glass transition, i.e. 393-415 K depending on film thickness [14], significant changes on the morphology are not expected at lower temperatures. For this reason, and for keeping the approach computationally inexpensive, all charge transport simulations were conducted on the 300 K morphologies, neglecting the density variation with temperature.

Charge transport model

Charge transport in amorphous TFB is described with a hopping model. Specifically, Marcus formalism is employed for hopping rates, [32] (see Chapter 2, for details). The charge transfer rate in the semi-classical Marcus formalism reads:

k ij = 2π h J 2 ij √ 4πλk B T exp - (∆G ij + λ) 2 4λk B T (5.1) with ∆G ij = j -i + e E • R ij (5.2)
where e is the elementary charge, i and j are the energies sites i and j involved in the charge transfer, ∆G ij is the energy difference between the initial and final states, J ij is the electronic coupling (also known as transfer integral) λ = λ i + λ e is the reorganization energy, k B is the Boltzmann constant and T is the temperature.

In the presence of an applied electric field ( E), ∆G ij also includes the contribution e E • R ij , R ij being the distance between hopping sites. As positions of the hopping sites, the center of the Mulliken charges evaluated from single point AM1 calculations for positively charged oligomers at their instantaneous MD geometries is used. The reorganization energy is constituted of (i) an internal contribution (λ i ) associated to the change in geometry upon charge transfer and (ii) an external contribution (λ e ) which reflects the structural changes in the surrounding medium. The intramolecular part was set λ i = 0.1 eV on the basis of DFT calculations carried out on oligomers of increasing length (see Appendix 2), a value very similar to the one reported for the well-known P3HT polymer, [33] and for the TPA moiety. [34] The external contribution, a challenging quantity to calculate exactly, [35] in the majority of the calculations was set to λ e = 0.2 eV, a value expected for a low dielectric constant medium according to classical continuum model calculations [36], and in line with theoretical estimations for oligoacenes crystals. [35] For all molecules of the four MD samples the following conditions were applied: (i) for quantum chemistry calculations hydrogen atoms, absent in the MD model, were added to simulated molecular structures using geometrical criteria; (ii) hole transfer integrals J ij were calculated as one-electron HOMO-HOMO couplings using the dimer projection method [37] at the ZINDO level of theory with the Mataga-Nishimoto potential [38,39]; (iv) periodic boundary conditions were employed to select the pairs of molecules in contact, using 6 Å as the interatomic cutoff distance between atoms belonging to different molecules; (v) site energies i , j were approximated with ZINDO HOMO energies, including the effect of distortion of the molecular geometry in realistic conformations. Although the energetic disorder introduced by intermolecular electrostatic interactions, due to permanent and induced dipoles, can play a key role in improving charge separation at heterojunctions, [40,41,42] and charge-carrier localization [43], it was not included in this work because of the nonpolar nature of TFB, which implies low energetic disorder. as confirmed by Malliaras et al. [17].

Charge Propagation

The propagation of charges is simulated with the Kinetic Monte Carlo (KMC) technique. The charge is initially located at a randomly chosen site i. The hopping rates k ij from donor site i to any potential acceptor site j, computed using equation 5.1, are employed to generate instantaneous hopping times τ ij in the framework of the "first reaction" method [44] : (refer Chapter-02, Theoretical Methodology for details)

τ ij = 1 k ij ln 1 ξ (5.3)
where ξ is a random number drawn from a uniform distribution within the interval 0 and 1. The destination site i is selected from the set of available neighbor sites j as the one having the smallest reaction time. The simulation progresses by following the subsequent hopping steps during charge propagation and the simulated elapsed time is increased by τ ii at each hopping step. All the KMC simulations were run assuming a linear voltage profile across the sample, by varying the magnitude of the electric field from E = 1 • 10 4 V/cm to 30 • 10 4 V/cm. In each simulation, a single charge is propagated, using 3D periodic boundary conditions, until it covers a fixed distance d = 4 µ along the direction of the electric field vector. Since only one charge is present at a time, the simulation targets the limit of very low charge densities, following experimental evidence of charge density-independent mobility for a polymer very similar to TFB. [45] For each charge propagation run k, mobility is then estimated as:

µ k = d E 1 τ k (5.4)
where τ k is the sum of all hopping times τ ii , i.e. the time needed for a charge propagating in the field direction to cover the distance d. Unless explicitly stated, mobility values reported in this chapter are averages over 500 up to 750 KMC runs performed with the electric field vector parallel to x, y and z Cartesian directions, and over the four MD samples. Average values are reported either as arithmetic means, µ = ( k µ k )/N , or as logarithmic means µ ln = exp( k ln µ k /N ), as the latter is a better estimator of the mobility in small samples. [46] It was noticed that for the sample sizes and propagation distances considered, the two means are almost coincident, indicating a satisfactory convergence of the KMC simulations.

Analysis of charge mobilities

According to Bässler [47], at intermediate fields and in conditions of thermally-assisted hopping, the dependence of mobility on electric field modulus can be described by the following empirical relation:

µ(E) = µ 0 exp β √ E (5.5)
where β > 0 is the Poole-Frenkel factor. Using zero field, temperature-dependent mobilities µ 0 (T ) derived from equation 5.5, the "apparent" diagonal energetic disorder σ A and the infinite temperature mobility (µ ∞ ) may in turn be obtained by extrapolating the temperature dependence of µ 0 with:

µ 0 (T ) = µ ∞ exp - 4 9 σ A k B T 2 (5.6)
The Poole-Frenkel factor depends on both apparent diagonal (σ A ) and off-diagonal (Σ A ) disorders:

β(T ) = C σ A k B T 2 -Σ 2 A (5.7)
with C being an empirical proportionality factor. The combination of equations 5.5, 5.6, and 5.7 provides an approximate [48] but universal equation for the dependence of mobility on temperature, field, and energetic disorder. This theoretical framework is used to compare calculated/simulated data with experimental results obtained by Malliaras and coworkers [17,49].

In the original Gaussian disorder lattice model (GDM) [47] σ A and Σ A correspond to the standard deviations of the distributions of site energies (the density of states, DOS) and transfer integrals. In this these, two parameters are labeled as "apparent", because they are effective quantities extracted by fitting the mobility obtained from the off-lattice KMC simulations with the equations above, that are strictly valid for the Gaussian disorder model. The "apparent" tag serves to distinguish the phenomenological disorder parameters obtained by such a fit from the microscopic parameters employed in the KMC simulations. Although the GDM model prescribes a 1/T 2 dependence of the logarithm of the mobility, it has been shown experimentally that often the temperature dependence of the mobility of organic semiconductors at low fields, [50] and for exponential densities of states, [51] can be well reproduced with an Arrhenius law, a behaviour well known also for inorganic disordered semiconductors [48]:

µ 0 (T ) = µ ∞ exp - ∆ k B T (5.8)
The latter equations are therefore used to fit the temperature dependence of the zero field mobility, obtained by extrapolating simulation and experimental data with equation 5.6. The infinite temperature mobility in equation 5.8 is primed because it is different (higher) with respect to the one obtained by fitting with equation 5.6.

Results and Discussions

Charge transport parameters

Initial analysis is focused on the fluctuations of site energies in the samples, which generate the so-called diagonal (energetic) disorder, known to strongly influence charge transport [47,52,53,54,55,43]. The distribution of site energies (HOMO levels) calculated at MD geometries has an approximately Gaussian shape with standard deviation σ = 50.2 meV, as presented in Fig. 5.2, consistent with the assumptions of the Gaussian disorder model. This value is also not far from the effective disorder σ A = 65.9 meV reported by Fong et al. [17], and was hence used in the KMC simulations. Other fundamental parameters in determining the transfer rate, besides the reorganization energy discussed previously, are the magnitude of the electronic couplings between the different sites (the nodes of the charge percolation network) and the average number of pathways linking one node to the others, i.e. the number of neighboring sites with non-zero electronic couplings.

Examining the calculated transfer integrals it is interesting to notice (Fig. 5.3) the very low values of the couplings, with an average value of 0.40 meV, when for instance values of the order of ∼10 meV are found for crystalline and amorphous fullerenes at short intermolecular distances [43]. In addition, as it happens also for crystalline systems along the directions with small couplings [56], the standard deviation σ J = 1.43 meV is higher than the mean value, suggesting an important contribution of dy-namic fluctuations on the mobility values. [57] Despite the low couplings, the number of neighbors with non zero coupling (J > 0.1 meV) ranges between 10 and 20 per oligomer, ensuring a high number of percolation pathways for charge transport. Such couplings most often correspond to TPA-TPA contacts, while FLU-FLU contacts are less frequent because of the presence of the bulky dioctyl substituent on one side of the indeno group (cf Fig. 5.1). The chemical design principle of attaching the solubilizing side chains on the electron acceptor (FLU) is of course beneficial for charge transport, since the hole is mostly localized on the TPA unit (87% according to AM1 calculations), and then the HOMO-HOMO contacts are maximized if this region of the backbone is not encumbered by the alkyl chains.

Mobility calculated for static networks

In the following sections the relation between charge transport simulation parameters and calculated mobility as a function of temperature and electric field are analyzed in detail, utilizing the comparison with available experimental data as a quality test of the simulation predictions. In doing this, it is implicitly assumed that an oligomeric representation of the polymer is sufficient to describe the charge transport for longer chains; this entails that the conjugation length of the charge is assumed to be shorter than the oligomer length, and that inter and intrachain charge transport mechanisms are equally described by oligomer-oligomer hopping and have identical characteristic times. [5] A low conjugation length is expected for an amorphous and relatively flexible polymer like TFB, which was verified at the AM1 level observing that for the four MD structures the extra charge is indeed localized on average within 1.5 units only. A similarly short localization length was reported by Sancho-García et al. for polyfluorene and TFB. [26] Working at the low molecular weight limit implies that trapping at the chain ends is maximized, [7] hence it is expected to slightly underestimate experimental mobilities for long polymer chains, which for TFB are of the order of 10 -2 cm 2 /(Vs) at room temperature. [17] As a starting point, KMC simulations are performed on the original MD samples (denoted as R 1 in Table 5.1) with site energies and transfer integrals calculated for each sample at the ZINDO level ( Q , J Q ). In these simulations, charges are therefore propagating in time in a charge transport network whose properties are static in time, corresponding to single snapshots of the samples. As discussed in the literature, [58,33,59] transfer integrals fluctuate over timescales often shorter than the typical hopping times and taking a static picture eventually leads to an underestimation of materials mobility. The KMC output at different electric fields are reported in Fig. 5.4. Apart the underestimation of the experimental mobilities of one or two orders of magnitude, what is striking in Fig. 5.4 is the absence, at least at the higher temperatures, of the Poole-Frenkel behaviour (ln µ ∝ E 1/2 ) reported experimentally. To demonstrate that this offset is not originated by an unrealistically high value of the reorganization energy, also reported in the same panels is the mobilities calculated with λ=0.1 eV (filled symbols), which can be considered a lower limit for this parameter. The calculated mobilities for λ=0.1 eV are of course larger than the ones obtained at λ=0.3 eV (empty symbols), consistently with equation 5.1, but they are still one order of magnitude lower than the experimental values and show a decreasing trend of mobility with field at high temperature.

Effect of sample size and of dynamic diagonal disorder

In order to understand the factors that contribute to the large offset between simulated and experimental mobilities, KMC simulations were carried out by considering different conditions: (i) sample size, varied by replicating in space the original MD box (R 1 ) two or three times along each Cartesian direction (R 2 and R 3 respectively);

(ii) diagonal disorder, with site energies assigned from quantum chemistry calculations ( Q ), or randomized by extracting site energy values from a Gaussian distribution having the same standard deviation of the calculated HOMO energies, either at the beginning of each KMC simulation (static diagonal disorder, S G ) or at each KMC step (dynamic diagonal disorder, G ); (iii) calculated and fixed (J Q ) or dynamic off-diagonal disorder (J G ). While the fixed case for transfer integrals has an exact correspondence with the fixed site energies, in the dynamic case disorder is added on top of the quantum chemistry values, using random numbers drawn from a Gaussia distribution with σ J =1.43 meV. Finally, (iv) KMC simulations are performed in which off-diagonal disorder, that obtained at T M D =300 K, is assumed to be dynamic and is rescaled to the simulation temperature T using a factor (T /T M D ) 1/2 ( T G , J T G ), according to the electron-phonon coupling model [60,39]. Furthermore, expected in-0 100 200 300 400 500 600 crease of the mean square transfer integrals is considered in an average way by adding the variance of transfer integrals (σ 2 J ) to the absolute value of QM transfer integrals (J Q ) [61,62]. The meaning of the symbols corresponding to different simulations conditions are summarized in table 5.1. HOMO energies drawn from a Gaussian distribution at the beginning of each KMC simulation, with σ =50.2 meV (static diagonal disorder) J Q

E 1/2 (V 1/2 /cm 1/2 )
Electronic couplings fixed at their quantum mechanics value. J G Electronic couplings fixed at their quantum mechanics value, + random fluctuations x extracted at every KMC step from a Gaussian distribution with σ J =1.43 meV : J = J Q + x J T G Electronic couplings augmented according to a temperature-dependent disorder,

J 2 = J 2 Q + σ 2 J • (T /T M D )
The limitations of calculating mobilities from small-sized systems, yielding to unrealistically high values and excessively broad distributions of single charge velocities, are well discussed in the recent literature. [53,27,46,63] Among others, Andrienko and coworkers [53,64] underlined the relationship with the calculated mobility dispersion, energetic disorder and sample size in amorphous systems, and showed that for large disorder values (σ > 3k B T ), only samples containing millions of sites can lead to correct results. The problem originates from the fact that the thermally accessible, low energy tail of the density of states, which becomes populated if the charge propagation experiments are adequately long, is not well represented by a discrete distribution of site energies (Fig. 5.2), in particular if the site energies are fixed in time and the distribution is broad (i.e. the disorder is high). In order to alleviate sample size effects and test their importance as a possible cause of the mismatch with experimental mobilities, (i) a static diagonal disorder is introduced by randomizing the site energies at the beginning of each KMC run, (ii) the simulations are performed for periodic replicas of the original MD boxes, so as to have a larger number of sites and then a smoother density of states, and (iii) the mean logarithmic mobilities were calculated as suggested by Bobbert and coworkers. [46] 0 100 200 300 400 500 600 By comparing the simulated mobilities for replicated samples R 2 and R 3 , it is observed that the differences between mobilities for the two systems is minimal. In fact for TFB, sample size is not expected to have significant impact on mobility due to the small energetic disorder present in these systems (σ < 3k B T ) at every temperature studied. The convergence with system size is confirmed by the similarity between arithmetic and logarithmic averages (points and lines in Fig. 5.5, respectively), and consequently in the remainder of this chapter only logarithmic mobilities are reported. More interestingly, by comparing figures 5.4 and 5.5, it turns out that the static randomization of site energies scales down the mobility by a factor of about five. This reduction can be explained by the removal of (spurious) static correlations between the site energies which is instead present if Q values are used; in fact, as shown for instance by Kordt and Andrienko, [65] uncorrelated disorder generally leads to lower mobilities.

E 1/2 (V 1/2 /cm 1/2 )
To investigate the existence of spatial correlation in R 1 | Q |J Q simulations, a closer look into the mobility anisotropy for one of the four original MD samples (Fig. 5.6, filled symbols) is considered. Quite strikingly, there are marked differences along the three directions of application of the electric field, which in principle are not expected for an amorphous material. This behavior can be attributed to the relatively small size of the simulation samples (with respect to the oligomer size), since a few hopping events are sufficient for a charge to travel across the whole box. The randomization of site energies (grey-shaded and empty symbols) seems to heal the problem and improve the realism of transport landscape by yielding to more isotropic mobility values, albeit two order of magnitudes lower than the experimental ones and decreasing with increasing electric field. The mobilities obtained with frequent randomization of site energies (dynamic diagonal disorder, G ) are as expected slightly larger (and more isotropic) than the ones with static (random) site energies (static diagonal disorder, S G ).

Mobility for dynamic networks

In the discussion so far, it was elicited that despite adopting reasonable approximations and parameters, calculated mobilities for static networks are distant from experimental time-of-flight measured values. One possible source of discrepancy could be the adoption of the semiempirical ZINDO Hamiltonian in the calculation of the electronic couplings. ZINDO is in fact known to underestimate the couplings with respect to DFT calculations, however the scaling factor ranges approximately from one to two depending on the material [66,43], which is not sufficient to recover two orders of magnitude in predicted mobilities. The other possible source of discrepancy with respect to experimental observations is the neglect of time fluctuations of the couplings, i.e. the static description of the percolation network. Cornil and coworkers elegantly showed how including fluctuations in KMC simulations increases the calculated mobility only when the average value of the couplings is smaller than their standard deviation ( J /σ J < 1) [57]. The existence of these fluctuations, essentially generated by low frequency vibrations and with typical correlation times in the subpicosecond range, were observed by many authors [67,68,58,69,33,59].It then appears a sensible choice to introduce dynamic disorder on transfer integrals, as an effective way of considering small fluctuations of atomic positions around their equilibrium value, the so-called positional disorder [47]. The electronic couplings calculated from the four MD samples indeed provide an instantaneous photograph of the network which magnifies the hindering effects on mobility due to trap sites.

In order to improve the realism of the simulation model, the transfer integrals are reshuffled at each KMC step by adding to each ZINDO-calculated coupling a random number extracted from a Gaussian distribution with σ J = 1.43 meV. It is worth noting that this scheme (J G conditions) (i) assumes that the correlation time of transfer integral fluctuations is faster than the average hopping time; (ii) preserves the topology of the network, but allows for a time-dependent probability of the different pathways; (iii) accounts for the distance dependence of the transfer integrals; (iv) presumes that spatial fluctuations are a measure of the fluctuations in time, as it is expected for a crystal, and (v) neglects the temperature dependence of the disorder. A more rigorous way of accounting for dynamic off-diagonal disorder, at least in the framework of Marcus theory, is to augment all the squared transfer integrals of the dynamic variance [61], wherein the off-disorder takes into account the dynamic fluctuation on transfer integrals, and to consider that this variance is proportional to temperature [39] (J T Q conditions).

In Fig. 5.7, the simulation results from KMC simulations performed on replicated sample R 3 are plotted, and compared with the experimental data on the right panel.

Clearly, introducing dynamic fluctuations on the electronic coupling increases significantly the simulated mobility and is key to improve the agreement with experimental data, which confirms the importance of including dynamic fluctuations in charge transport studies of amorphous systems characterized by small values of electronic couplings. ). On the rightmost panel are the mobilities from experimental time-of-flight measurements (adapted from references [17] and [15]). The second set of measurements at 295 K, denoted by a star, is reported for completeness; however it comes measurements where an excessive dispersion of the signal was registered, probably due to the presence of impurities acting as traps. [15] The last step in this systematic review of the KMC simulation parameters concerns the temperature dependence of the energetic and positional disorder. So far an assumption is made implicitly that the disorder is independent on temperature, but it is worth testing the opposite case and how much this would affect the conclusions. It can be seen from Fig. 5.7 that actually the impact of the temperature-dependent disorder is quite weak, confirming that the key factor for recovering the experimental order of magnitude is indeed the introduction of the fluctuations of the electronic couplings, with the S G |J T Q setup giving the best agreement with the experimental results. Moreover, it can be seen how adding dynamic fluctuations on the electronic coupling, which has the physical meaning of incorporating spatial fluctuations in the simulation, inherently present in experiments, the expected Poole-Frenkel behavior is eventually recovered.

E 1/2 (V 1/2 /cm
E 1/2 (V 1/2 /cm 1/2 ) 10 -3

Analysis of field and temperature dependence of mobilities

Analytical equations providing relatively simple formulae to fit the variation of mobility with field and temperature [47], as well as charge density, [70,48] are often used in the experimental literature, with the double purpose of validating theoretical models and extracting intrinsic material parameters under the tacit assumption of the validity of the models. Here the same exercise is performed using simulated mobilities, in order to compare directly with the parameters extracted from experiment, and to verify the internal consistency of the various theories developed to describe charge transport, since some of the microscopic parameters, and notably the magnitude of the energetic disorder, are known beforehand.

In Fig. 5.8 the temperature dependence of zero field mobilities were analyzed utilizing either equation 5.6, which is appropriate for hopping in a Gaussian density of states (left panel), or equation 5.8, typical of an exponential DOS, or of a Gaussian DOS in the presence of extrinsic charges (right panel). [50] Even if the fitting with the Gaussian disorder model is more accurate as expected, the Arrhenius equation gives rather satisfactory results too, although not in the whole temperature range. The two fits produce very different infinite temperature mobilities (µ ∞ µ ∞ in Table 5.2) as a consequence of the different temperature dependence (refer to Eqs. 5.6 and 5.8), however the apparent energetic disorders σ A and the activation energies ∆ are rather consistent, with ∆ 2σ A for all KMC setups and experiments. Changing the reorganization energy from λ = 300 meV to 200 meV does not alter significantly the activation energy, consistent to what was reported in reference [71] for oriented polymers and fields perpendicular to the direction of chain alignment, i.e. inter-chain transport. The variations of a few meV of σ A upon changing the simulation parameters are too small for drawing conclusions about their origin, but it is interesting to notice that all simulations produce similar apparent energetic disorder σ A , in line with the experimental value of 66 meV, and only slightly larger than the intrinsic energetic disorder σ = 50.2 meV. The most significant change with KMC conditions, as shown already in figures 5.4, 5.6, and 5.7, occurs to infinite temperature mobilities µ ∞ ; in particular the best match with the experiments of Malliaras et al. [17] is achieved with the static diagonal disorder and temperature dependence of dynamic off-diagonal disorder (R 3 | S G |J T Q ). Also plotted in Fig. 5.8 are the zero field mobilities as obtained with R 3 | G |J G conditions but with a lower reorganization energy (λ = 0.2, green pentagons), with the double intention of demonstrating the robustness of the results versus the partly arbitrary choice of this parameter, and of showing how the experimental data are well bracketed by the KMC results with λ = 0.2 and λ = 0.3 (red squares).

To conclude this discussion, emphasis is put on the material parameters that are extracted from fitting the temperature dependence of the Poole-Frenkel factors β with Eq. 5.7 (Fig. 5.9). The analysis could be performed only for some simulation setups because not all of them showed a Poole-Frenkel behaviour (Eq. 5.5), and it consists in utilizing previously derived σ A values in Eq. 5.9, obtaining the parameter C from the slope of log β vs (σ A /k B T ) 2 , and finally obtaining the square of the positional disorder Σ A as the intercept of the fitting line with the x-axis. Fig. 5.9 shows that experiment and simulations follow closely the behaviour predicted by Bässler [47], however the angular coefficient obtained from simulations is about half the experimental one (see also Table 5.2). If the diagonal disorder is instead assumed as completely static, indeed it is observed from Fig. 5.9 and Table 5.2 that larger values of C are obtained, very close to experimental results. The positional disorder parameter Σ A , which is a measure of the spread of the transfer integrals, shows instead a more consistent agreement between experiment and simulations, indicating that the different simulation setups, captures rather well the physical nature of the actual system. To provide a birds-eye view of the effect of different simulation conditions on the absolute value of mobilities and the parameters that are subsequently extracted with respect to Table 5.2, simulated mobilities for all the simulation conditions as per Table 5.1 is presented in Fig. 5.10 as a function of applied electric field with temperature set to 295 K.

β × 10 -3 (cm/V) 1/2 exp G |J ( * ) G G |J G S G |J G S G |J T Q Figure 5
.9: Calculated (R 3 supercell) and experimental [17] Poole-Frenkel factors plotted against the square of reduced diagonal disorder. Dashed lines correspond to leastsquare fits with Eq. 5.7. Fitting parameters are reported in Table 5.2. An asterisk * indicates simulation results for λ = 0.2 eV; λ = 0.3 eV otherwise.

Table 5.2: Physical parameters extracted by fitting experimental and simulated mobilities with equations 5.5, 5.6, and 5.7. Mobilities are given in units of cm 2 /(Vs), σ A and ∆ in meV, β and C in (cm/V) 1/2 . Σ A is a dimensionless quantity. a) Quantity determined at T=295 K. b) Simulation results for λ = 0.2 eV; λ = 0.3 eV otherwise.

Quantity exp [17] exp 0 100 200 300 400 500 600

[15] R 1 | Q |J Q R 3 | G |J Q R 3 | S G |J G R 3 | G |J G R 3 | G |J b G R 3 | G |J T Q µ 0 • 10
E 1/2 (V 1/2 /cm 1/2 ) 10 -4 10 -3 10 -2 µ (cm 2 /Vs) exp a exp b G|J ( * ) G G|JG S G |JG S G |J T Q Q|JQ D G |JQ S G |JQ

Conclusions

A theoretical study of the charge transport properties of TFB fluorene-triphenylamine copolymer, a system with remarkably high hole mobility among amorphous organic semiconductors, is presented. KMC simulations based on atomistic input, i.e. atomistic molecular dynamics to build the morphologies and quantum chemistry to extract charge hopping rates, are employed to calculate the temperature and electric-field dependence of the hole mobility. Detailed analysis allowed rationalizing why amorphous semi-conducting systems, characterized by poor intermolecular packing and small charge transfer integrals, can have reasonably high mobilities, highlighting the crucial role of dynamic fluctuations on charge transport. The introduction of dynamic off-diagonal disorder boosts mobility by two orders of magnitude, allowing the achievement of a semi-quantitative agreement with experimental mobility values and the recovery of the so-called Poole-Frenkel behavior of the mobility against the electric field. These results provide a strong evidence of the importance of dynamic energetic disorder for the effective simulation of charge transport in the hopping regime, and introduce a relatively simple and cheap protocol for the computational screening of amorphous organic semiconductors.

CHAPTER 6 ELECTRO-MECHANICAL RESPONSE IN RUBRENE SINGLE CRYSTALS

Stretchable and flexible electronics constitute the core elements for a variety of cuttingedge applications, ranging from mechanical sensors to wearable or biocompatible devices [1,2]. While a controlled response to mechanical deformation is at the basis of the operation of devices such as pressure sensors, flexibility is a key feature for many other applications such as foldable displays or photovoltaic panels, and most notably for all diagnostic devices interfaced with the human body [3,4]. The latest development in the field include electronic skin [5], bio-integrated circuits [4], implants, as well as wearable electronics. Organic polymeric or molecular semiconducting/conducting materials are particularly suited for the fabrication of flexible architectures, owing to their light weight, reduced cost, and the ease of processing; it is thus fundamental to understand the effect of mechanical deformation (strain) and pressure (stress) on their electrical response [6,7]. Remarkably, despite the large number of studies dealing with the performance of organic electronic devices under deformation, to date only a limited number of studies have focused on single crystals and investigated the microscopic relationship between strain and charge transport.

In addition to the insufficiency of specific literature, the peculiarities of each system and the differences in the adopted experimental setups, often make it difficult to draw generalized conclusions [8].

Among crystalline organic semiconductors, rubrene is one of the most studied, not only because of its record hole mobilities (reaching several tens of cm 2 /V s) but also due to the relative ease of growing large, defect-free single crystals [9,10], qualities that make this material particularly suitable for proof-of-principle experiments [11].

A fundamental aspect of the research on rubrene-based devices has regarded their performance reproducibility in terms of their electrical and structural stability. However, although the impact of many factors (i.e. humidity, light, oxygen) on rubrene field effect transistors (FETs) was investigated to a great extent [9,11], the influence of mechanical strain received only limited attention. In a first study, Briseno et al.

showed that the performances of flexible rubrene FETs upon mild bending could be superior or comparable to more conventional rigid devices [12]. More recently, Briseno and co-workers carried out a deeper investigation by exploiting wrinkling instability to apply local strains of different magnitudes along the conducting channel of the transistor. They showed that mechanical compression along the a axis (the principal π-stacking direction, see : Fig 6 .1-[a]) correlates with an increase in mobility, while expansion diminishes it [13]. In another complementary study, Batlogg and coworkers [14] measured threshold voltage and mobility changes along a axis upon application of strain along the two in-plane directions a and b, pointing out that the effect of strain along a and b is similar with respect to measured mobility, contrary to physical intuition that a positive ab Poisson ratio indicates a contraction of the lattice along a axis when b axis is expanded, and then an increase of mobility along a axis is expected, and vice versa. Available studies on rubrene indeed reveal the complexity of the relationship between strain and charge carrier mobility: the simple association of compressive (tensile) strain to a decrease (increase) of the intermolecular distance and of their fluctuations [15], and in turn to the rise (decrease) of the electronic couplings between the molecules, does not seem sufficient to rationalize the experimental results. For instance, the application of hydrostatic pressure is generally believed to increase mobility [6,16,17,18], however in some cases an opposite effect is measured [19]. Compressive strain or bending typically brings an increase of mobility for oligoacenes, but the effect is not general and negative responses to strain have been reported [20,21,22]. In addition, the sensitivity changes from material to material: it has been shown that the mobility of TIPS-pentacene hardly changes upon bending [23,24], while the expected behaviour is achieved upon lattice deformation [22,25] or shear [26]. To complicate further the picture, both the sensitivity and the maximum measurable strain after which irreversible structural changes occur depend not only on the semiconducting material but also on the underlying substrate [6,27,28].

Indeed the complexity of the problem suggests that, for the time being, only punctual studies on a given material can provide the desired structure (strain)-property (mobility) relationship. Here, such an attempt is described for rubrene, where modern computational techniques [26,29] are coupled to a novel experimental setup. Since previous experimental studies on pentacene demonstrated that not only the variation of cell parameters but also grain morphology can affect the response to strain [30,31], emphasis is put on single crystal [32] devices and simulations as a way to eliminate the latter effect. A multiscale simulation protocol is used to predict both the mechanical properties of rubrene from stress-strain curves and the corresponding variations of charge carrier mobility. Theoretical results are compared to experimental values obtained using field effect transistors (FET) produced by lamination of single crystals on top of a cantilever beam [33,34,35] (cantiFET), and recording their electrical response to strain in real time.

Simulating mechanical strain

Molecular dynamics (MD) simulations were performed with the Generalized Amber Force Field (GAFF) and the NAMD [36] software. The dihedral angle between the tetracene and phenyl moieties was parameterized against DFT calculations at PBE0/def2-TZVP level carried out on phenyltetracene using Orca 3.03 [37]. Molecular Dynamics (MD) simulations were carried on a rubrene supercell that exhibits orthorhombic phase, consisting of 256 molecules with periodic boundary conditions applied along the three crystallographic directions a, b and c. MD simulations were performed at room temperature and constant volume. The force field was then validated simulating a bulk rubrene crystal supercell. The initial coordinates and lattice parameters were taken from ref. [10]. A first MD simulation run in the NPT ensemble was carried out for 10ns at 300K, followed by relaxation of 20ns in the NVT ensemble.

The average values extracted from the simulation are reported in Table 6.1. Subsequently, to characterize the mechanical response to applied strain of a bulk rubrene crystal, two sets of simulations were carried out at constant volume and temperature set to 300 K. These two simulations correspond to uniaxial strain or uniaxial stress conditions (see Fig. 6.1 and Fig. 6.3), wherein the crystallographic axes of the supercell were selectively rescaled with respect to their zero strain average values at atmospheric pressure.

Uniaxial Strain

For the uniaxial strain condition, the supercell and the intermolecular distances were rescaled in order to compress/expand the cell along a, b or c axis independently, which entails that only one lattice constant of the crystal was rescaled, while the other two were kept constant. Uniaxial strains to the crystal lattice were imposed by simultaneously modifying one of the supercell box sides and the intermolecular distances along that direction while leaving the other ones unaltered, with values ranging from -0.4% to +0.4% in steps of 0.1%. After equilibration, for each value of applied strain along every direction (crystallographic axis), a 10 ns-long production run was carried out at constant volume, while recording the pressure profile along the three sides of the box. A schematic representation of uniaxial strain conditions along with the resulting crystal lattice parameters for strains applied along a and b axes is presented in Fig. 6.1.

[a]

[b] For any material, stress (σ) and strain ( ) are related by the stiffness tensor, given by Eq. 6.1, where i, j, k, l run over Cartesian axes, which coincides with crystallographic axis in the present case.

[c] [d]
σ ik = k l C ijkl kl (6.1)
In the absence of shear strain, and by adopting Voigt notation, it is possible to reduce the stiffness tensor to a symmetric 3×3 matrix where σ ii = k C ik kk , so that only six independent elements can be identified. Then, the three corresponding stress-strain equations for the a, b and c directions, for uniaxial strain condition become,

  σ aa σ bb σ cc   = C   aa 0 0   ;   σ aa σ bb σ cc   = C   0 bb 0   ;   σ aa σ bb σ cc   = C   0 0 cc  
The stiffness tensor elements were computed by fitting each stress-uniaxial strain plot (Fig 6 .2), with a line having intercept at (0, 0). The angular coefficients of the stress calculated along the directions of strain correspond to the three diagonal elements C ii of the tensor. Off-diagonal elements were derived in a similar fashion, but taking into account the symmetry of the tensor by averaging C ij and C ji . For instance, C ba was calculated as σ bb / aa from the simulation where the strain was applied along a axis, and correspondingly C ab = σ aa / bb when strain was applied along b. From uniaxial strain simulations the stiffness tensor of crystalline rubrene was calculated, The order of magnitude of the values is in line with bulk and Young moduli reported for crystalline oligoacenes [18,38]. From the matrix elements of C, the reduced elastic coefficients for in-plane stress

C =   16
C ii = C ii -C 2 ic /C cc (with i = a, b
) and the anisotropy ratio C aa /C bb were also calculated and compared with recent measurements and calculations [39,40]. The values reported in Table 6.2 show a broad agreement with published theoretical and experimental values, substantiating the reliability of the predictions made through our computational protocol. It is to be noted that, calculations from this work, as well as the ones in reference [39] and the measurements in reference [40], indicate two negative Poissons ratios, calculated as follows, and reported in Table 6.3. The higher magnitude of ν ab and ν ba clearly denotes the high inter-dependency in the strain behavior existing within the ab crystal plane. 6.2: Comparison between reported reduced elastic constants (in units of GPa) and anisotropy ratios for rubrene single crystals, from experiments and simulations. Note that a different notation for crystal axes was used in references [39,40] with respect to this work . Term Exp [39] Exp [40] Sim [39] Sim [40]

C -1 =       1 Ea -ν ba E b -νca Ea -ν ab Ea 1 E b -ν cb E b -νac Ec -ν bc Ec 1 Ea       Table

Uniaxial Stress

Uniaxial strain condition represents one limiting case of a crystal laterally confined, clamped, or any other experimental setup preventing a molecular rearrangement along the cell vectors normal to strain. In a second series of simulations, the other limiting case was investigated, that of uniaxial stress, which ideally corresponds to a rubrene crystal free to respond to applied stress with a lateral deformation according to its Poisson ratio elements. This was done, by rescaling the crystallographic axis directions for a strain applied in one direction using Poisson ratio elements. The following equations are employed, using the calculated stiffness tensor elements to remove the residual off-diagonal strain, thus achieving uniaxial stress conditions, where σ ii values were chosen to match the corresponding ii obtained from uniaxial strain conditions. A schematic representation of uniaxial stress conditions along with the resulting crystal lattice parameters for strains applied along a and b axes is presented in Fig. 6.3 and the stress-strain curves, calculated for uniaxial stress conditions are presented in The main difference when comparing the variation of the in-plane intermolecular distances (|r a |, |r b |) for uniaxial strain and uniaxial stress conditions from Figs. 6.1 and 6.3, respectively, is the variation of the components of the intermolecular distances orthogonal to the strain direction. While in the case of uniaxial strain such terms show little to no dependence on strain, for uniaxial stress especially along a and b crystallographic directions, the stress along one axis affects all intermolecular distances. Conversely the variation of |r c | is quite modest in the considered interval when stress is applied along a or b axis in line with the low magnitude of the corresponding stiffness matrix components C 13 and C 23 . Further, the linear variations of all the crystalline parameters (plotted quantities) confirm that the chosen interval of ±0.4% strain is located within the elastic regime. The rescaled cell parameters for uniaxial stress conditions along with respective stress values is reported in Table . 6.4.

Electro-Mechanical coupling in rubrene single crystals

Charge transport in rubrene single crystals is described to be band-like at low temperature, and becomes dominated by activated hopping at higher temperatures [11,41,42]. Further, charge carrier mobility in Rubrene is strongly anisotropic [43,44,45,46], as can be surmised from the observation of its crystal structure (refer Fig. 6.1/6.3). Molecular arrangement in the ab plane is herringbone like and highest electronic coupling (J a ≈ 95 meV [41,45]) is found along the a-axis, which corresponds also to the π stacking direction and to the facet of the fastest crystal growth [10]. J b , correspond to the electronic coupling between the molecules in the herringbone like packing along b-axis, is considerably lower, around 10 meV . J c is the transfer integral of the two closest out-of-plane neighbours (C-axis, neighbors perpendicular to ab plane) is of the order of 1 meV [3].

In order to evaluate the electro-mechanical response in rubrene single crystals, dimers were extracted from MD morphologies of Rubrene simulated super cell, for both uniaxial strain and uniaxial stress conditions, as a function of strain, and transfer integrals are evaluated. Two different kinds of dimers configurations were extracted, (i) the average coordinates of all dimers in the supercell were averaged every 100 ps over a MD production run of 10 ns, (ii) for a total of 34 -300 ps-spaced from 10 ns simulation time, MD morphologies were extracted and dimers were extracted from each morphology, to account for instantaneous fluctuation in atomic positions and molecular conformational space that result from thermal fluctuations [47,48]. The two cases would henceforth be referred to as "average dimers" and instantaneous dimers", respectively.

For all the dimers, the transfer integrals were calculated in the one electron approximation as J = φ A HOM O | H|φ B HOM O , where A and B represent the two adjacent rubrene molecules that correspond to a dimer involved in the hole transfer, using the projective method at the ZINDO level of theory [29,49,50] (see Chapter 2 for details).

Distribution of transfer integrals

Distribution of transfer integrals was computed for instantaneous dimers, for strains applied along a, b and c crystallographic directions. The corresponding instantaneous transfer integrals along the specified directions (a, b, & c) are henceforth called as J t a , J t b and J t c , respectively, so as distinguish them from the electronic couplings calculated from average dimers, which are henceforth referred to as J a , J b and J c , respectively. The distributions of the instantaneous electronic couplings at strain values of 0, +0.4% and -0.4% (extreme strains considered in this study) along a, b and c directions are reported for MD simulations performed in both uniaxial strain and uniaxial stress in Fig. 6.5 and 6.6 respectively.

For both uniaxial strain and uniaxial stress conditions, the distributions J t a and J t b show a skewed normal distribution (or log-normal distribution) whereas J t c looks more like a folded normal distribution which is also skewed, whose standard deviation: σ(J t a ), σ(J t b ) and σ(J t c ), remains mostly unchanged upon positive or negative strains. Changes in the mean values of the distribution of instantaneous transfer integrals: < J t a >, < J t b > and < J t c >, as a function of applied strain are similar to the transfer integral values obtained from average dimers: < J a >, < J b > and < J c >, therein reflecting the relative mobility trends as a function of variation of transfer integrals upon strain being similar, irrespective of the choice of instantaneous or average dimers. The distribution of instantaneous transfer integrals with respect to the mean value, standard deviation (σ), skewness and kurtosis and η =< J t > /σ for both uniaxial strain and uniaxial stress conditions are reported in Appendix 3: Tables A3-1,A3-2 and A3-3 for uniaxial strain conditions and in Tables A3-4, A3-5, A3-6 for uniaxial stress conditions. Skweness of the distribution is a measure of asymmetry while kurtosis is the measure of 'tailedness' of the distribution, which entails that both of these values describe the shape of the distribution. Also reported in the tables are the values obtained from average dimers : < J a >, < J b > and < J c >, to exhibit the fact that the trend for both type of configurations (average and instantaneous dimers) is similar, despite instantaneous averages being (considerably) more noisy. In the present case, the major role of the fluctuations of intermolecular distances (recently suggested for explaining the increase of charge carrier mobility of a naphthobenzodithiophene upon compressive strain [15]) can be ruled out for both uniaxial strain and uniaxial stress conditions, since the standard deviations of |r a |, |r b | and |r c | vary less than 1% in the considered strain interval (See Fig. A3-1). To fully understand the effect of the intermolecular fluctuations on the modulation of the transfer integrals, it is useful to observe the values of the parameter η i = J i /σ i , in absence of strain, where σ i is the standard deviation of the J i distribution, and i = a, b, c. The smaller the value of η, the higher the impact of lattice dynamics on the transfer integral value [47,51,52]. The three rubrene nearest neighbors are characterized by three distinct values of η: from the calculations which are 3.4, 2.0 and 1.6 for dimers along a, b and c axis, respectively. All the values of η are large enough to indicate a weak or negligible dependence of the mobility on the amplitude of intermolecular vibrations, hence it is not expected that a strain-induced variation of σ i could significantly affect mobility values. In addition, imposing compressive/tensile strains on the rubrene crystal does not significantly modify the value of σ i . As highlighted in Figs. 6.5 and 6.6 (also see Appendix 3, Tables A3-1,A3-2, A3-3 and Tables A3-4, A3-5, A3-6), the variation of the average values of J i actually originates from a shift of the whole distribution rather than a broadening. Further, in the attempt of rationalizing the strain-induced variation of the electronic couplings J between the nearest neighbors along the cell vectors, relevant intermolec-ular degrees of freedom were characterized as a function of strain for 1000 snapshots taken from MD simulations along a, b or c axes, in both uniaxial strain and uniaxial stress conditions. Parameters like lateral shift of one molecule along one of its three axes, the rotation of one molecule with respect to its three symmetry axes, etc., were analyzed. However, there is not a specific intermolecular mode that can be held solely responsible for the observed variation of mobility with strain. Details of the analysis are presented in Appendix 3.

Relative mobility

The relative magnitude of mobility along the crystal axis i = a, b, c is calculated at zero field and neglecting fluctuations of the ionization potential, by simultaneously taking into account the variations of transfer integrals and intermolecular distances vectors with strain [45,53], employing the transfer integrals obtained from "average dimers". Mobility, in such condition is given by Eq. 6.2.

µ i = µ 0 k J 2 k [ r k • e i ] 2 (6.2)
Where the summation runs over three types of neighbors : two neighbors along the a axis (π stack direction, refer Fig. 6.1/6.3 -[c]), four in the ab plane (herringbone like packing, refer Fig. 6.1/6.3 -[d]) and four neighbors along the c axis, that are perpendicular to both π stacking direction and ab plane, with intermolecular distance vector r k , e i is the cell axis unit vector and µ o is a constant. These ten neighbors are considered to compute the relative mobility. Transfer integral value squared in Eq. 6.2 is consistent both with a Marcus-like localized hopping regime and with a partially delocalized picture, more appropriate for the intralayer transport in the ab plane [53,54]. The relative changes in mobility, as a function of applied strain, is calculated using Eq. 6.3, where i = a, b, c and strain applied along j = a, b, c, with respect to mobility computed at zero strain (µ 0 i at jj = 0).

∆µ j i = µ j i -µ 0 i µ 0 i (6.3)
Relative mobilities for both uniaxial strain and uniaxial stress (from Eq. 6.3) are presented in Fig. 6.7. A linear regression fit is employed to fit the ratio between the mobility in the presence and absence of strain and the results of each strain-mobility calculation are condensed into a single empirical parameter m ij [6], given by Eq. 6.4, where i and j indicate the direction of the electrical measurement and the direction of applied strain. Values of m ij are reported in Table. 6.5. Uniaxial strain and uniaxial stress conditions give rise to striking differences, albeit in both cases the response is linear in the considered strain interval, and antisymmetric for tensile or compressive strains. Considering the relative variation of mobility measured along the higher mobility axis, a, arising from the strain along the same axis, when uniaxial strain is applied, ∆µ a a decreases by about 5% for a tensile strain of 0.4%, which can be rationalized by taking into account the change in the intermolecular distances r a , since the electronic coupling (here J a ) is known to decrease exponentially with the increase of the intermolecular distance [25,49]. Unexpectedly, when uniaxial stress conditions are applied (orange circles in Fig. 6.7 -[aa]), a marginal increase in ∆µ a a is noticed for tensile strain, yielding to a positive m ij coefficient which cannot be related to the changes of the intermolecular distances along a axis, as their average values along this direction are unchanged with respect to the uniaxial strain case -see green lines in Fig. 6.1-(c) and Fig. 6.3-(c) for comparison.

µ j i µ 0 i = 1 + m ij ij (6.4)
The trend is not repeated when the strain is applied along b axis (Fig. 6.7 -[ab]). The two types of deformation have similar effect on ∆µ b a , with the expected decreasing relative mobility profile at increasing tensile strain (positive values of m ij ). For the mobility probed along the second in-plane direction (b axis), for deformations/strains applied along a and b axis, (∆µ a b & ∆µ b b ) strain-mobility profile is qualitatively the same as the strain-mobility along a axis -see 6. The variation of mobility along c axis (four neighbors that are perpendicular to ab plane), appears to be the most sensitive when compressive/tensile strain along c is applied. ∆µ c c in this case reaches 60% (Fig. 6.7 -[cc]) and follows the expected decreasing trend with strain, as reflected by the large m ij values. Although ∆µ c c 60%, the absolute values of mobility along the c axis is much lower in comparison to those along a axis (π stack direction 95 meV) or along the b axis (ab plane 10 meV). In case of deformation along a and b axes, the mobility measured along c direction (Fig. 6.7 -[ca],[cb]), is less perturbed (in comparison to Fig. 6.7 -[cc]), and in particular it hardly changes upon a stress applied along a. This can be related to the corresponding Poisson ratios: since ν cb and ν ca are relatively small, deformations along a and b axes are not expected to produce large effects on the arrangement of rubrene molecules along the c axis. The effect of strain applied along c axis on the mobility along b axis (Fig. 6.7 -[bc]) shows the expected decreasing trend for tensile strains with a response very similar to the one where stress is applied along b axis (Fig. 6.7 -[bb]). Conversely, for mobilities along a axis, with tensile strains applied along c axis, an increase in relative mobility is noticed (∆µ c a ), in particular for uniaxial stress conditions, and the same trend is also noticed for mobility measured along c upon stress application along a (∆µ a c ) : Fig. 6.7 -[ac] and [ca] respectively, in record with the negative Poissons ratios : ν ac and ν ca .

Experimental Measurements

The electromechanical response of rubrene single crystals under mechanical strain was electrically characterized, in a transistor configuration, and the relative variations of drain-current with strain were compared to the calculated mobility changes discussed above. Different approaches exist to characterize the electrical properties of single crystals under stress and strain. Here a new method based on an Organic Field Effect Transistor (OFET) embedded in a suspended organic micro-cantilever was applied, henceforth called cantiFET, first described by Rao et al [33] and improved recently by Thuau et al., who fabricated highly sensitive electro-mechanical transducers using piezoelectrically-gated OFETs [35]. In the present work, single crystal rubrene FETs were integrated in organic MEMS plastic cantilevers. By bending the flexible cantilever, controlled and uniform strains can be applied to the embedded transistors, thus enabling the accurate evaluation of the electromechanical properties of organic semiconductors (Fig. 6.8). The triangular shape was chosen to induce a uniform longitudinal stress in the crystal when a force was applied at the free-end of the cantilever [55]. The rubrene crystal was positioned near the clamping region of the cantilever, in order to minimize bending-induced strain along the cantilever width direction. Then the applied force at the cantilever free-end induces uniaxial tensile strain in rubrene crystals along the cantilever length direction since in the width direction it is almost clamped. As a consequence, the intermolecular distance along a specific axis on rubrene single crystals is changed in a controllable fashion, as schematized in Fig. 6.8 -[a]. While monitoring the drain current of the transistor for different applied strain values, the electromechanical response of rubrene is characterized accurately and reversibly (Fig. 6.8-[d]).

In practice, different mechanical and electrical configurations were tested experimentally. Fig. 6.9 shows the results obtained for each device configuration: for all setups, the drain current, directly proportional to the mobility, decreases when increasing tensile strain, as predicted by simulations in uniaxial strain conditions (refer to Fig. 6.7 and Table 6.5 to compare with uniaxial stress conditions). The experimental results on the one hand validate the trends obtained by molecular dynamic simulations and on the other hand demonstrate the impact of strain on the mobility of rubrene using an original approach based on OFET-embedded MEMS configurations. Experimental mobility with respect to calculated mobility for uniaxial strain conditions, is presented in Table. 6.6. Table 6.6: Parameter m ij obtained from the linear fit of calculated (uniaxial strain) and experimental mobilities versus strain with

µ j i µ 0 i = 1 + m ij ij m ij MD uniaxial strain Experimental aa bb aa bb µ a -12.4 -16.8 -6 -21 µ b -4.2 -19.2 -11 - 25
From the values of m ij , a good qualitative agreement could be observed for µ a and µ b between experimental and simulated results, with respect to sign of the slope for m ij . Relative variation of mobility along b axis is higher with respect to variation along a axis. This qualitative agreement of simulated mobility variation for uniaxial strain conditions with experimental results, with respect to those for uniaxial stress conditions, indicate that the experimental measurements are influenced by the mechanical adhesion between rubrene single crystals and underlying substrate. However, as the simulated values does not take into account experimental factors like the influence of substrate stiffness and mismatch between the Poisson ratios of rubrene and underlying substrate, the quantitative discrepancy between the values of m ij for simulated and experimental conditions can be attributed to these effects.

Conclusion

Controlled uniaxial tensile strain is applied to rubrene single crystals, simultaneously measuring the induced variation of current, using an original setup based on an organic field effect transistor embedded in a suspended polymeric micro-cantilever [35] (Fig. 6.

-[a] & [b]

). Considering its reduced cost and ease of fabrication, this setup represents an excellent mechanical platform to characterise accurately the electrical properties of organic materials. Contrary to previous reports of an elastic-to-plastic transition occurring via thermal expansion on a poly(dimethylsiloxane) substrate at a strain of 0.05% [27], the registered response is elastic and reversible in the strain range considered (up to 0.18%, Fig. 6.8-[c & d]). Quantitatively, the relative variation of current is negative, coherently with recent experiments exploiting wrinkling instabilities [12], and more pronounced when measured along the b axis (Fig. 6.9 -

[b] & [d]
) than along a-axis (the direction ı stacking and of higher mobility) (Fig. 6.9 -[a] & [c]). Interestingly, the effect is qualitatively independent on the direction of the application of the strain, revealing a strong coupling between the two in-plane directions of rubrene crystals and a microscopic mechanism more elaborate than the mere reduction of the electronic couplings caused by the increase of the intermolecular distances along the direction of application of strain. Rubrene crystals then appear to be more sensitive to strain when compared to TIPS-pentacene [56] or alkylthiophenes [15], as their higher electromechanical sensitivity and stiffness suggest.

The variation of charge carrier mobility as a function of compressive and tensile strain is also predicted by calculations combining molecular dynamics, quantum chemistry, and charge diffusion models. The computational results are validated through the comparison of the stiffness tensor with previously published experimental and theoretical data [39,40], notably indicating negative Poissons ratios when stress and strain couple the a and c (out-of plane) crystal axes. The relative mobility changes under uniaxial strain conditions (Fig. 6.7, blue lines) match semi-quantitatively with measurements, with larger electromechanical sensitivities when strain is applied along the b-axis ( 20) with respect to a-axis ( 10), independently on the direction of the measurement of the current (see Table . 6.5).

The calculations extend the scope of the OFET experiments by considering also the condition of uniaxial stress, corresponding to an ideal case in which the two dimensions of the crystal perpendicular to the applied stress direction are let free to adapt so as to minimize lateral stresses according to the corresponding Poissons ratios (Fig.

6.1-[b]

). While the predicted change in mobility is very similar for uniaxial stress and uniaxial strain applied along b (Fig. 6.7 -[ab] & [bb]), quite surprisingly the response changes of sign for uniaxial stress along a, i.e. the mobility increases for tensile strains and decreases for compressive strains (Fig. 6.7 -[aa] & [ba]), despite intermolecular distances behave oppositely. This result underlines the importance of considering the exact mechanical setup in interpretation of electro-mechanical experiments, because i) the effect of mechanical deformation might not always be readily translated into a predictable mobility change, and ii) the actual conditions of the experiment, such as the substrate mechanical properties, adhesive behavior, or clamping, might strongly affect the magnitude and sign of the electric response of the device.

Concerning the microscopic nature of the electrical response, for rubrene it originates from the shift of the whole electronic coupling distribution upon strain (refer Fig: 6.5), and not from a change in their standard deviation caused by the suppression or enhancement of intermolecular vibrations as recently put forward for explaining strain-mobility trends for a benzodithiophene derivative [15]. Despite the strong dependence of the electronic coupling on the intermolecular distances, their variation alone is sufficient to rationalize the electric response only for uniaxial strain conditions (with mobility measurements performed along the direction of applied stress), and not for uniaxial stress ones or other directions of measurements. In the latter case, although the unambiguous variation of mobility is clearly generated by the rearrangement of the molecules upon strain, it appears to be the complex result of tiny yet collective variations of several inter/intra-molecular parameters that influence the electronic couplings and ultimately, the measured current.

[55] Hussein Nesser, Cédric Ayela, Isabelle Dufour, and Hélène Debéda. Highly deformable printed organic trapezoidal micro-beams for vibration energy harvesting. Flexible and Printed Electronics, 2(1):015001, mar 2017.

[56] B. Blülle, A. Troisi, R. Häusermann, and B. Batlogg. Charge transport perpendicular to the high mobility plane in organic crystals: Bandlike temperature dependence maintained despite hundredfold anisotropy. Physical Review B, 93

:035205, jan 2016.

CHAPTER 7

STRAIN-MOBILITY TRENDS IN CRYSTALLINE ORGANIC SEMICONDUCTORS: A PERSPECTIVE FROM BAND DISPERSION & TRANSFER INTEGRALS

In the previous chapter (Chapter 6), strain-mobility trends of orthorhombic crystalline rubrene were obtained by using Molecular Dynamics (MD) simulations to induce mechanical strain and electro-mechanical responses were computed as a function of variations of transfer integrals and intermolecular distances. The present chapter is aimed at employing an alternative computational protocol with respect to MD, wherein, plane-wave density-functional-theory (PW-DFT) calculations are employed to report the strain induced relative mobility variations as a function of the effective masses. Polymorphs of Rubrene (triclinic and monoclinic) as well as derivatives of [1]benzothieno [3,2-b]benzothiophene (BTBT) family (BTBT and C 8 -BTBT), are considered for this study, to evaluate the impact of the structural anisotropy on the resulting electro-mechanical response. Furthermore, transfer integrals are computed for the strain induced structures and relative mobility variation is also presented as a function the transfer integrals, therein providing a perspective for variation of mobility in crystalline organic semiconductors from both band and hopping transport regimes.

Introduction

The effect of anisotropy of the electro-mechanical response of orthorhombic rubrene, originating from the inherent anisotropy of molecular packing, has been studied in the previous chapter. However, two additional polymorphs were reported for crystalline rubrene, along with the widely studied orthorhombic phase [1,2,3,4,5,6]: triclinic [7,8] and monoclinic [9,10] phases. While the base-centered orthorhombic phase of rubrene is obtained from physical vapor deposition method, triclinic and monoclinic phases are obtained from precipitation or reprecipitation methods [9]. Further, it has been demonstrated that the rate of precipitation and type of solvent employed modulate the shape of micro-crystals of triclinic and monoclinic rubrene structures [9,11].

The different polymorphs of rubrene along with the π stacked arrangement is shown in Fig. 7.1 and the crystallographic parameters are presented in Table 7.1. Both triclinic and monoclinic structures of rubrene differ in molecular packing, relative orientation of the molecules and π-stacked arrangement with respect to each other and also with respect to the orthorhombic phase. The triclinic phase of rubrene exhibits a face-to-face slip stack arrangement between the two neighboring tetracene backbones that constitutes the π-stacking direction, with an inter-molecular distance between the backbone similar to that of the orthorhombic phase [8]. However, the mobilities of triclinic rubrene were reported to be lower than that of orthorhombic one [8]. This decrease in mobility is attributed to the absence of herringbone disposition of the molecules and to the decreased density of molecular packing along the axis perpendicular to π stacking direction [7,8,9]. In contrast to the triclinic and orthorhombic phases, monoclinic rubrene exhibits minimal π stacked interactions, leading to a further decrease in mobility [9]. This entails that, mobility of orthorhombic rubrene is greater than that of triclinic, which in turn is greater than monoclinic rubrene. Table 7.1: Crystallographic parameters for rubrene polymorphs. Space group is represented in Hermann-Mauguin notation and Z corresponds to the number of formula units per unit cell.

Polymorph

Cell paramaters Reference - Although rubrene is one of the most widely studied crystalline organic semiconductors, thienoacene-based derivatives are particularly interesting, since they show high mobility with relatively high air stability [12], due to the highly delocalized electronic structure and low-lying highest occupied molecular orbital (HOMO). In addition to this, strong non-bonded interactions between thiol groups (S-S) and intermolecular π-π interactions in the solid state (crystalline material) promote large intermolecular orbital overlap [13]. Of particular relevance in the thienoacene family is [1]benzothieno [3,2-b][1]benzothiophene (BTBT), which in turn has been at the core for several derivatives [14,15], mainly driven by the objective to find best p-type organic semiconductors from the family of BTBT [16,17]. Takimiya et al had reported the synthesis of series of n(2,7-dialkyl)-BTBT derivatives (C n -BTBT) where n ranges from 5 to 14, wherein the solubility of these derivatives in chloroform at room temperature increases as a function of n, thereby easing the solution processing to obtain a b c crystalline thin films [16]. Subsequently these BTBT derivatives were incorporated into organic thin film transistors and an increase in mobilities recorded as a function of n (C n -BT BT ) was observed [18]. Lemaur et al reported a series of computational studies, wherein the transfer integrals for the family of BTBT and C n -BTBT were calculated and KMC simulations were performed to obtain the in-plane hole mobilities. An increase in calculated mobilities as a function of n in C n -BTBT derivatives is reported [19,20,21], in line with the experimental observations from Takimiya et al [18]. The family of C n -BT BT derivatives for n being equal to zero (BTBT), eight (C 8 -BTBT) and twelve (C 12 -BTBT), with their crystallographic arrangement is presented in Fig. 7.2, along with the crystallographic parameters in Table . 7.2. For BTBT, the π stacked arrangement is along b-axis which correspond to the maximum mobility direction and for C 8 -BTBT and C 12 -BTBT the this arrangement is along a-axis, corresponding, to the maximum mobility direction [19,20,21].

a ( Å) b ( Å) c ( Å) α • β • γ • Volume ( Å3 )
- a ( Å) b ( Å) c ( Å) α • β • γ • Volume ( Å3 ) Symmetry Z - BTBT 8

Computational Details

PW-DFT calculations were performed on rubrene (triclinic and monoclinic phases), as well as BTBT and C 8 -BTBT derivatives at vdw-df-c09 [23] level of theory employing the Quantum Espresso software [24]. Geometry relaxation (cell optimization) was performed starting from the experimental crystallographic structures, as reported in Tables 7.1 and 7.2, employing a k-point sampling grid of 3 × 3 × 3. The kinetic energy and charge density cut-of are set to 693 eV (50 Ry) and 4854.5 eV (350 Ry) respectively. Force, energy and self-consistent-field (SCF) convergence thresholds were set to 5.0 × 10 -4 , 1.0 × 10 -6 and 1.0 × 10 -8 , respectively. Subsequently, compressive and tensile strains were applied along the three crystallographic directions on the optimized structures (up to 0.8% in increments of 0.2%) by rescaling the crystallographic axes. This procedure is identical to the uniaxial strain case studied for orthorhombic rubrene in the previous chapter. For every strain condition (compressive and tensile) all internal degrees of freedom were relaxed self-consistently. Subsequently, band structure calculations are performed at the same level of theory as that of SCF calculations (cell/ion relaxation -PBE/vdw-df-c09). A Gaussian smearing of 27.2 meV (0.002Ry) was employed, similar to the computational procedure employed by Zhang et al [25].

Geometries of the optimized structures (zero strain and strained induced structures) were used to extract the dimers along the principle crystallographic directions (similar to Chapter 6) and transfer integrals were computed at PBE/GGA/DZ level of theory employing ADF package [26] using the projection method. Effective mass is calculated by employing a second order polynomial fit to the band dispersion, from the Γ point to the direction of interest (cf. to Figs. 7.4 and 7.5 for details). In Chapter 6 it has been demonstrated that transfer integrals vary as a function of applied strain reflecting the crystalline anisotropy. Accordingly, the relative change in mobility as a function of applied uniaxial strain is reported in two ways: (i) change of mobility as a function of transfer integrals (J), where in (mobility) µ ∝ J 2 , which holds valid under the hopping-like picture of charge transport (Marcus/Marcus-Levich-Jortner formalisms, see Chapter 2 for details) and (ii) change of mobility as a function of effective mass (m * ), where µ ∝ 1/m * , which holds valid within the Drude approximation (see Eq. 7.2 for details) of band-like transport.

Results & Discussion

Crystallographic parameters, obtained from PW-DFT calculations of the optimized crystalline structures of rubrene polymorphs as well as BTBT derivatives are reported in Table . 7.3.

Table 7.3: PW-DFT optimized structural parameters of rubrene polymorphs (O:orthorhombic, T:triclinic and M:monoclinic) and BTBT derivatives with respect to experimental parameters. ∆V correspond to the percentage change in volume for the PW-DFT optimized structure with respect to the experimental data. 2). Furthermore, earlier studies show that a decrease in intermolecular distance along the π stack direction and/or along the direction transverse to the π stacking is possible [27,18], attributing this to the intermolecular van der Waals interactions along the alkyl group, called the "zipper effect" [28] bringing a tight molecular packing in the solid state.

Type Experimental Theoretical ∆V - a ( Å) b ( Å) c ( Å) α • β • γ • a ( Å) b ( Å) c ( Å) α • β • γ • - Rubrene: O 26.

Band Structure

Band dispersion was calculated for all rubrene polymorphs (orthorhombic, triclinic and monoclinic), and for BTBT derivatives (BTBT, C 8 -BTBT, and C 12 -BTBT). However, band dispersion for the orthorhombic rubrene is obtained by reducing the orthorhombic symmetry to hexagonal symmetry to reduce the computational cost. The band dispersion along with k-path used to compute the band dispersion for rubrene derivatives is presented in Fig. 7.4 and for BTBT derivatives in Fig. 7.5. As mentioned in the earlier sections, orthorhombic phase of rubrene exhibits better charge transport characteristics in comparison to triclinic, which in-turn is better than monoclinic phase. Accordingly, a decrease in band dispersion can be observed for triclinic in comparison to orthorhombic along the π stack direction (Fig. 7.4 : Γ → B for orthorhombic and Γ → A for triclinic) and also along the herringbone like packing direction (Fig. 7.4 : Γ → AB for triclinic). Comparatively, band dispersion in monoclinic phase is greatly suppressed, with respect to the first occupied band as well as the overall band dispersion. As mentioned in previous chapters of this dissertation, charge transport in crystalline organic semiconductors can be described to be band-like at low temperature, and becomes dominated by activated hopping at higher temperatures [3,29,30]. In crystalline organic semiconductors, a charge carrier can be considered to be localized leading to the polarization of electronic and nuclear subsystems of the surrounding lattice, leading to the formation of the quasi-particle, "polaron" [31]. Based on the interaction time (τ ) of the charge carrier with the surrounding medium, which inturn is related to the relaxation time (τ R ) necessary for the formation of polaron, τ can be classified as (i) fast and inertialess electronic residency time (τ e , of the order of 10 -16 s to 10 -15 s) or (ii) slow and inertial residency time (τ h , of the order of 10 -14 s to 10 -13 s) [31]. While τ e is related to the electronic polarization cloud formed by induced electric dipoles around the charge carrier, τ h is related to the localization of the charge carrier on the spatial lattice. This localization of the charge carrier on a spatial lattice, also leads to the interaction of the charge carrier with the intra-molecular vibrational modes, leading to the formation of lattice polarons with a characteristic time τ L . Furthermore, for a hopping-like transport, the carrier residence time also determines the mean time interval between hopping events. This entails that, the nature of hopping transport can be characterized by carrier residence time, τ , when τ = τ h [31] given by Eq. 7.1 [32].

τ (s -1 ) = h W = 2 3 • 10 -15 W (eV ) (7.1)
where W is the full effective band width. From Eq. 7.1 it follows that, as a general rule of thumb, an effective band width of the order of 0.1-0.2 eV is required to ensure that carrier residence time is smaller than typical electron-phonon coupling such that the molecule does not allow enough time to geometrically relax and trap the charge [32]. Further, in the Drude approximation, band mobilities can be obtained from Eq. 7.2, where e is the elementary charge, τ R is the relaxation time and m * is the effective mass obtained from the band dispersion energy at maximum energy points (Γ point) [33,34]

µ = eτ R m * (7.2)
In line with the above, band widths (W ) in the π stack direction along with the effective mass and transfer integrals for rubrene-polymorphs and BTBT derivatives are presented in Table 7.4. The band widths of rubrene orthorhombic and triclinic and all BTBT derivatives considered in this study is greater than 0.2 eV along the π stack direction. This entails that, band transport prevails in these materials, at appropriate temperature.

Table 7.4: Band width (W , eV), effective mass (m * ) and transfer integrals (J) (meV) of the rubrene polymorphs (O:orthorhombic, T:triclinic and M:monoclinic) and BTBT derivatives along the π stack direction. 7.4, it can be observed that the transfer integrals (effective mass) along the π-stacked direction decreases (increases) for rubrene polymorphs moving from orthorhombic to triclinic and to monoclinic phases, whereas for BTBT derivatives transfer integrals (effective masses) increase (decrease) along the π-stacked direction from BTBT to BTBT derivatives as a function of n in C n -BT BT . This entails that orthorhombic rubrene shows better charge transport when compared with other polymorphs as discussed in the previous sections, while for BTBT derivatives the charge transport ability increases as a function of n in C n -BT BT , which can be attributed to both the increase in transfer integrals and effective molecular packing of C n -BT BT as a function of n.

Strain-mobility trends

Strain-mobility trends of triclinic and monoclinic phases of rubrene, BTBT and C 8 -BTBT are reported in the subsequent sections as a function of transfer integrals and effective mass, for strains applied along the three crystallographic directions. Relative variation of mobility as a function of transfer integrals is obtained using [(J 2 i( ) /J 2 i ) -1] × 100 (reported as percentage variation), where J i( ) is the transfer integrals obtained when mechanical strain is applied along crystallographic axis i (a, b, c) and J i is the corresponding transfer integral when no strain is applied. Similarly, relative variation of mobility as a function of effective mass is obtained using [(m * i /m * i( ) )-1]×100. Relative variation with respect to transfer integrals and effective mass, will be, henceforth, referred to as ∆µ J and ∆µ m , respectively. However, as reported in Chapter 2, m * ∝ 1/J, and that the mobility as a function of transfer integrals depends on temperature and the choice of the semi-classical approach and the corresponding model [35], the relative mobility trends as a function of transfer integrals (∆µ J ) are also reported as [(J i( ) /J i ) -1] × 100, providing a qualitative comparison with that ∆µ m , which are reported in Appendix A4, for completeness.

In the previous chapter it was reported that, for orthorhombic rubrene, the relative variation of mobility is the highest when absolute values of transfer integrals are small. An extremely large relative variation does not indicate a relative mobility trend if the absolute values themselves are negligible. In order to avoid such spurious and often misleading observations, strain mobility trends are reported only along those directions where the absolute magnitude of transfer integrals is large enough, ie., in the present case, only values of J > 1 meV are considered to report observed strain mobility trends. Following this choice, for BTBT derivatives (BTBT and C 8 -BTBT) strain mobility trends are reported for strains applied along a and b axes, since the relative distance between molecules along c-axis is comparatively large. For instance, for both BTBT and C 8 -BTBT, the absolute values of transfer integrals along this direction is very small (≈ 10 -3 meV , for zero applied strain condition). Subsequently, strain-mobility trends as a function of effective mass are compared with those obtained as a function of transfer integrals. Again, for the sake of completeness, the transfer integral values as a function of applied strain along c-axis and the relative variations along c-axis for strains applied along a and b-axis, for BTBT, are reported in Appendix A4.

Rubrene: triclinic

Transfer integrals and effective mass of triclinic phase of rubrene at zero strain are reported in Fig. 7.6. Transfer integrals and effective mass along the π stack direction are lower than that of orthorhombic phase. However, transfer integrals along herringbone like packing -[110] direction are similar to that observed in orthorhombic phase. Relative variation of mobility for strains applied along a, b, c, measured along the aaxis (π stack direction) is reported in Fig. 7.7 and along [110] direction in Fig. 7.8. For strains applied along all crystallographic directions a relative increase (decrease) in transfer integrals is observed along π stack direction for compressive (tensile) strains, however differing in the magnitude. As expected, the opposite trend is observed for relative effective mass, since m * ∝ 1/J. Strain applied along the π stack direction (a-axis), induces a maximum variation of ≈ 13 % when transfer integrals (∆µ J ) are considered and ≈ 7 % when effective mass (∆µ m ) is considered. Relative variation is lower when strain is applied along c-axis (9 % and 4 % for ∆µ J and ∆µ m respectively) and further lowered for strains applied along b-axis (3.5 % and 2 % for ∆µ J and ∆µ m respectively). For strains applied along a and b crystallographic directions an increase (decrease) in relative mobility is observed along the [110] (herringbone) direction for compressive (tensile) strains, whereas a decrease (increase) of relative mobility is observed for compressive (tensile) strain along c-axis (see Fig. 7.8). For compressive strains along a and b axes the relative increase of ∆µ J is found to be 30 % and 17 % respectively, whereas that of ∆µ m is found to be 17 % and 14 % respectively. For compressive strain along c-axis a relative decrease of ∆µ J is found to be 7 % and that of ∆µ m is found to be 4 %. Further, for strains applied along a and b-axis, a deviation from linearity could be noticed for relative transfer integrals along [110] direction, indicating the sensitivity of the transfer integral variations as a function of distance [36].

Both transfer integrals and effective mass show similar trends for strain-mobility, with ∆µ J being larger than ∆µ m . The large value of ∆µ J with respect to ∆µ m , is not surprising, since relative variation with respect to transfer integrals is considered as a squared value of the relative transfer integrals, whereas the corresponding effective mass values are not squared. When the relative transfer integral are not squared, such that ∆µ J = J i( ) /J i ( , is the strain applied), the relative variations of both ∆µ J and ∆µ m , follow (almost) similar values, as reported in Appendix A4. In this chapter, the squared value of transfer integral is considered for ∆µ J and effective masses are not squared for ∆µ m , under the assumption that, ∆µ J reflect the hopping-like charge transport, whereas ∆µ m reflects the band-like transport, within the Drude approximation (Eq. 7.2). However, a qualitative similarity between ∆µ J and ∆µ m indicates the validity of both the approaches in capturing the strain-mobility trends.

Rubrene: Monoclinic

Transfer integrals and effective mass of monoclinic phase of rubrene at zero strain is reported in Fig. 7.9. Transfer integrals and effective masses along the π stack direction (along a axis) are lower than that of orthorhombic or triclinic phase. No herringbone-like arrangement is observed for monoclinic phases. Also, as the transfer integrals along [010] and [110] directions are lower than other polymorphs (see Fig. 7.9), strain mobility trends are reported only along π stack direction for strains applied along a, b and c axes.

For strains applied along a and c crystallographic directions an increase (decrease) in relative mobility is observed along a-axis (π stack : [100] direction), for compressive (tensile) strains, whereas almost no variation is observed for strains applied along b-axis. For compressive strains along a the relative increase of ∆µ J and ∆µ m along the a-axis is found to be 13 % and 7 % respectively, whereas for strain applied along c-axis the relative increase of ∆µ J and ∆µ m is found to be 6 % and 4 % respectively (see Fig.

BTBT

The principal packing plane for BTBT derivatives, corresponding to maximum inplane mobility, is the ab plane [19,20]. The contribution to mobility from molecules along the c-axis is negligible, as the transfer integrals along this direction is ≈ 10 -3 meV. Transfer integrals and effective masses along the π stack direction (b-axis) and along the herringbone packing ([110] direction) for zero strain condition (optimized crystal structure) are presented in Fig. 7.11. Strain mobility trends in BTBT are reported for strains applied along a and b (π stack direction). As mentioned earlier, the contribution to mobility from molecules in c direction is negligible (see Appendix 4) and henceforth, the impact of mechanical strain on relative variation of mobility is also neglected.

The variation of ∆µ J and ∆µ m along b-axis is presented in Fig. 7.12. For compressive (tensile) strains applied along π stack direction (b-axis) an increase (decrease) in relative mobility is observed along π stack direction whereas no variation is observed for strains applied along a-axis. For compressive strain applied along b-axis, relative increase of ∆µ J and ∆µ m along b-axis is 6 % and 2.5 % respectively.

The variation ∆µ J and ∆µ m along [110] direction is presented in Fig. 7.13. Similar to variation of mobility along π stack direction for strain applied along b-axis, for compressive (tensile) strains applied along b-axis, an increase (decrease) in relative mobility is observed along the [110] herringbone direction. For strain applied along a-axis, no variation in relative mobility is observed. For compressive strain applied along b-axis, the relative increase of ∆µ J and ∆µ m along b-axis is of 20 % and 15 % respectively. For compressive (tensile) strains applied along π stack direction (a-axis) and along the b-axis an increase (decrease) in relative mobility is observed along π stack direction. The variation of mobility as a function of ∆µ J and ∆µ m along a-axis is presented in Fig. 7.15. For compressive strain applied along a-axis, relative increase of ∆µ J and ∆µ m along a-axis is 6 % and 2.5 % respectively and that along b-axis is 2.5 % and 2 % respectively.

For compressive (tensile) strains applied along a-axis, an increase (decrease) in relative mobility is observed along the [110] herringbone direction, where as a decrease in relative mobility is observed when strain is applied along b-axis. The variation of mobility as a function of ∆µ J and ∆µ m along [110] direction is presented in Fig. 7. 16.

For compressive strain applied along a-axis, the relative increase of ∆µ J and ∆µ m along the [110] direction are of 7 % and 2 % respectively, where as for compressive strains along b-axis, decrease in ∆µ J and ∆µ m along [110] direction are of 2 % and 1 % respectively. Further, comparing BTBT and C 8 -BTBT, C 8 -BTBT gives higher (lower) transfer integrals (effective masses) along the π stack and also along the [110] directions (see Figs. 7.11 and 7.14). However, the relative variation of mobility along the π stack direction, with respect to both transfer integrals and effective masses, is similar for C 8 -BTBT and BTBT for strain applied along the π stack direction, whereas the relative mobility variation along the [110] direction is lower for C 8 -BTBT, for strains applied along a and b-axis.

Summary

Magnitude of relative mobility changes for triclinic and monoclinic rubrene, as well as BTBT and C 8 -BTBT derivatives obtained as a function of transfer integrals and effective masses, along the π stack direction for compressive strain in that direction is summarized in Table 7.5. Also summarized are the relative mobility variations in the herringbone-like packing, for applied strain along the crystallographic axis which generates maximum increase of relative mobility. It can be observed from Table 7.5 that, both increase in relative mobility along the π stacked direction as well as the the herringbone-like packing direction are due to compressive strains applied in the crystallographic direction that correspond to the π stacking, suggesting a strong electro-mechanical coupling between these directions.

Table 7.5: Summary of relative mobility variation along the π stack direction for triclinic and monoclinic rubrene, as well as BTBT and C 8 -BTBT derivatives for compressive strain applied along π stack direction (d π ). Relative variation of mobility along the [110] direction for strain applied along the crystallographic axis that generates maximum variation (d [110] ) is also reported. Transfer integrals (J) are in units of meV and relative mobility: ∆µ J and ∆µ m are reported as maximum percentage increase for compressive strains (0.8 %) with respect to the corresponding values (J and m * ) when no strain is applied.

Valence Band Shift (VBS)

In the previous sections, anisotropic strain-mobility trends are reported as a function of mechanical strain. Furthermore, Frisbie et al [37] reported that, mechanical strain also induces a variation in work function at device scale. Frisbie et al [37] showed that an increase (decrease) in work function is observed in rubrene (orthorhombic) crystals under compressive (tensile) strain. A qualitative agreement between theoretically calculated work function shift (∆W F ) as a function of shift in valence band maximum (∆V BM ) is also reported [37], entailing that ∆W F ∝ ∆V BM . In line with the above, valence band shift (VBS), ie., the shift in valence band maximum with respect to Fermi level, is computed for rubrene: triclinic and monoclinic phases, as well as BTBT and C 8 -BTBT crystals. This is done by calculating the energy difference between Fermi level and energy of the valence band at Γ point, given by E F -E Γ , for all the crystalline structures for no strain and strain conditions. The relative shift in valence band maximum (∆V BS) with respect to no strain condition is reported in Fig. 7.17. For every strain condition (no strain, compressive and tensile strains) a Gaussian smearing of 27.2 meV (0.002 Ry) is employed to extract the Fermi level. ∆V BS for the crystalline semi-conductor considered in this study, also exhibits an anisotropic response similar to that noticed with respect to relative mobility trends as a function of transfer integrals and effective mass. Subsequently, it follows that workfunction shows a similar trend, as ∆W F ∝ ∆V BS. Furthermore, the ∆V BS values for the compounds considered in this chapter, are very similar to the calculated values reported by Frisbie et al [37], however, these values can be much lower than the experimentally measured values [37] due to several factors like, presence of impurities, imperfections, relative workfunction shift of the material with respect to electrodes etc.

Conclusion

Strain-mobility trends for rubrene polymorphs (triclinic & monoclinic), BTBT and C 8 -BTBT are obtained as a function of transfer integrals and effective mass employing PW-DFT calculations as an alternative computational protocol to MD simulations to induce mechanical strain. A good qualitative agreement observed between relative mobility obtained as a function of transfer integrals (∆µ J ) and as a function of effective masses (∆µ m ) indicates the validity of both the approaches in capturing the strain-mobility trends. Relative increase in mobility along the π stack and the [110] directions for strain applied long thr π stack direction are the highest for triclinic rubrene with respect to the others. BTBT and C 8 -BTBT show similar mobility variations (∆µ J and ∆µ m ) along the π stack direction for strain applied along the same direction, whereas the relative mobility variation along the [110] direction is lower for C 8 -BTBT, for strains applied along both a and b-axis.

Relative strain-mobility trends reflect the inherent anisotropy present in the crystalline structures. Also, valence band shift shows a similar anisotropy that in-turn translates to the anisotropy in workfunction shift for strains applied along different directions. Both relative mobility and workfunction shifts are to be considered to derive conclusions with respect to experimental observations, as a relative increase in transfer integrals for compressive strains associated with a relative decrease in workfunction may not translate into the relative increase of mobility.

CHAPTER 8

CONFORMATION AND PROTONATION EFFECTS ON SULFONAMIDE FLUORINATED POLYMERS AS REPLACEMENT OF PSS IN PEDOT:PSS TRANSPARENT ELECTRODES

Introduction

Transparent and flexible conductive films are a crucial component of various organic electronic devices with promising applications such as transparent and flexible electrodes in Organic Light Emitting Diodes (OLEDs) [1,2], Organic Photovoltaics (OPVs) [3,4], and as active material in electrochromic devices [5], bio-compatible materials [6], electrochemical transistors [7] and as anti-static coatings. Due to its printability, non-toxicity and its reasonably high opto-electronic performance, one of the most used material for these applications is the conducting polymer complex based on Poly(3,4-ethylenedioxythiophene) (PEDOT) and poly(styrenesulfonate) (PSS) [1,2,3,4,5]. Despite several efforts made to improve the performance and understanding of PEDOT:PSS, a real breakthrough has not been achieved yet. PEDOT:PSS conducting complexes show hygroscopicity, acidity and poor wetting behavior on organic layers, limiting the scope of application. As an alternative to PSS, is the poly[4styrene sulfonyl trifluoromethyl sulfonyl] imide (PSTFSI) polyelectrolyte complex, showing, in blends with PEDOT, excellent printablity paving way towards printable (inkjet printed) conducting thin films as reported by Hoffmann et al [8].

In complexes with PEDOT, the anionic polyelectrolyte acts as counter ion and dopant. PSTFSI is known from its application in lithium batteries and has been shown to complex and dope PEDOT efficiently in aqueous dispersion [8]. Even if PSTFSI has the same styrenic backbone as PSS, it behaves differently in aqueous medium, due to the different nature of the acid functionality (sulfonate vs bissulfonate).

Polymeric PSTFSI is synthesized prior to the addition of the polyelectrolyte to EDOT, which polymerizes in its presence. During the polymerization reaction, the protons cleaved from the EDOT monomer during the polymerization lead to a drastic decrease of the pH in the synthesis medium. PEDOT is then further, partially, oxidized into the positively charged, so called "doped" conducting form of PEDOT. Therefore, the characteristics of the polyelectrolyte are expected to be decisive for the chemical-physical properties of the resulting PEDOT:polyelectrolyte complex. Thus, it is of great interest to study the charge, acid-base behavior or conformation of the anionic polyelectrolyte, and how those characteristics are related to the interaction with PEDOT.

Experimental studies were carried out to obtain a fundamental understanding of the complexation and doping mechanism of PEDOT:PSTFSI. In order to elucidate the relation between the inherent characteristics of PSTFSI polyelectrolyte and the resulting properties of PEDOT:PSTFSI complex such as absorption and conductivity, series of experimental investigations were done by Hoffmann et al., such as XPS, UV/Vis, electrophoresis, electron microscopy and rheological analysis. Experimental observations indicate the presence of both protonated and de-protonated states of nitrogen atoms in PSTFSI. Subsequently, it was found that, when the conducting PEDOT:PSTFSI complex is synthesized using PSTFSI with completely protonated nitrogen atoms, conductivity of the resulting PEDOT:PSTFSI complex is increased, with respect to the complex when PSTFSI with completely de-protonated nitrogen atoms is used. Important experimental observations and the subsequent findings, in brief, are as follows.

• The absorption coefficient of PEDOT increases linearly with the increase of molar concentration of PEDOT, whereas the conductivity increases up to a molar EDOT/STFSI ratio of ≈ 0.9 and then decreases.

• The amount of free PSTFSI in PEDOT:PSTFSI dispersions was reduced with increasing PEDOT concentration. At a molar EDOT/STFSI ratio of 1.10, almost no free PSTFSI was detected by capillary electrophoresis

• XPS measurements revealed that PSTFSI stabilizes in a mixed charged form, bearing 60 to 70% protonated nitrogen states in the PSTFSI repeating units and 30 to 40% negatively charged nitrogen states (peak at 398 eV)

Computational approach

The protonation state of PSTFSI can be considered as one of the decisive parameters controlling the morphology and physical properties of PEDOT:PSTFSI complexes.

In line with the above experimental findings, a series of theoretical and computational protocols were employed to understand the inherent characteristics of PSTFSI and to provide insights to experimentalists. First, DFT calculations were performed to estimate the acidity constant (pK a ) of STFSI monomer and oligomers (dimers and trimers) using implicit solvation models. Subsequent to the observation that hydrogen bonds are present in STFSI oligomers, Nitrogen 1s binding energies were extracted employing natural-bonding-orbital (NBO) analysis to measure the shift of 1s binding energy as a function of hydrogen bonding with respect to experimental XPS analysis. Further, experimental IR and Raman spectra of PSTFSI were compared with the corresponding theoretical spectra of monomer, dimer and trimer STFSI units as a function of protonation state, to estimate the percentage in experimental samples and to complement the XPS protonation. Molecular Dynamics (MD) simulations were then performed on oligomers of PSTFSI, treating water explicitly, to study the nature of hydrogen bondings. Subsequently, free energy of conformation of STFSI oligomers were analyzed to quantify the impact of protonation on polymer/oligomer conformation.

Estimation of acidity constants (pK a ) of STFSI oligomers

The amide group in STFSI can exist in three forms depending on the pH, namely N H, N -, N H + 2 . Neglecting the later which is typically a superacid (cal) the protonation reaction of R-N H, R-N -is studied, wherein, R-N H aq → R-N (-) aq + H (+) aq

The equilibrium constant for this reaction is the acidity constant,

K a = [R-N -][H + aq ] R-N H , with, pK a = -log 10 K a (8.1)
From the thermodynamic definition of equilibrium constant ∆G * = -RT ln K a , it follows that,

pK a = ∆G deprot(aq) 2.303RT (8.2)
where,∆G deprot(aq

) = G aq,R-N (-) -G aq,R-N H + G aq,H (+)
In order to establish a reliable computational methodology, DFT calculations were first performed on molecules with experimentally known pK a values at the PBE/def2-TZVP level of theory with continuous-polarized continuum (CPCM) [9,10], an implicit solvation model. Water is used as implicit solvent to compute Gibbs free energies (G * ). Frequency calculations are performed to account for the zero point vibrational energy (ZVPE), thermal and entropies corrections at the standard temperature (T) and pressure (P) conditions, with T = 298.15 K and P = 1 bar. Theoretical pK a is computed, employing the "direct method" [11,12] using Eq (8.

2) [12], in which G aq,R-N (-) and G aq,R-N H are calculated Gibbs free energies with zero point energy correction and G aq,H (+) is the hydration enthalpy of hydrogen in water (solvent) obtained from experimental data (-265.9 kcal/mol [13]). The calculations were performed for a series of compounds with known pK a and chemically similar to STFSI, as shown in Table : 8.1.

Computed pK a values differ from experimental ones [11,12,13], due to inherent approximations introduced in the DFT functionals, as well as the continuum solvation models, which does not account for the specific interactions such as H-bonding or London dispersions. This may lead to systematic errors in the estimation of free energy of solvation, that may be corrected by simple rescaling schemes. Linear regression between calculated and experimental pK a was employed to correct the error through two adjustable parameters, using Eq. 8.3.

pK a (Corrected) = pK a (T heory) -C S (8.3)
CF 3 -SO 2 -N H -SO 2 -C 6 H 5 -CF 3 4.54 3.75 [20] CF 3 -CF 2 -SO 2 -N H -SO 2 -CF 2 -CF 3 -3.65 -3.0 [19]
The acidity constant of PSTFSI monomer, dimer and trimer were then calculated, applying the appropriate correction obtained from the benchmark calculations. The pK a values reported in Table 8.2 suggest that the monomer is not a strong acid. Further, from Table 8.2 it appears that, with the increase in oligomer length and of deprotonation, further deprotonation of the system becomes increasingly difficult (increase in pK a with increase in oligomer length). It subsequently follows that, complete removal of hydrogens from the system is difficult and in a polymer STFSI would retain a fraction of protonated nitrogens, in line with polyelectrolyte theory [21,22,23]. However, the values of pK a 1 , pK a 2 and pK a 3 reported in The presence of intra-chain hydrogen bonding between NH group of a given STFSI neutral monomer and one of the S = O group of neighboring STFSI unit is different from what noticed by Vaden et al [20], wherein a N -bridged hydrogen bond (N -• • • H + • • • N -) was observed, but is reminiscent of the H-bond network present in TFSI crystal structure [19].

Since experimental XPS spectra indicate the presence of different protonation states in the system, experimental IR and Raman spectra were compared with theoretical calculations, to address the coexistence of protonated and de-protonated nitrogens. IR and Raman spectra are computed for monomers and trimers of STFSI at the PBE/def2-TZVP level of theory with CPCM implicit solvation model. This is done by identifying the vibration modes that arise from the protonation state of nitrogen in the monomer/oligomer. Theoretical IR & Raman spectra of the monomer are plotted in Fig. 8.3 and Fig. 8.4, respectively. IR and Raman (N H and N -) spectra are linearly combined to simulate vibration spectra as a function of percentages of protonated and de-protonated states. When these spectra are compared with experiments, a rough estimate of percentage of protonated nitrogen states in PSTFSI can be obtained. In IR spectra, presence of peaks around ≈ 750 cm -1 corresponds to the protonated nitrogen (symmetric stretching modes of S -N ), whereas the corresponding band red shifts to ≈ 1000 cm -1 when deprotonated nitrogen state is present. This could be attributed to the S -N bond length variation as a function of protonation state of nitrogen.

Similarly, in Raman spectra, presence of peaks around ≈ 760 cm -1 corresponds to the protonated nitrogen (symmetric stretching modes of S -N ), whereas this peaks blue shifts to ≈ 735 cm -1 when deprotonated nitrogen state is present. is obtained when STFSI oligomer is de-protonated at the terminal nitrogen atom (single de-protonation). In summary, calculations and experiments consistently indicate that only about one third of the N H groups in the polyelectrolyte is dissociated at acidic pH, a higher fraction with respect to the standard poly(styrenesulfonate):PSS [24,25].

The presence of protonated nitrogens in PSTFSI, is established and the possibility of hydrogen bonding as observed from DFT calculations (refer to Fig. 8.2) is considered.

To study the impact of hydrogen bonding on energy density of nitrogen atom, the energy of the 1s orbital of the nitrogen atoms is extracted using natural bonding orbital analysis (NBO) from DFT optimized geometries (at the PBE/dft2-TZVP, implicit solvent water/CPCM) of trimer units of STFSI, as a function of the protonation state of nitrogen. The shift in calculated (experimental) 1s orbital energy of nitrogen between protonated (N H) and de-protonated (N -) nitrogen is 1.4 eV (1.5 eV) whereas between protonated (N H) and doubly protonated (N H + 2 ) nitrogen is 1.5 eV (1.8 eV). Shift in 1s orbital energy of nitrogen that participates in hydrogen bonding (N H • • • O) with respect to that which does not participate in hydrogen bonding (N H) is 0.9 eV. Calculations indicates that presence of hydrogen bonding shifts the 1s energy of nitrogen towards higher binding energy by 0.9 eV , entailing that the additional 1s nitrogen energy peaks from experimental XPS spectra (with a shift of 0.8 eV) can be attributed to presence of hydrogen bonding in experimental samples. Although a direct comparison between experimental shift in binding energy and shift in 1s energy of nitrogen atom from DFT/NBO analysis is not fully quantitative, shift in binding energy can be attributed to the effect of neighboring environment on corelevel atoms of nitrogen atom. The presence of hydrogen bonding, unambiguously leads to a decrease in electron density of nitrogen atoms participating in hydrogen bonding. This reduction of electron density leads to the variation of energy of corelevel atoms due to reduced screening of nuclear charge [26], entailing that a linear relation between theoretical shift in 1s energy of nitrogen atoms and experimental shift in binding energy of 1s nitrogen orbitals can indeed be established as an effect of neighboring environment and attributed to hydrogen bonding [26].

Molecular Dynamics Simulations

Molecular dynamics simulations are performed with NAMD [27] under constant pressure (1 bar) and temperature (298 K) on STFSI oligomers treating, water explicitly to better understand the impact of hydrogen bonding on the oligomer conformation as a function of protonation state. Three different sets of simulations are performed with each system consisting of one STFSI decamer and 5000 water molecules. STFSI was described with a customized general AMBER force field [28] where for both protonated (N H) and deprotonated (N -) forms, the atomic charges were calculated at PBE0/6-311G * level. The potentials for the four torsions (

• • •) along the phenyl • • • SO 2 • • • N -(H + ) • • • SO 2 • • • CF 3 side
group were calculated at the same level of theory and inserted in the force field following the scheme detailed in reference [29] and for water the well-known and reliable "modified TIP3P" force field was used [30]. The three STFSI decamer systems differ in the nature of oligomer with (i) all de-protonated nitrogen atoms, henceforth referred to as (N ) 10 , (ii) de-protonated and protonated nitrogen state of alternating STFSI units, henceforth refereed to as (N -N H) 5 and (iii) all protonated nitrogen states, henceforth refereed to as (N H) π-stacked-like regions are present in STFSI oligomers for MD simulations treating water molecules explicitly, similar to that noticed from DFT calculations of STFSI dimer and trimer with implicit solvation scheme (refer to Fig. 8.2). π-stacked conformations are stabilized by intra-chain hydrogen bondings between protonated nitrogen of STFSI unit and oxygen of neighboring STFSI unit, and protonated nitrogen of STFSI unit and de-protonated nitrogen of neighboring STFSI unit, although the latter is less frequently noticed. No π-stacking is noticed when only de-protonated nitrogen states are present (N 10 ) as the phenyl units are drifted apart due to the higher repulsive forces between the negatively charged nitrogens.

To provide a quantitative assessment of the oligomer conformation, parameters like oligomer length, π-stack distance between center of mass of phenyl units and average number of hydrogen bonds per oligomer per MD simulation time step (computed with a cutoff distance of 3.0 Å between donor and acceptor moieties and a cut-off angle of 25 PSTFSI oligomers with protonated nitrogen states show longer end-to-end distance with stacked regions along the chain, while complete removal of amide protons results in the collapse of oligomer backbone and no stacked regions are noticed. This outcome can be easily attributed to the strong repulsive forces between the electron rich STFSI units with de-protonated nitrogen atoms resulting in phenyl units drifting away from each other to maximize the distance between the negative charges. When protonated nitrogen is present, π stacked regions are noticed between the phenyl units, strengthened by the intra-chain hydrogen bonding as reported in Hydrogen bonds (Dynamic)

( Å) r E-E ( Å) n N H•••O n N H•••N n N H•••ODW n ODW •••N r initial r f inal (N ) 10 7
NH • • • On (a) n N H•••ODW (b) n N H•••O MD simulation time (25ns)
Hydrogen bonds (Dynamic)

HO -H • • • N - MD simulation time (25ns) 0 2 4 6 8 10
Hydrogen bonds (Dynamic) Hydrogen bonds (Dynamic)

NH • • • N - n (c) n ODW D •••N - (d) n N H•••N -
NH • • • O n (a) n N H•••ODW (b) n N H•••O Figure

Conformational free energy difference in solution

In order to estimate the relative stability of the conformers (where π-stacking is noticed to depend on protonation state), free energy of interaction is computed between STFSI units using the Adaptive Biasing Force (ABF) method [31] and MD simulations treating water explicitly. Simulations were performed using NAMD [27] under constant pressure (1 bar) and temperature (298 K) for 100 ns simulation time with 1500 water molecules used per timer. Force field parametrization is similar to that reported in Sec. 8.4. Free energy is computed for four different sets of identical STFSI trimers possessing, (i) all protonated amide nitrogen states : (N H) 3 , (ii) all de-protonated nitrogens : (N ) 3 , (iii) two protonated nitrogen on the terminal STFSI units and one de-protonated nitrogen in the central unit : N H-N --N H and (iv) two de-protonated nitrogen on the terminal STFSI units and one protonated nitrogen in the central unit : N --N H-N -. Free energy, reported in Fig. 8.11, is computed as a function of distance between interacting nitrogen atom of the central STFSI units, and the sulfur atom of one neighboring STFSI units, with both the atoms defining the collective variables to obtain the free energy of interaction. Stable states corresponding to free energy minimum for all the four configurations from Fig. 8.11 is presented in Fig. 8.12.

For neutral trimer, the minimum of free energy is found at ≈ 4 Å. This is also the situation when central de-protonated nitrogen present (N H-N --N H). For trimer with all de-protonated nitrogens minimum of free energy is found at a distance of ≈ 12 Å, which is also the situation when central protonated nitrogen (N --N H-N -), however with an additional minimum is present at ≈ 6 Å. From the free energy properties it follows that, the presence of protonated amide groups stabilizes conformations in which two consecutive phenyl-TFSI side groups are close in space and forms N • • • O = S hydrogen bonds between them. Consequently it follows that, the presence of protonated states in PSTFSI allows for the formation of π stacked regions between phenyl units, wherein the oligomers exhibit a linear π-stacked like conformation.

Conclusion

The acidity constant (pK a ) of STFSI monomer is ≈ 0.34 pK a units, but STFSI is not a strong acid and de-protonated form can attract hydrogen atoms. Increase of pK a as a function of oligomer size, in line with polyelectrolyte theory, indicates that PSTFSI can retain protonated nitrogen states and complete de-protonation of nitrogen is unlikely.

Superimposing IR and Raman spectra of STFSI monomer as a function of protonated and de-protonated nitrogen atoms, as well as the corresponding spectra of trimers indicate that the percentage of protonated nitrogen atoms in PSTFSI is ≈ 60 %, which is in line with the observation from experimental XPS spectra. Subsequently, this entails that, ≈ 60 % of protonated nitrogen states can be considered as the saturation limit of PSTFSI polymer with the remaining ≈ 40 % of de-protonated nitrogen states available for doping PEDOT units.

From MD simulations of PSTFSI oligomers as a function of protonation state of nitrogen atoms, wherein water is treated explicitly, it emerges that,

• The presence of protonated nitrogen atoms leads to the formation of π-stacked arrangements of the phenyl units of neighboring STFSI units, therein resulting in the oligomer conformation to be linear.

• Oligomers with completely de-protonated nitrogen atoms do not show any πstacked like arrangement, with the neighboring STFSI units that are drifted apart due to the higher repulsive forces between the negatively charged amide nitrogens.

• The conformational free energy analysis also confirms that, the structure with linearly stacked (π stacked-like) arrangement is the lowest energy conformation, when all the nitrogen atoms are protonated or when the percentage of protonated nitrogen atoms is higher than that of de-protonated ones.

It can be then hypothesized that, starting PEDOT polymerization from all-protonated PSTFSI (starting from acid form: PSTFSI-H), results in a conformation with regions of linearly stacked phenyl functional groups, stabilized by higher number of intra-chain hydrogen bonds, which in-turn influence the PEDOT conformation, this linearly arranged conformation could act as template in favor of formation of long, straight PEDOT chains. In such system the intermolecular hopping is modulated and CHAPTER 9

CONCLUSIONS & PERSPECTIVES

This dissertation addresses the relation between the physical and chemical properties of organic materials like chemical structure, molecular conformation, supra-molecular arrangement and the resulting charge-transport characteristics. The results of the dissertation can be split into three main areas, wherein (i) charge carrier mobilities of crystalline tetravalent Group-IV-phthalocyanines and amorphous TFB samples were estimated by a home built Kinetic Monte Carlo code, were presented in chapter 4 and 5 (ii) Electro-mechanical response of crystalline organic semiconductors with two different methodologies to induce mechanical strain, either by Molecular Dynamics simulations or by Plane-Wave DFT calculations, and the resulting strain-mobility trends of rubrene polymorphs and BTBT derivatives were presented in chapters 6 and 7 and (iii) the protonation state of PSTFSI, a poly-anion that acts as a dopant when combined with PEDOT to obtain PEDOT:PSTFSI conducting complexes, and their impact on the polymer conformation and also on the resulting conducting complexes was discussed in chapter 8.

Charge carrier mobilities-KMC simulations: Charge transport characteristics in crystalline tetravalent Group-IV-phthalocyanines was studied employing the Marcus-Levich-Jortner rate formalism, to establish the structure-charge carrier mobility trends of different compounds that vary in central metalloid core, number/position of fluorine atoms on the axial phenoxy group or as a function of methyl / iodine atoms on the phenoxy group. These changes in functional groups induce subtle variations in crystalline packing, that in-turn lead to large variations in charge transport characteristics. Ambipolar and dimensionality of charge transport of these materials was presented, with specific emphasis on the contribution of energetic disorder, that is inherent to organic semiconducting materials, to simulated mobilities. Subsequent to this study, charge transport in amorphous semiconducting fluorene-triphenylamine (TFB) copolymer was studied, employing in this case the simpler Marcus formalism, and explicitly treating the different energetic disorder conditions which can be static or dynamic in time. A systematic comparison to available experimental results was made to highlight the impact of energetic fluctuations on charge carrier mobilities. A semi-quantitative estimate of calculated mobilities is obtained with respect to experimental measurements, only when the energetic fluctuations are included in the simulations. It evolves that, the inherent energetic disorder that is present in the organic electronic materials is crucial in governing the charge carrier mobilities and by sensibly including these fluctuations into simulations models a quantitative agreement between simulations and experimental observations can be obtained. A methodology is proposed to include the thermal fluctuations towards providing a semi-quantitative estimate of charge-carrier mobilities, with respect to available experimental observations. This modeling methodology can be used for material screening of different materials for organic electronics, at reduced computational cost.

Electro-mechanical response in crystalline organic semiconductors: Crystalline rubrene and its polymorphs, as well as BTBT derivatives (well studied high mobility organic materials) were subjected to mechanical strain and their electronic response was analyzed. Two different methodologies were employed to induce mechanical strain, (i) Molecular Dynamics simulations and (ii) Plane-Wave DFT calculations. Compressive and tensile strains were applied to the super-cell of orthorhombic rubrene and relative change in mobility as a function of applied strain with respect to the mobility under zero strain was obtained. Calculated relative mobility variations were compared to experimental measurements that were performed by mounting the crystalline rubrene on a novel cantilever set-up, wherein mechanical strain was induced in rubrene single crystals by applying a force on the cantilever support.

A semi-quantitative agreement between theoretical calculations and experimental measurements is obtained. By a systematic comparison between calculations and measurements, it evolves that, (i) strain-mobility trends are anisotropic in nature, wherein strain applied along different crystallographic directions generates different strain-mobility trends and (ii) thermal fluctuations do not impact the strain mobility trends. Subsequent to this study, polymorphs of rubrene (triclinic and monoclinic forms) and derivatives of BTBT (BTBT and C 8 -BTBT) were subjected to mechanical strain by PW-DFT calculations and strain-mobility trends were obtained both as a function of transfer integrals as well as effective masses, thereby providing a perspective from both hopping transport and band-like transport.

Analyzing the strain-mobility trends of rubrene polymorphs and BTBT derivatives, it emerges that strain mobility trends in crystalline organic materials cannot be generalized as the change in crystalline structure (crystallographic packing) also changes the strain-mobility trends. This variation acts on top of the anisotropic strain mobility trends that are inherent for a material class, which are a function of strain applied along different crystallographic directions. For orthorhombic rubrene, where in MD simulations were employed to induce mechanical strain, thermal fluctuations hinder the identification of parameters that are responsible for inherent strain-mobility anisotropy of the material. PW-DFT calculations can be used to identify to param-eters that generate this anisotropy, as thermal fluctuations are absent. Although thermal fluctuation do not impact strain-mobility trends in orthorhombic rubrene, the impact of thermal fluctuations on strain-mobility trends of other class of materials cannot be generalized. The complex interplay of different forces in action and the interconnectivity of electronic and mechanical response of the material which in-turn originates from the crystalline packing and the elastic response of the material along different crystallographic directions, results in the generalization of strain-mobility trends to be very complicated.

Protonation-dependent conformation of polyelectrolyte and its role in governing the conductivity of polymeric conducting complexes: Protonation state of polymeric bis(sulfonyl)imide polystyrenes, that are currently employed as counter-ions and dopants for conducting poly (3,4-ethylenedioxythiophene), PEDOT, resulting in PEDOT-polyelectrolyte conducting complexes, was studied by employing MD simulations and DFT calculations. The inherent characteristics of the polyelectrolyte like its acid-base behavior, protonation state and conformation of the polymer driven by the protonation state, were analyzed in conjunction with available experimental data. This Molecular Dynamics simulation study, although focused only on the intra-molecular interactions that are inherent to the PSTFSI as a function of protonation state and the inter-molecular interactions with explicit water molecules, reveals the impact of subtle interactions within the poly-electrolyte materials that can induce large variations on the conductivity of resulting conductive complexes. The work can be extended by considering more PSTFSI oligomers to consider the intermolecular interactions between the independent PSTFSI units. A further extension of the work is to include the PEDOT oligomers to obtain a blend of PEDOT:PSTFSI conducting complex and study the electro-static interactions of PEDOT units as a function of protonation state of PSTFSI.

This dissertation addressing the relation between structure-charge transport properties of different class of organic electronic materials, with specific emphasis on the impact energetic interactions (that are strongly coupled and vary in terms of magnitude, strength and type of interactions), on the charge transport nature of organic electronic materials, demonstrate the subtlety of charge transport nature of these class of materials, emphasizing the need of computational studies to better understand these materials. By systematic comparisons to experimental observations, the role of these energetic interactions can be de-coupled, thereby providing additional insights to experimentalists. These insights from the application of Computational Chemistry to study organic electronic materials can propel the field of organic electronics towards providing stable and highly efficient organic electronic devices.

: GROUP IV PHTHALOCYANINES

Inter-Molecular distances

Inter-molecular distances between the molecules (in units of Å) along the crystallographic planes is presented in the following. [4,5]. Atomic charges were calculated on the optimized structure with the electrostatic potential (ESP) fitting method from B3LYP/cc-pVDZ calculations. The reproduction of the molecular dipole moment was imposed in the fit. The geometry optimization and the calculation of ESP charges are performed with the Gaussian09 software [6].

In our classical model of TFB we impose to all fluorene and triphenylamine comonomers to have the same atomic charges, so that the polymer/oligomer is a sequence of chemically bounded identical units. Atoms of co-monomers were labelled as shown in the upper panel of figure A2-2 and the charge on each atom type was averaged over the central co-monomers of the oligomer (4 fluorenes, 3 triphenylamines), excluding the terminal units. The bottom panels of Fig. A2-2 show ESP charges (red crosses) and the average charge (black squares) for each atom type of fluorene (left) and triphenylamine (right). The average value is representative for the charge of each atom type and it is therefore used as a good approximation for the simulation of oligomers or polymers. United atoms of fluorene octyl chains, corresponding to aliphatic CH, CH 2 , and CH 3 groups, are assumed to be electrically neutral. A:18

Analysis of inter-molecular degrees of freedom

In the attempt of rationalizing the strain-induced variation of the electronic couplings J between the nearest neighbours along the a cell vector, relevant intermolecular degrees of freedom were characterized as a function of strain for 1000 snapshot taken from MD simulations along a, b or c, in both uniaxial strain and uniaxial stress conditions. Those included the lateral shift of one molecule along one of its three axes (as shown in Fig. A3-2, A3-3, A3-4), and the rotation of one molecule with respect to its three symmetry axes (see Fig. A3-5). In addition, the transfer integral of a dimer taken from the crystallographic structure was calculated while varying the above-mentioned intermolecular coordinates. The comparison between the trends allowed thus to assess the relative importance of each coordinate with respect to the transfer integral variation upon strain.

As far as the rotation angles are concerned, it is shown from Fig. A3-5 : [b] and [c] that, despite a certain dependence of angles θ and χ from strain in the case of uniaxial strain can be found, the absolute variations is always of the order of 10 -2 degrees. The most sensitive angle with respect to the transfer integral variation is χ, but while J a increases with χ, χ increases with positive strain, which in turn was found to decrease J a . Instead, transfer integral shows a monotonic decrease with increasing ϕ and θ, although only θ increases slightly with positive strain in the case of uniaxial strain. This brings us to the conclusion that there is not a specific intermolecular mode that can be held solely responsible for the observed variation of mobility with strain.

A:22 A:27

Table A3-1: Calculated mean : < J >, < J t >,, ratio of < (J t ) 2 > / < J t > 2 , standard deviation (σ), η, skewness and kurtosis of the distribution of transfer integrals along a axis -(π) stacked direction, for uniaxial strain applied along a, b and c axis : aa , bb and cc , respectively. All values, except η are in meV. η is dimensionless quantity.

(%) < J > < J t > <(J t ) A3-3: Calculated mean : < J >, < J t >, ratio of < (J t ) 2 > / < J t > 2 , standard deviation (σ), η, skewness and kurtosis of the distribution of transfer integrals along c axis, for uniaxial strain applied along a, b and c axis : aa , bb and cc , respectively. All values, except η are in meV. η is dimensionless quantity.

(%) < J > < J t > <(J t ) A3-4: Calculated mean : < J >, < J t >, ratio of (< (J t ) 2 > / < J t > 2 , standard deviation (σ), η, skewness and kurtosis of the distribution of transfer integrals along a axis -(π) stacked direction, for uniaxial stress applied along a, b and c axis : aa , bb and cc , respectively. All values, except η are in meV. η is dimensionless quantity.

(%) < J > < J t > <(J t ) A3-6: Calculated mean : < J >, < J t >, ratio of < (J t ) 2 > / < J t > 2 , standard deviation (σ), η, skewness and kurtosis of the distribution of transfer integrals along c axis, for uniaxial stress applied along a, b and c axis : aa , bb and cc , respectively. All values, except η are in meV. η is dimensionless quantity.

(%) < J > < J t > <(J t ) A:39 A:42

Figure 1 . 1 :

 11 Figure 1.1: Schematic representation of OPVs and OLEDs depicting several layers of specific functionality. Left(a): Device architecture and working principle of organic photovoltaics device and Right(b) : Device architecture and working principle of organic light emitting diode (adapted from references [27, 28])

Figure 2 . 1 :

 21 Figure 2.1: Temperature dependence of mobility and classification of different transport regimes as a function of strength (strong: g 2 P >> 1 or weak: g 2 P << 1) of local electron-phonon coupling. Adapted from Ref [4, 5]. E C , E V and E T correspond to energy of conduction band minimum, valence band maximum and energy of trap states, respectively.

Figure 2 . 2 :

 22 Figure 2.2: Schematic representation of internal reorganization energy from the fourpoint adiabatic potential (AP) approach.

Figure 2 . 4 :

 24 Figure 2.4: Schematic representation of electronic splitting between HOMO/LUMO levels of two pentacene molecules superimposed in cofacial configuration. The splitting is two times the electron coupling (J HOM O/LU M O )

Figure 2 . 5 :

 25 Figure 2.5: Schematic representation of energy dispersion in 1D TBA approach.

Figure 2 . 7 :

 27 Figure 2.7: Schematic representation of energy terms involved in molecular mechanics. Adapted from Ref. [49].

Figure 3 . 1 :

 31 Figure 3.1: Discrete Gaussian distribution of site energies obtained with σ = 2.0 k B T = 51.4 meV for (a) Left: 100 lattice sites, (b) Middle: 2000 lattice sites and (c) Right: 8000 lattice sites.

Figure 3 . 2 :

 32 Figure 3.2: Charge carrier mobilities issued from KMC simulations obtained in this dissertation (empty dots), as a function of diagonal energetic disorder: σ (in units of k B T ) and comparison to reference results from [1] (dashed lines).

3 . 4 .Figure 3 . 4 :

 3434 Figure 3.4: Evolution of the site energy distribution as a function of the simulation end distance (d E ), for a lattice model of 8000 hopping sites with σ set to 2 k B T . E Eq. = -σ 2 /k B T correspond to the expectation value of mean energy of site energy when equilibrium is attained.

Figure 3 . 5 :

 35 Figure 3.5: KMC simulations, with d E = 2 µm and N = 250, on a lattice model with 8000 hopping sites as a function of σ . Lines correspond to µ A and dots correspond to µ L . Deviation between µ A and µ L can be noticed for σ = 2 k B T indicating that simulations have not attained complete equilibrium.

Figure 3 . 7 :

 37 Figure 3.7: Hole mobility of the pentacene crystal as a function of the electric field orientation in the ab plane. Radial grid corresponds to mobility values in units of cm 2 /V s. Left: this work and Right: adapted from reference [8].

Figure 3 . 8 :

 38 Figure 3.8: Hole mobility obtained for the pentacene crystal as a function of electric field orientation in the ab plane when no diagonal disorder is considered (black line) and when static ( S G , red) and dynamic ( G , blue) diagonal disorder, with σ = 2 k B T , are considered. Radial grid correspond to mobility values in units of cm 2 /V s.

Figure 3 . 10 :

 310 Figure 3.10: Hole mobility of pentacene crystal as a function of electric field orientation in ab plane when no off-diagonal disorder is considered (black-squares) and when diagonal and off-diagonal energetic disorders are introduced as (i) S G |J σ G : bluedashed curve, (ii) S G |J G : red-solid curve, (iii) G |J σ G : green-hexagons curve, (iv) G |J G : magenta-pentagons curve. Radial grid correspond to mobility values in units of cm 2 /V s.

Figure 4 . 1 :

 41 Figure 4.1: Chemical structure of metal (metalloid) -substituted phthalocyanines studied in this chapter.

Table 4 . 1 :

 41 Molecular reference and corresponding crystal structure parameters. Z represent number of formula units in the cell, and M,T,O represent Monoclinic, Triclinic and Orthorhombic symmetries respectively. Space Group (SG) is represented in Hermann-Mauguin (H-M) notation. Molecule Cell Parameters Ref. Unit a( Å) b( Å) c( Å) 10.31 11.52 72.5 70.1 83.8 1 T:P1 5 8.41 10.16 11.51 72.5 69.5 81.3 1 T:P1 6 8.44 13.60 15.23 90.0 97.4 90.0 2 M:P 21/c 7 10.20 17.69 10.79 90.0 117.6 90.0 2 M:P 21/n 8 9.63 19.77 10.04 90.0 118.6 90.0 2 M:P 21/c 9 10.25 16.56 11.51 90.0 115.8 90.0 2 M:P 21 10 12.64 19.58 7.54 90.0 103.2 90.0 2 M:P 21/c

4. 3 . 1

 31 Internal reorganization energies, Transfer integrals & Huang-Rhys parameters Experimental and theoretical (calculated at DFT/B3LYP/631-G(d) level) Electron Affinity (EA) and Ionization Potential (IP), along with calculated HOMO and LUMO energies of all the molecules considered in this chapter (see Fig. 4.1) are reported in Table

Figure 4 . 2 :

 42 Figure 4.2: Compound 1 molecule and crystallographic packing along the π stack direction (π stack direction is along a axis)

  .6 and the results obtained from the parameters are reported in Fig. 4.3. Table 4.6: Simulation conditions employed to study the effect of energetic disorder with respect to Fig. 4.3. 0 Identical site energies, ∆ ij = 0 S G Gaussian distribution of site energies with σ = 51.4 meV , static in time: disorder generated at the beginning of every KMC simulation run D G

Figure 4 . 3 :

 43 Figure 4.3: Compound 1 mobility: Hole mobility (at 300 K) along the π stacked direction (a axis, with electric field orientation along the a axis) computed using Marcus (empty symbols with dashed lines) and MLJ formalisms (filled symbols), as a function of variation of electric field for different energetic disorder configurations (as per Table 4.6). Mobilities obtained for reorganization energies, left(a): λ=0.1 eV and right(b): λ=0.2 eV are reported, for comparison.

5 x 2 JFigure 4 . 4 :

 5244 Figure 4.4: Left(a): Hole mobility of compound 1 along π stack direction, as a function of variation of percentage of disorder on J π and (b) Right:J 2 π + σ 2 as a function of increasing value of fluctuation on J π , σ J .

  .

Figure 4 . 5 :

 45 Figure 4.5: Contribution of energy from Gaussian distribution of percentage disorder on J π for hole mobility of compound 1. (a) Left: Variation in mobility as a function od reorganization energy (λ S ) and Right: Variation of mobility as a function of electric field (E in units of V /cm).

Figure 4 . 6 :

 46 Figure 4.6: Compound 1 mobility: Hole and Electron mobility (at 300 K) along the π stacked direction (a axis, with electric field orientation along the a axis) computed using MLJ formalisms, as a function of variation of electric field for different energetic disorder configurations (as per Table 4.6). Mobilities obtained for, left(a): hole transport and right(b): electron transport are reported, for comparison.

Figure 4 . 7 :

 47 Figure 4.7: Electron and hole transport of compound 1 (using MLJ formalism and λ S = 0.2 eV) as a function of different simulation conditions: Left(a): ab plane, Middle(b): ac plane and Right (c): bc plane. a b and c represent crystallographic axis direction. Top: images with non-zero transfer integrals (in meV), Middle: hole mobilities and Bottom: electron mobilities. Refer to Fig. 4.3 for the symbols corresponding to different simulation conditions

Figure 4 . 8 :

 48 Figure 4.8: Compound 4 mobility along π stack direction (a axis): Hole (Black-Circle) and Electron (Red-Hexagon) mobility computed using Marcus (filled symbols) and MLJ formalisms (empty symbols).

  4.9, 4.10 and 4.11 respectively. Compound 1 exhibits π stacked arrangement along the a axis whereas compounds 2 and 3 exhibit π stacked arrangement along the b axis. However, compounds 1, 2 and 3 show similar mobility of ≈ 0.75 cm 2 /V s along π stack direction. Comparing compounds 2 and 3, both show similar crystallographic packing and relative orientation of molecules as reported in Table-4.1 but slightly differing in relative inter-molecular distances (refer to Appendix 1 for details). This variation in inter-molecular distances leads to a marginal variation in transfer integrals and thereby to variation in absolute values of mobility, with the mobility of compound 2 being larger (along a, b and c axis) than that of compound 3.

Figure 4 . 9 :

 49 Figure 4.9: Compound 1, electron and hole transport. Left(a): ab plane, Middle(b): ac plane and Right (c): bc plane. a b and c represent crystallographic axis direction. Top: Images with non-zero transfer integrals and Bottom: in-plane mobilities

Figure 4 . 11 :Category- 2 :Figure 4 . 12 :Figure 4 . 13 :Figure 4 . 14 :

 4112412413414 Figure 4.10: Compound 2, electron and hole transport. Left(a): ab plane, Middle(b): ac plane and Right (c): bc plane. a b and c represent crystallographic axis direction. Top: Images with non-zero transfer integrals and Bottom: in-plane mobilities.

Figure 4 . 15 :Figure 4 . 16 :

 415416 Figure 4.15: Compound 7, electron and hole transport. Left(a): ab plane, Middle(b): ac plane and Right (c): bc plane. a b and c represent crystallographic axis direction. Top: Images with non-zero transfer integrals and Bottom: in-plane mobilities.

Figure 4 . 17 :

 417 Figure 4.17: Compound 9, electron and hole transport. Left(a): ab plane, Middle(b): ac plane and Right (c): bc plane. a b and c represent crystallographic axis direction. Top: Images with non-zero transfer integrals and Bottom: in-plane mobilities.

Figure 4 . 18 :

 418 Figure 4.18: Compound 10, electron and hole transport. Left(a): ab plane, Middle(b): ac plane and Right (c): bc plane. a b and c represent crystallographic axis direction. Top: Images with non-zero transfer integrals and Bottom: in-plane mobilities.

Figure 5 . 1 :

 51 Figure 5.1: Left: Chemical structure of poly[(9,9-dioctylfluorenyl-2,7-diyl)-co-(4,4 -(N-(4-sec-butylphenyl)diphenylamine)] (TFB, top), and example of an individual oligomer conformation as produced by MD simulations at united atom level; FLU and TPA units are shown in blue and green, respectively while aromatic rings are highlighted in red (bottom). Right: Snapshot of a single sample after equilibration at 300 K and 1 atm, under periodic boundary conditions.

Figure 5 . 2 :

 52 Figure 5.2: Distribution of ZINDO site energies of the four MD samples and corresponding Gaussian fit with a standard deviation σ =50.2 meV and a mean values of -6.82 eV.

Figure 5 . 3 :

 53 Figure 5.3: Average transfer integral as a function of the minimum interatomic carbon-carbon distance between atoms belonging to different oligomers. The rapid decay with distance is well reproduced by an exponential function J 0 e -d/d 0 with J 0 = 6.05 eV and d 0 = 0.37 Å.

Figure 5 . 4 :

 54 Figure 5.4: Average mobility at two different reorganization energies, as a function of electric field and temperature (given in units of K on the left side of the graph): λ =0.1 eV (empty symbols) and 0.3 eV (filled symbols), obtained from R 1 | Q |J Q KMC simulations (see table 5.1 for details). Mobilites are in quantitative and qualitative disagreement with experimental data.

Figure 5 . 5 :

 55 Figure 5.5: Simulated mobility of replicated samples as a function of electric field, with randomized static diagonal disorder with σ = 50.2 meV and fixed transfer integrals, for 3 × 3 replicas of the MD samples (R 3 | S G |J Q conditions). Points correspond to arithmetic averages, dashed lines to logarithmic averages of mobility.

Figure 5 . 6 :

 56 Figure 5.6: Mobility as a function of the magnitude and direction (x, y, z) of the electric field for a single MD sample, obtained with fixed, calculated site energies (filled symbols, R 1 | Q |J Q conditions), with site energies randomized at the beginning of every KMC simulation (empty symbols, R 3 | S G |J Q ), and randomized at every KMC step (gray-shaded symbols, R 3 | G |J Q ).

Figure 5 . 7 :

 57 Figure 5.7: Calculated and experimental mobilities as a function of the electric field and at different temperatures (K units). From the left to the right, in sequence: static diagonal disorder and dynamic off-diagonal disorder (R 3 | S G |J G ); dynamic diagonal disorder and dynamic off-diagonal disorder (R 3 | G |J G ); static diagonal disorder and temperature-dependent dynamic off-diagonal disorder (R 3 | G |J T G). On the rightmost panel are the mobilities from experimental time-of-flight measurements (adapted from references[17] and[15]). The second set of measurements at 295 K, denoted by a star, is reported for completeness; however it comes measurements where an excessive dispersion of the signal was registered, probably due to the presence of impurities acting as traps.[15]

Figure 5 . 8 :

 58 Figure 5.8: Temperature dependence of simulated (R 3 supercell) and experimental zero field mobilities, extrapolated with Eq 5.5 for fields ranging from 150 2 to 450 2 V/cm, and corresponding fitting lines with Eqs. 5.6 (left panel) and 5.8 (right panel). Fitting parameters are reported in Table 5.2. An asterisk * indicates simulation results for λ = 0.2 eV; λ = 0.3 eV otherwise.

Figure 5 . 10 :

 510 Figure 5.10: Effect of simulation parameters on the mobility at 295 K. Filled and empty circles dots correspond to experimental time of flight measurements for references [17] (a) and [15] (b) , respectively.

Figure 6 . 1 :

 61 Figure 6.1: (a) Crystal lattice orientation of rubrene according to Witte's notation [9]; (b) Schematic representation of mechanical strain for uniaxial strain conditions; (c) Top: Graphical representation of strain applied along a axis (π stack direction) with arrows indicating compression or expansion and bottom: Calculated intermolecular distances between first neighbors along the crystalline directions |r a |; (d) Top: Graphical representation of strain applied along b axis (herringbone like packing) with arrows indicating compression or expansion and bottom: Calculated intermolecular distances between first neighbors along the crystalline directions |r b |.

Figure 6 . 2 :

 62 Figure 6.2: Applied strain and the corresponding stress values (in GPa) along the three cell axes, for simulations performed at uniaxial strain conditions.

FigFigure 6 . 3 :

 63 Figure 6.3: (a) Crystal lattice orientation of rubrene according to Wittes notation [9]; (b) Schematic representation of mechanical strain for uniaxial stress conditions; (c) Top : Graphical representation of strain applied along a axis (π stack direction) with arrows indicating compression or expansion and Bottom : Calculated intermolecular distances between first neighbors along the crystalline directions |r a |; (d) Top : Graphical representation of strain applied along b axis (herringbone like packing) with arrows indicating compression or expansion and Bottom : Calculated intermolecular distances between first neighbors along the crystalline directions |r b |.

Figure 6 . 4 :

 64 Figure 6.4: Applied strain and the corresponding stress values (in GPa) along the three cell axes, for simulations performed at uniaxial strain conditions.

Figure 6 . 5 :

 65 Figure 6.5: Normalized distributions of the instantaneous values of the transfer integrals along a, b and c axis : J t a , J t b and J t c , calculated for = 0 (black line), +0.4% (blue line) and -0.4% (red line), in uniaxial strain conditions. From left to right, strain applied along a, b and c directions respectively.

Figure 6 . 6 :

 66 Figure 6.6: Normalized distributions of the instantaneous values of the transfer integrals along a, b and c axis : J t a , J t b and J t c , calculated for = 0 (black line), +0.4% (blue line) and -0.4% (red line), in uniaxial stress conditions. From left to right, strain applied along a, b and c directions respectively

  7 -[ba] and 6.7 -[bb].

Figure 6 . 7 :

 67 Figure 6.7: Calculated relative mobility variations as a function of strain magnitude and direction, calculated along the a (top), b (middle), and c (bottom panels) crystallographic axes. Insets: variation of intermolecular distances |r a |, |r b | and |r c | between dimers along the same direction of the calculated mobility variation, as a function of strain. Error bars were estimated as the difference between the mobility parallel to the cell vectors and with a misalignment of 5 • .

Figure 6 . 8 :

 68 Figure 6.8: Field effect transistors (FET) produced by lamination of single crystals on top of a cantilever beam (cantiFET), a) 3D representation of a cantiFET. The top scheme represents the cantilever at rest, while at the bottom the cantilever is bent, inducing a uniaxial tensile strain along the b-axis of the crystal. b) Architecture of the cantiFETs: a plastic substrate made of polyethylene naphthalate (PEN) is used as support layer of the suspended cantilever. The integrated OFETs consist in an aluminum (AL) bottom gate electrode, a polystyrene (PS) gate dielectric, and gold (Au) source/drain top electrodes. c) Typical forward and reverse transfer curve of a rubrene cantiFET at Vds=-50 V along the a-axis, recorded in absence of strain. d) Real time measurement of drain current along the b-axis in a device submitted to different tensile strains applied along the b-axis at V ds and V gs equal to -50 V, showing the reversible behaviour of the cantiFET response.

Figure 6 . 9 :

 69 Figure 6.9: Experimental relative variation of drain current as a function of strain magnitude and direction of application and measurement. The inset pictures schematize the different setups: a) strain applied along a-axis and drain current measured along the same axis; b) strain applied along b and current measured along a; c) strain applied along a and current measured along b; d) strain applied along b and current measured along b. Horizontal error bars represent the standard deviation on the determination of strain from the deflection profiles, while vertical error bars correspond to the standard deviation of the current calculated over three repeated measurement cycles.

Figure 7 . 1 :

 71 Figure 7.1: Crystalline packing of polymorphs of rubrene, (a) orthorhombic, (b) triclinic and (c) monoclinic. Top: Crystalline packing, Middle: π stacked arrangement and Bottom: top view of π stacked arrangement. Similarity between orthorhombic and triclinic phases for π stacked arrangement can be noticed from (a) & (b) -bottom images. For the sake of clarity, hydrogen atoms are not shown.

Figure 7 . 2 :

 72 Figure 7.2: Crystalline packing of BTBT derivatives: BTBT, C 8 -BTBT, and C 12 -BTBT along with the intermolecular distance and herringbone angles as a function of n in C n -BTBT. For the sake of clarity, hydrogen atoms are not shown.

  789 7.170 14.211 90.0 90.0 90.0 26.452 7.154 13.677 90.0 90.0 90.0 -0.51 % Rubrene: T 7.019 8.543 11.948 93.0 105.5 96.3 6.835 8.311 11.861 93.0 105.2 96.1 -0.57 % Rubrene: M 8.739 10.125 15.635 90.0 90.9 90.0 8.528 9.890 15.448 90.0 90.9 90.0 -0.58 % BTBT 8.102 5.893 11.907 90.0 106.4 90.0 7.589 5.804 11.728 90.0 105.4 90.0 -0.86 % C 8 -BTBT 5.927 7.880 29.18 90.0 92.4 90.0 5.779 7.151 29.129 90.0 92.4 90.0 -1.17 % C 12 -BTBT 5.864 7.740 37.91 90.0 90.58 90.0 5.737 7.072 38.094 90.0 90.5 90.0 -1.01 % Interestingly, for BTBT derivatives (BT BT, C 8 -BT BT, C 12 -BT BT ), PW-DFT optimized structures show a decrease in π stacking distance and decrease in herringbone angle, with the presence and the subsequent increase in chain length (n in C n -BT BT ). The herringbone angle decreases with the increase in chain length, as presented in Fig. 7.3, indicating the variation on packing along the [110] direction from herringbone-like to a slip stack-like arrangement. The validity of decrease in π stacking distance as a function of different vDW-DFT functionals is not verified in this dissertation. However, the decrease in intermolecular distances (along the π stack and transverse to π stack [110] directions) and a slight decrease in herringbone angle can also be observed for C 12 -BTBT in comparison to C 8 -BTBT experimental structures (refer to Fig. 7.

Figure 7 . 3 :

 73 Figure 7.3: Crystalline packing of BTBT derivatives obtained after PW-DFT oprimization : BTBT, C 8 -BTBT, and C 12 -BTBT along with the intermolecular distance and herringbone angels as a function of n in C n -BT BT . For the sake of clarity, hydrogen atoms are not shown.

Figure 7 . 4 :Figure 7 . 5 :

 7475 Figure 7.4: Band dispersion of rubrene: (a) orthorhombic , (b) triclinic and (c) monoclinic, along with k-path used to compute the band dispersion (middle). Band dispersion in orthorhombic phase is obtained from optimized structure with reduced symmetry : hexagonal representation.

Figure 7 . 6 :

 76 Figure 7.6: Transfer integrals (left) and effective masses (right) obtained for PW-DFT optimized structure, without any applied strain for triclinic rubrene.

Figure 7 . 7 :

 77 Figure 7.7: Strain mobility trends in triclinic rubrene along a-axis (π stack direction), for strain applied along a, b and c axes. Left: absolute values of transfer integral (top) and effective mass (bottom) as a function of applied strain. Right: relative variation of mobility obtained from transfer integrals (top) and effective masses (bottom).

Figure 7 . 8 :

 78 Figure 7.8: Strain mobility trends in triclinic rubrene along the ab direction [110], for strain applied along a, b and c axes. Left: absolute values of transfer integral (top) and effective mass (bottom) as a function of the applied strain. Right: relative variation of mobility obtained from transfer integrals (top) and effective masses (bottom).

Figure 7 . 9 :

 79 Figure 7.9: Rubrene: monoclinic. Transfer integrals (left) and effective mass (right) obtained for PW-DFT optimized structure, without any applied strain.

Figure 7 . 10 :

 710 Figure 7.10: Strain mobility trends in monoclinic rubrene along a-axis (π stack direction), for strain applied along a, b and c axes. Left : absolute values of transfer integral (top) and effective mass (bottom) as a function of applied strain. Right: relative variation of mobility obtained from transfer integrals (top) and effective masses (bottom).

Figure 7 . 11 :

 711 Figure 7.11: BTBT. Transfer integrals (left) and effective masses (right) obtained for PW-DFT optimized structure,without strain.

Figure 7 . 12 :

 712 Figure 7.12: Strain mobility trends in BTBT crystal along b-axis (π stack direction), for strain applied along a and b axes. Left: absolute values of transfer integral (top) and effective mass (bottom) as a function of applied strain. Right: relative variation of mobility obtained from transfer integrals (top) and effective masses (bottom).

Figure 7 . 13 :

 713 Figure 7.13: Strain mobility trends in BTBT crystal along ab (110) direction, for strain applied along a and b axes. Left: absolute values of transfer integral (top) and effective mass (bottom) as a function of applied strain. Right: relative variation of mobility obtained from transfer integrals (top) and effective masses (bottom). 150

C 8 Figure 7 .Figure 7 . 15 :

 87715 Figure 7.14: C8-BTBT. Transfer integrals (left) and effective masses (right) obtained for PW-DFT optimized structure, without strain.

Figure 7 . 16 :

 716 Figure 7.16: Strain mobility trends in C 8 -BTBT crystal along ab direction (110), for strain applied along a and b axes. Left: absolute values of transfer integral (top) and effective mass (bottom) as a function of applied strain. Right:relative variation of mobility obtained from transfer integrals (top) and effective masses (bottom).

Figure 7 . 17 :

 717 Figure 7.17: Relative valence band shift (∆V BS) for (a) Top-Left -rubrene triclinic, (b) Top-Right -rubrene monoclinic for strains applied along a, b and c crystallographic axis and (c) Bottom-Left -BTBT , (d) Bottom-Right -C 8 -BTBT for strains applied along a and b crystallographic axis.

Figure

  Figure 8.1: pK a values computed with DFT-PBE/def2-TZVP/CPCM in comparison to experimental values. Filed circles: correspond to theoretical/calculated values and empty squares: to values corrected with Eq. 8.3. Lines is a guide to the eye to show linear correction between experimental and calculated values, obtained using Eq. 8.3.

8. 3

 3 Protonation state: Hydrogen bonding & XPS spectra DFT calculations are performed on Dimer and Trimer of STFSI units at the PBE/dft2-TZVP level of theory, treating water implicitly using CPCM continuum solvation model. DFT optimized structures show the presence of π-stacked-like or at least eclipsed phenyl units depending on the protonation of amide nitrogens. π stackedlike conformation is stabilized by intra-chain hydrogen bonding between N H group of a given STFSI neutral monomer unit and the oxygen atoms of the S = O group belonging to neighboring STFSI units. The presence of hydrogen bondings in DFT optimized geometries is indicated with red dashed lines in Fig. 8.2.

Figure 8 . 2 :

 82 Figure 8.2: Presence of hydrogen bonding and π stacked regions in N H -N H (left) and N H -N H -N H (right) oligomers, as observed in DFT calculations/optimized structures.

Figure 8 . 3 :Figure 8 . 4 :

 8384 Figure 8.3: IR spectra of STFSI monomer. Left(a): Superpositioned spectra as the percentage contribution of protonated and de-protonated spectra, Middle(b): Comparison with PSTFSI -KOH treated experimental spectra and Right(c): Comparison with PSTFSI -acid treated, experimental spectra.

3 Figure 8 . 5 :

 385 Figure 8.5: IR (left) and Raman spectra (right) of STFSI trimer as a function of protonation state of nitrogen. Legend correspond to protonation state of nitrogen in the timers. STFSI monomer units with protonated and de-protonated nitrogens are represented by N H and N , respectively.

Figure 8 . 6 :

 86 Figure 8.6: IR (left) and Raman spectra (right) of STFSI trimer with terminal deprotonated nitrogen, in comparison to experimental spectra for PSTFSI.

Figure 8 . 7 :

 87 Figure 8.7: Final conformations of oligomers after MD simulations. STFSI units with de-protonated nitrogens are shown in red while units with protonated nitrogens in green. Regions of π stacking can be noticed when protonated nitrogen states are present.

Figure 8 . 8 :

 88 Figure 8.8: Dynamic hydrogen bonds for PSTFSI -(N ) 10 conformation with all deprotonated states: Hydrogen bonding between surrounding water molecules and deprotonated nitrogen atoms (n ODW •••N -)

Figure 8 . 9 :

 89 Figure 8.9: Dynamic hydrogen bonds for the PSTFSI -(N H -N ) 5 conformation with alternating protonated nitrogen states. See Table. 8.3 for details.

EnergyFigure 8 . 11 :Figure 8 . 12 :

 811812 Figure 8.11: Free energy computed for trimers as a function of distance between neighboring nitrogen and sulpher atoms.

  1: ab plane (b) 1: ac plane 1: 6: ab plane (b) 6: ac plane 6: bc plane A:2

Figure A2- 2 :Figure A2- 5 :

 25 Figure A2-2: Top: Definition of charge-atom types for the two TFB co-monomers FLU (left) and TPA (right). Bottom: ESP charges on the central co-monomers of each atom types (red crosses) and respective average values (black squares).

Figure A2- 6 :

 6 Figure A2-6: Distribution of the number of neighbours with transfer integral greater than 0.1 meV along the TFB oligomer chain. Odd and even numbers correspond to FLU and TPA units, respectively.

Figure A3- 5 :

 5 Figure A3-5: a) Visual representation of ϕ, θ and χ. The rotation direction with respect to the symmetry axes (dashed lines) is shown by the black arrows. b) Variation of J a as a function of ϕ, θ and χ, calculated from a dimer extracted from the crystal unit cell in the absence of strain. c) Variation of ϕ, θ and χ as a function of strain, extracted from MD simulations.

Figure A4- 1 :

 1 Figure A4-1: Strain mobility trends in Rubrene : Triclinic crystal along a axis (π stack direction), for strain applied along a, b and c axes. Left : Absolute values of transfer integrals (top) and effective mass (bottom) as a function of applied strain and Right : relative values of transfer integrals (top) and effective mass (bottom).

Figure A4- 2 :Figure A4- 3 :

 23 Figure A4-2: Strain mobility trends in Rubrene : Triclinic crystal along ab direction (110), for strain applied along a, b and c axes. Left : Absolute values of transfer integrals (top) and effective mass (bottom) as a function of the applied strain and Right : relative values of transfer integrals (top) and effective mass (bottom).

Figure A4- 4 :

 4 Figure A4-4: Strain mobility trends in BTBT crystal along b axis (π stack direction), for strain applied along a and b axes. Left : Absolute values of transfer integrals (top) and effective mass (bottom) as a function of applied strain and Right : relative values of transfer integrals (top) and effective mass (bottom).

Figure A4- 5 :Figure A4- 6 :

 56 Figure A4-5: Strain mobility trends in BTBT crystal along ab (110) direction, for strain applied along a and b axes. Left : Absolute values of transfer integrals (top) and effective mass (bottom) as a function of applied strain and Right : relative values of transfer integrals (top) and effective mass (bottom).

C 8 Figure A4- 7 :Figure A4- 8 :

 878 Figure A4-7: Strain mobility trends in C8-BTBT crystal along a axis (π stack direction), for strain applied along a and b axes. Left : Absolute values of transfer integrals (top) and effective mass (bottom) as a function of applied strain and Right : relative values of transfer integrals (top) and effective mass (bottom).

  Band Transport . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 Kinetic Monte-Carlo . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 Molecular Mechanics . . . . . . . . . . . . . . . . . . . . . . . . . . .

1 Introduction 1.1 Materials and Devices . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Modeling: State of the art . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Challenges & Motivation . . . . . . . . . . . . . . . . . . . . . . . . . 2 Theoretical Background 2.1 Charge transport : Methodology & Models . . . . . . . . . . . . . . . 2.1.1 Hopping picture of charge transport . . . . . . . . . . . . . . . 2.1.2 Hopping Mechanism & Models . . . . . . . . . . . . . . . . . . 2.2 Reorganization Energy (λ) . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Electron-Phonon coupling & Huang-Rhys parameters . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Transfer Integrals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 ix CHAPTER 1

Table 4 .

 4

	Molecule	IP (eV)	EA (eV)	HOMO (eV) LUMO (eV) Ref.
		Exp Calc. Exp. Calc.	Calc.	Calc.
	1	-	6.074	-	1.829	-5.116	-2.954	-
	2	5.8 6.356 4.0 2.216	-5.291	-3.163	[13]
	3	-	6.251	-	2.075	-5.426	-3.310	-
	4	5.7 6.161 3.8 2.110	-5.264	-3.168	[13]
	5	5.7 6.066 4.3 2.012	-5.171	-3.073	[15]
	6	5.8 6.174 4.4 2.085	-5.281	-3.160	[15]
	7	5.9 6.269 4.5 2.200	-5.372	-3.252	[15]
	8	5.4 6.120 4.3 1.911	-5.109	-2.995	[15]
	9	-	6.005	-	1.908	-5.098	-2.944	-
	10	-	6.140	-	2.091	-5.005	-3.106	-

2: Caculated ionization potential (IP) and electron affinity (EA) along with available experimental values in eV (with the corresponding reference) and calculated HOMO and LUMO energies (eV) of the molecules considered in this dissertation. Internal reorganization energies and Huang-Rhys parameters corresponding to Neutral (N), Charged(C) and Total (sum of neutral and charged states) for Holes and Electrons, are reported in

Table 4.3 and 4.4 respectively. The internal reorganization energies and effective Huang-Rhys parameters for electron transport are henceforth referred to as λ e i , hω e E and S e E , respectively, and those for hole transport as λ h i , hω h E and S h E .

Table 4 . 3 :

 43 Huang-Rhys parameters (hω h and S h ) of the molecules for neutral (N) and positively charged state (C) computed following Ref.[17] for hole transport, λ h

	Molecule	λ h i	S h	S h E	hω h	hω h E
		eV	N	C Total	N	C Total
	1	0.036 0.068 0.068 0.137 0.262 0.263 0.261
	2	0.046 0.107 0.103 0.211 0.213 0.223 0.217
	3	0.040 0.089 0.085 0.174 0.224 0.235 0.229
	4	0.064 0.325 0.319 0.644 0.098 0.100 0.099
	5	0.061 0.258 0.270 0.529 0.118 0.113 0.115
	6	0.058 0.135 0.133 0.268 0.215 0.217 0.215
	7	0.064 0.275 0.269 0.544 0.116 0.119 0.117
	8	0.063 0.250 0.245 0.495 0.126 0.128 0.127
	9	0.054 0.097 0.100 0.197 0.278 0.269 0.272
	10	0.051 0.103 0.103 0.206 0.248 0.247 0.246

i is the hole reorganization energy (internal). λ h i and hω h are in units of eV.

Table 4 .

 4 4: Huang-Rhys parameters (hω e and S e ) of the molecules for neutral (N) and negatively charged state (C) computed following Ref.[17] for electron transport, λ e i is the electron reorganization energy (internal). λ e i and hω e are in units of eV.

	Molecule	λ e i	S e	S e E	hω e	hω e E
		eV	N	C Total	N	C Total
	1	0.206 0.774 0.785 1.559 0.133 0.131 0.132
	2	0.216 0.948 0.970 1.917 0.114 0.111 0.112
	3	0.220 0.988 0.936 1.924 0.111 0.1178 0.114
	4	0.240 1.125 1.195 2.319 0.107 0.100 0.103
	5	0.239 1.122 1.144 2.265 0.107 0.1051 0.105
	6	0.253 1.003 0.961 1.963 0.126 0.1321 0.128
	7	0.244 1.012 0.984 1.996 0.120 0.1243 0.122
	8	0.236 1.049 1.071 2.121 0.112 0.1105 0.111
	9	0.272 0.718 0.673 1.135 0.190 0.202 0.195
	10	0.228 0.503 0.507 1.010 0.227 0.225 0.225

Table 4 .

 4

	Mol. Crystal direction J h J e	Mol. Crystal direction J h J e
	1	(1, 0, 0)	42 31	5	(1,0,0)	40 41
		(0.5, 0.5, 0.5)	08 18		(0.5,0.3,1)	11 12
		(0,1,0)	03 02		(0.3,0.7,1)	12 05
	2	(0, 1, 0)	33 29	6	(1,0,0)	19 23
		(1, 1, 0)	11 12		(0,1,0)	01 01
		(0.7,0.3,1)	05 08		(1,0.5,1)	08 04
		(0.2,0.3,1)	02 18	7	(1,0,0)	37 37
		(1,0,0)	01 04		(0.5,0,1)	03 14
	3	(0,1,0)	33 30		(0.25,1,0.5)	09 10
		(1, 1, 0)	13 11	8	(0.5,0,1)	54 40
		(0.7,0.3,1)	06 19		(0.4,0,1)	03 17
		(0.2,0.3,1)	03 18		(0.25,1,0.5)	02 05
		(1,0,0)	01 04	9
	4	(1,0,0)	29 43
		(0.5,0.3,1)	11 12
		(0.3,0.7,1)	10 05

5: Transfer integrals (absolute value of J, in units of meV) between a reference molecule and its first neighbors (labeled in crystal units) for holes (h) and electrons (e) in the investigated phtalocyanine crystals.

  Gaussian disorder generated at every hopping event during the KMC simulation run. x is the random number drawn from the Gaussian distribution centered at J Q and standard deviation of σ J .Effect of reorganization energy and nature of energetic disorder are presented in Fig.4.3. For the same value of the internal reorganization energy, MLJ formalism (filled) systematically gives lower mobility with respect to Marcus formalism (empty symbols with dashed lines), attributed to the contribution of Huang-Rhys parameters, eliciting the impact of coupling the quantum treatment of effective vibrational mode with the charge transport characteristics. Ref. Fig.4.3. However, a good qualitative agreement between both Marcus and MLJ formalisms is obtained.

Table 5 . 1 :

 51 List of symbols utilized for defining the KMC simulations conditions.

	Symbol Definition
	R 1	Original four MD samples (125 sites each)
	R 2	2×2×2 replica of MD samples (1000 sites)
	R 3	3×3×3 replica of MD samples (3375 sites)

Q

HOMO energies calculated at ZINDO level at the MD geometry, static in time with σ =50.2 meV G HOMO energies drawn at every KMC step from a Gaussian distribution with σ =50.2 meV (dynamic diagonal disorder) S G

Table 6 .

 6 1: Comparison between experimental density and unit cell axes and MD simulation results obtained for a 4 × 2 × 1 bulk rubrene supercell at 300K. Unit cell angles were held fixed to the experimental value of 90 degrees.

	Term	a( Å) b( Å) c( Å) Density (gm/cm 3 )
	Experiment	7.17 14.43 26.81	1.27
	This work (300K) 7.32 14.15 27.12	1.26
	% Difference	2.0	-2.1	1.1	-0.7

Table 6 .

 6 4: Summary of applied strain and the corresponding stress values (in GPa) along the three cell axes, together with the dimensions of the unit cell ( Å), for uniaxial stress simulations. . .010 7.311 14.162 27.115 0.006 7.323 14.136 27.123 0.012 7.318 14.157 27.092 -0.2 0.018 7.304 14.174 27.112 0.011 7.326 14.121 27.129 0.024 7.317 14.164 27.064 -0.3 0.027 7.297 14.187 27.110 0.016 7.33 14.107 27.135 0.036 7.316 14.171 27.036 -0.4 0.036 7.290 14.199 27.107 0.021 7.333 14.093 27.140 0.049 7.315 14.177 27.008

	(%)	σ aa	a	b	c	σ bb	a	b	c	σ cc	a	b	c
	0.4 -0.036 7.348 14.100 27.128 -0.021 7.304 14.206 27.095 -0.049 7.323 14.122 27.226
	0.3 -0.027 7.341 14.113 27.125 -0.016 7.308 14.192 27.100 -0.036 7.322 14.129 27.200
	0.2 -0.018 7.334 14.125 27.123 -0.011 7.311 14.178 27.106 -0.024 7.321 14.136 27.172
	0.1	0.010 7.326 14.137 27.120 -0.006 7.315 14.164 27.112 -0.012 7.32 14.143 27.144
	-0.1 0

Table 6 .

 6 5: Parameter m ij obtained from the linear fit of experimental mobilities for uniaxial strain and uniaxial stress conditions, obtained with Eq. 6.4.

	m ij MD uniaxial strain MD uniaxial stress
		aa	bb	cc	aa	bb	cc
	µ a -12.4 -16.8	8.1	3.3 -12.5 14.7
	µ b	-4.2 -19.2 -16.8 11.6 -13.2 -7.8
	µ c 18.6	0.6 -148.9 3.2 24.8 -147.5

Table 7 .

 7 2: Crystallographic parameters for BTBT, C 8 -BTBT and C 12 -BTBT derivatives. Space group is represented in Hermann-Mauguin notation and Z corresponds to the number of formula units per unit cell.

	Derivative	Cell paramaters	Reference

Table 8 .

 8 8.1: pK a values computed with DFT-PBE/def2-TZVP/CPCM in comparison to experimental values. Filed circles: correspond to theoretical/calculated values and empty squares: to values corrected with Eq. 8.3. Lines is a guide to the eye to show linear correction between experimental and calculated values, obtained using Eq. 8.3.From the linear regression, an R 2 of 0.95 (with S = 0.946 and C = 0.824) and mean standard error between the experimental and corrected pK a values of 0.1 pK a units were obtained, indicating that the scheme is appropriate, and it can be used to estimate the pK a values with good acurracy. (See Fig.8.1 and Table:8.1 for details) 1: Experimental vs calculated pK a obtained with Eq. 8.3.

	Molecule

Table 8 .

 8 2: Calculated pK a (after benchmark correction) as a function of oligomer length. pK a 1 refer to pK a when one amide hydrogen is removed (first de-protonation), pK a 2 when second amide hydrogen is removed from pK a 1 state (2 nd de-protonation of dimer and trimer) and pK a 3 when amide hydrogen is removed from pK a 2 state (3 rd de-protonation of trimer).Oligomer pK a 1 pK a 2 pK a 3

	Monomer 0.34	-	-
	Dimer	0.37	1.63	-
	Trimer	0.41	0.86 1.72

  Table. 8.2 for dimer and trimer are obtained by removing the amide hydrogens starting from one end of the oligomer (end STFSI unit). These values may change (slightly) if, the sequence of removing the hydrogens is changed, i.e. if the first amide hydrogen is removed from the central STFSI unit, subsequently followed by the removal of hydrogens from the end STFSI units.

  • ) are reported in Table. 8.4, and described and detailed in Table. 8.3. Average number of intra-chain hydrogen bonds between protonated nitrogen atom of any given STFSI unit and deprotonated nitrogen atom of neighboring STFSI unit. n N H•••ODW Average number of hydrogen bonds between protonated nitrogen atom of any given STFSI unit and oxygen atom from explicit water molecules.n ODW •••N -Average number of hydrogen bonds between de-protonated nitrogen atom of any given STFSI unit and oxygen atom from explicit water molecules.

		Table 8.3: Nomenclature employed in the analysis.
	Symbol	Definition
	r π	Average distance between center of mass of neighboring
		phenyl units
	r E-E	End-to-End distance of the oligomer computed as the average
		distance of terminal carbon atoms.
	n N H•••0	Average number of intra-chain hydrogen bonds between pro-
		tonated nitrogen atom of any given STFSI unit and oxygen
		atom of neighboring STFSI unit.
	n N H•••N -
	Table 8.4: Parameters extracted from MD simulations of STFSI decamers. Data
	corresponding to hydrogen bonds represent average number of hydrogen bonds formed
	per STFSI unit per MD simulation time step.
	Type	r π

  8.10: Dynamic hydrogen bonds for the PSTFSI -(N H) 10 , all protonated nitrogen states conformation. See Table. 8.3 for details.

Table A1 -

 A1 1: Compound 1 (ϕ 1 in fig. A2-1) of 36 degrees. The triphenylamine block presents coplanar central NCCC atoms and a propeller-like arrangement of phenyl rigs with a CNCC dihedral angle (ϕ 2 in fig. A2-1) of 40 degrees

	Table A1-2: Compound 2 Table A1-4: Compound 4 Table A1-6: Compound 6 Table A1-8: Compound 8 Table A1-10: Compound 10 triphenylamine dihedral angle
	Neigh. Neigh. Neigh. Neigh. Neigh. Neigh.	J h J h J h J h J h J h	J e J e J e J e J e J e	PBC vector ( Å) PBC vector ( Å) PBC vector ( Å) PBC vector ( Å) PBC vector ( Å) PBC vector ( Å)	PBC scalar ( Å) PBC scalar ( Å) PBC scalar ( Å) PBC scalar ( Å) PBC scalar ( Å) PBC scalar ( Å)
		(meV) (meV) (meV) (meV) (meV) (meV) (meV) (meV) (meV) (meV) (meV) (meV)	x x x x x x	y	y y y y y	z z z z	z z	x 2 + y 2 + z 2 x 2 + y 2 + z 2 x 2 + y 2 + z 2 x 2 + y 2 + z 2 x 2 + y 2 + z 2 x 2 + y 2 + z 2
	1	1.22 12.89 0.26 0.02 53.68 1.62	3.26 -8.615 -13.344 11 -7.913 -8.686 0.4 -12.759 3.022 0.03 -8.448 -13.602 40.06 -4.822 0 -8.814 0 0 10.41 0 0.72 -10.919 0 -7.34	15.883 11.75 16.742 16.012 10.047 13.157
	2	0.01 0.16 28.86 18.37 0.12 4.84	0.01 -7.979 0.72 -6.124 -5.743 10.411 0 -10.367 42.98 -8.342 0 0 23.11 -8.448 0 0 0.58 -9.64 0 0 3.56 -12.643 0 0	13.082 13.374 8.342 8.448 9.64 12.643
	3	41.71 0.31 10.65 0.04 2.42 0.7	30.89 -8.615 4.1 -7.365 12.04 -4.417 0.09 -8.448 13.602 0 0 3.022 17.39 4.817 0 -8.814 0 0 10.41 0 2 -14.367 0 7.34	8.615 7.365 11.705 16.012 10.044 16.133
	4	2.49 4.78 0.44 1.87 2.42 121.36	1.54 6.27 -5.576 2.944 10.411 0 -13.344 0 0.63 -11.047 7.224 -10.411 1.44 0 -13.602 0 17.4 -4.817 0 8.814 86.26 1.724 0 -7.34	13.344 12.172 16.811 13.602 10.044 7.54
	5	2.49 0.03 0.31 1.87 0.12 121.36	1.54 0.03 -5.028 11.63 10.411 0 13.344 0.44 -7.123 10.247 1.44 0 13.602 0 0.58 9.64 0 0 86.26 -1.724 0 7.34	0 0	13.344 16.399 12.48 13.602 9.64 7.54
	6	41.65 1.21 1.67 0.04 53.68 0.7	30.93 1.17 -2.337 -11.63 -10.411 8.614 0 0 0.38 -1.218 -10.248 0 0.09 8.448 -13.602 0 40.06 4.822 0 8.814 2 14.367 0 -7.34	8.614 15.783 10.32 16.012 10.047 16.133
	7	0.01 33.01 9.33 18.37 0.03 4.84	0.01 29.48 -0.548 -8.686 7.978 5.3 2.706 -7.225 0 10.366 0 10.41 23.11 8.448 0 0 0.02 -7.231 -9.886 -4.407 3.56 12.643 0 0	13.081 8.703 12.957 8.448 13.017 12.643
	8	1.22 2.74 0.05 0.02 0.03 1.62	3.26 18.46 -1.789 -2.943 -10.411 8.614 13.344 0 0.43 -3.924 -3.023 -10.411 0.03 8.448 13.602 0 0.02 -12.048 -9.886 4.407 0.72 10.919 0 7.34	15.883 10.966 11.529 16.012 16.196 13.157
	9	1.28 2.74 0.05 7.06 0.03 0.14	1.33 -12.604 -6.672 -5.183 18.48 1.789 2.944 10.411 0.43 3.924 3.022 10.41 3.58 -7.46 -6.801 -7.555 0.02 -7.231 9.885 -4.407 0.79 -11.781 -9.794 -3.67	15.174 10.966 11.528 12.609 13.017 15.754
	10	1.28 33.07 9.38 7.06 0.03 0.01	1.33 -12.604 29.45 0.548 8.687 6.672 -5.183 0 5.31 -2.706 7.224 -10.411 3.56 -7.46 6.801 -7.555 0.02 -12.048 9.885 4.407 0.01 -13.505 -9.794 3.67	15.174 8.704 12.958 12.609 16.196 17.081
	11	7.5 1.21 1.67 0.99 1 0.14	18.31 -3.989 -6.672 -5.183 1.17 2.337 11.63 10.411 0.4 1.218 10.247 0 0.88 0.988 -6.801 -7.555 4.81 2.409 -9.886 -4.407 0.79 -11.781 9.793 -3.67	9.343 15.783 10.319 10.213 11.088 15.754
	12	0.07 0.03 0.27 1 1 0.01	0.04 -4.625 -6.672 0.03 5.028 -11.63 -10.411 5.183 0.44 7.124 -10.248 0 0.89 -0.988 -6.801 7.555 4.83 -2.408 -9.886 4.407 0.01 -13.505 9.793 3.67	9.632 16.399 12.481 10.213 11.088 17.081
	13	7.55 4.78 0.44 0.99 1.01 0.54	18.27 -3.989 6.27 5.576 -2.943 -10.411 6.672 -5.183 0.63 11.048 -7.225 10.41 0.89 0.988 6.801 -7.555 4.83 2.409 9.885 -4.407 1.8 0.862 -9.794 -3.67	9.343 12.171 16.812 10.213 11.088 10.494
	14	0.07 0.31 10.65 1 0.99 0.55	0.04 -4.625 4.1 7.365 12.04 4.418 -3.023 -10.411 6.672 5.183 0 0 0.89 -0.988 6.801 7.555 4.82 -2.408 9.885 4.407 1.8 -0.862 -9.794 3.67	9.632 7.365 11.707 10.213 11.088 10.494
	15	0.07 0.16 28.86 7.07 0.03 0.54	0.04 0.72 6.124 5.744 -10.411 4.625 -6.672 -5.183 42.99 8.342 0 0 3.59 7.461 -6.801 7.555 0.02 12.049 -9.886 -4.407 1.8 0.862 9.793 -3.67	9.632 13.375 8.342 12.61 16.196 10.494
	16	7.45 12.89 0.26 7.13 0.03 0.55	18.26 10.98 7.913 8.687 3.989 -6.672 0.42 12.76 -3.023 -10.411 5.183 0 3.57 7.461 6.801 7.555 0.02 7.231 -9.886 4.407 1.81 -0.862 9.793 3.67	9.343 11.751 16.744 12.609 13.017 10.494
		0.07 0.03 0.01	0.04 0.02 12.049 9.885 -4.407 4.625 6.672 -5.183 0.01 13.505 9.793 -3.67	9.632 16.196 17.081
		7.39 0.03 0.14	18.28 0.02 0.79 11.781 9.793 3.67 3.989 6.672 5.183 7.231 9.885 4.407	9.343 13.017 15.753
		1.28	1.34 12.604 -6.672 Table A1-3: Compound 3 5.183 Table A1-5: Compound 5 Table A1-7: Compound 7	15.173
		1.29	1.33 12.604	6.672	5.183	15.173
	Neigh. Neigh. Neigh.	J h J h J h	J e J e J e Table A1-9: Compound 9 PBC vector ( Å) PBC vector ( Å) PBC vector ( Å)	PBC scalar ( Å) PBC scalar ( Å) PBC scalar ( Å)
		(meV) (meV) (meV) (meV) (meV) (meV)	x x x	y y y	z	z z	x 2 + y 2 + z 2 x 2 + y 2 + z 2 x 2 + y 2 + z 2
	1 Neigh.	11.58 0.25 2.23 J h	10.45 -8.049 -8.649 0.47 -12.812 2.881 10.402 0 13.54 -5.196 0 -9.561 J e PBC vector ( Å)	11.815 16.753 10.882 PBC scalar ( Å)
	2	0.76 39.81 36.2 (meV) (meV) 0.55 -6.193 41.5 -8.417 36.56 -10.21 x	-5.78 10.506 0 0 0 0 y z	13.496 8.417 10.21 x 2 + y 2 + z 2
	3	0.26 10.42 1.68 3.42	4.32 -7.363 12.32 -4.394 0.79 5.013 7.63 -5.235	0 2.881 10.402 0 0 -9.561 0 -10.36	7.363 11.654 10.796 11.608
	4	5.87 0.5 1.68 18.59	7.56 -5.508 0.81 -10.914 0.8 -5.014 25.14 -10.256	2.868 10.506 7.172 -10.403 0 9.562 0 0	12.204 16.697 10.797 10.256
	5	0.01 0.68 36.2 1.52	0.02 -4.822 11.517 10.506 0.54 -6.891 10.053 0 36.56 10.21 0 0 0.86 5.021 0 -10.36	16.318 12.188 10.21 11.513
	6	1.4 2.01 2.22 1.53	1.38 -2.541 -11.518 -10.507 0.44 -1.525 -10.053 0 13.5 5.196 0 9.562 0.86 -5.021 0 10.359	15.796 10.168 10.883 11.512
	7	35.55 11.03 8.33 18.53	33.8 -0.686 -8.649 5.47 2.498 -7.172 10.402 0 9.9 -2.598 -8.849 -4.781 25.14 10.257 0 0	8.676 12.879 10.388 10.257
	8	1.95 11.03 0.18 3.41	18.87 -1.856 -2.869 -10.507 5.47 -2.497 7.172 -10.403 0.1 -7.612 -8.849 4.781 7.64 5.236 0 10.359	11.049 12.88 12.614 11.607
	9	1.95 2.01 8.34 7.01	18.88 1.855 0.44 1.526 10.053 2.868 10.506 0 9.9 -2.598 8.849 -4.781 9.89 -2.733 -8.283 -5.355	11.047 10.168 10.388 10.235
	10	35.41 0.68 0.18 0.61	33.78 0.685 0.54 6.892 -10.053 8.649 0.1 -7.612 8.849 4.781 0.14 -7.754 -8.283 5.005 0 0	8.676 12.189 12.614 12.401
	11	1.4 0.5 0.18 7.02	1.38 2.541 11.517 10.506 0.81 10.915 -7.172 10.402 0.1 7.612 -8.849 -4.781 9.91 -2.733 8.283 -5.355	15.795 16.697 12.613 10.235
	12	0.01 10.42 8.35 0.61	0.02 4.822 -11.518 -10.507 12.3 4.395 -2.881 -10.403 9.9 2.598 -8.849 4.781 0.14 -7.754 8.283 5.005	16.319 11.655 10.388 12.401
	13	5.87 39.87 0.18 0.68	7.55 5.507 -2.869 -10.507 41.44 8.418 0 0 0.1 7.612 8.849 -4.781 0.11 7.523 -8.283 -5.355	12.205 8.418 12.614 12.404
	14	0.26 0.21 8.3 8.02	4.32 7.363 0.47 12.813 -2.881 -10.403 0 0 9.88 2.598 8.849 4.781 11.35 2.502 -8.283 5.005	7.363 16.754 10.388 9.996
		0.76 0.68	0.56 6.192 0.11 7.523 8.283 -5.355 5.78 -10.507	13.496 12.405
		11.63 8.04	10.45 8.048 11.34 2.502 8.283 5.005 8.649 0	11.814 9.996
								A:7 A:8 A:9 A:10 A:11 A:12

4 6 8 10 12 14 16 Distance (Å)

ACKNOWLEDGEMENTS

"Elite and great are many, and I salute to thee..."

the charge transport is improved compared to a system consisting of short and randomly oriented PEDOT oligomers. Starting, instead, from all de-protonated PSTFSI (starting from potassium/lithium salt form: PSTFSI-K/Li, PEDOT could be aligned along the STFSI units that are drifted apart due to higher repulsive forces between the neighboring STFSI units, which in-turn can result in reduction of linearly arranged percolation pathways, thereby leading to a decrease in conductivity of PE-DOT:PSTFSI complexes. Indeed, subsequent experimental observations based on this theoretical study, confirmed that, when the acidic form of the PSTFSI-H polyelectrolyte was used instead of its potassium salt PSTFSI-K, an average conductivity of 327 S/cm for PEDOT:PSTFSI-H was obtained with respect to 135 S/cm for PEDOT:PSTFSI-K. This equals to a ≈ 2.4 fold increase as compared to PEDOT:PSTFSI-K, demonstrating that the inherent characteristics of polyelectrolyte, like its protonation-dependent conformation are decisive for the conductive properties of PEDOT:polyelectrolyte complexes.

Transfer Integrals

Transfer integrals of compounds 1 to 10, between first neighbors computed using ORCA software and employing projection method (refer Chapter-02: Theoretical Methodology for details) at the DFT/B3LYP/6-31G(d) level of theory by replicating the experimental crystal cell and using a spherical cut off of 20 Å are presented in Tables A1-1 to A1-10, wherein J h and J e correspond to transfer integrals for hole and electron transport respectively (in units of meV). Also reported, in Tables A1-1 to A1-10, are the PBC vectors along the three crystallographic directions and the corresponding PBC scalar distances.

A:6

: ENERGETIC FLUCTUATIONS IN AMORPHOUS SEMICONDUCTING POLYMERS

Development of a united atom force field for TFB TFB oligomers are modeled with a united atoms force field (ff). The ff is based on the AMBER parameterization [1] (AMBER united atom [2] for alkyl chains), and is complemented with quantum chemical calculations for atomic charges, torsional potentials and small adjustments of equilibrium values of bond distances and angles. Torsional disorder along polymer chains is very important since it is expected to affect the localization of excess charges on polymer chains. General purpose ff, such as AMBER, cannot adequately reproduce the torsional potential of the three dihedrals sketched in figure A2-1, that are here parameterized from DFT calculations. Quantomechanic torsional potential are plotted with black squares (and red circles for ϕ 2 ) in Fig. A2-3. Torsional potentials are obtained with relaxed geometry scans of the relevant dihedral angle, calculated at B3LYP/6-31G** level of theory. Note that the ϕ 2 potential depends on the scanning direction as already reported in the literature [4,5]. This dependence on the scanning direction results from the steric hindrance among phenyl rings arranged in a propeller-like structure.

The ab initio torsional potential is then introduced in the FF following a wellestablished procedure described in previous works [7] and here shortly summarized.

The torsional potential of the FF is the sum of an explicit dihedral potential plus an implicit contribution from other interactions in the ff, i.e. U ff (φ) = U Dih (φ) + U i (φ).

To correctly introduce the quantomechanic potential in the ff, we first measured the implicit potential, U i (φ), from a separate simulation with U Dih (φ) set to 0, performed at 400 K on a fluorene-triphenylamine-fluorene molecule with an inert gas to facilitate energy redistribution trough the sample (40 Ar atoms in a cubic box of side 50 Å). Unlike previous works, where the torsional potential was obtained from the dihedral probability distribution, here we use the adaptive biasing force (ABF) methods [8], allowing a more efficient sampling of the potential energy surface. Once U i (φ) is known, the explicit dihedral potential, U Dih (φ), is obtained by fitting the difference between the ab initio potential and U i (φ) with a series of cosines. The final FF parameters are reported in table A2-1. The reparametrized ff torsional potentials, shown as green lines in figure A2-4, carefully reproduce the ab initio potentials. The force field and the DFT optimized geometries are consistent within an accuracy of 0.005 Å for bond lengths, 3 degrees for bending angles and 5 degrees for dihedral angles. The curves reveal the absence of crystalline order in the sample, with the sharp peak at about 15 Åfor triphenylamine and the one ad about 10 Åfor terminal FLU corresponding to intra-chain correlation due to TPA units belonging to the same oligomer. Right: radial distribution function between aromatic carbons with on implicit hydrogens (CH). Occasional π-stacking between the aromatic rings, revealed by the scattered peaks between 3 and 4 Å, is present but not prevalent. Again the peak at 10 Å originates from an intra-chain correlation.

Table A2-2: Hole reorganization energies (λ in meV) for TFB oligomers, calculated with the four points method [9] with ORCA 3.0 [10]. In these calculations, octyl and butyl alkyl chains were omitted and replaced with hydrogen atoms. 0 100 200 300 400 500 600

Distribution of transfer integrals

Table A3-2: Calculated mean : < J >, < J t >, ratio of < (J t ) 2 > / < J t > 2 , standard deviation (σ), η, skewness and kurtosis of the distribution of transfer integrals along b axis, for uniaxial strain applied along a, b and c axis : aa , bb and cc , respectively. All values, except η are in meV. η is dimensionless quantity. A3-5: Calculated mean : < J >, < J t >, ratio of < (J t ) 2 > / < J t > 2 , standard deviation (σ), η, skewness and kurtosis of the distribution of transfer integrals along b axis, for uniaxial stress applied along a, b and c axis : aa , bb and cc , respectively. All values, except η are in meV. η is dimensionless quantity.

Relation between the relative variations of drain current and mobility

The drain current measured in transistor operating in linear regime (for V ds < V gs -V th ) can be described by the following equation,

For transistor operating in saturation regime (V d > V gs -V th ) we have instead

At fixed source-drain and gate-source voltages, and if the threshold voltage does not change with the applied stress, the first-order relative variation of drain current in both regimes can be written as,

where I is the drain current, W and L are respectively the width and the length of the transistor channel, C the capacitance and µ the charge mobility. For both configurations (electrodes placed in the cantilever length direction or in the cantilever width direction), variations of W and L can be neglected. In fact, in the former case, the variation of L will be zero and the variations of W small with respect to measured A:35 variations of mobility (the variation of W is indeed the strain, which is of the order of 0.2% while mobility changes are of the order of 10%). For the second case (electrodes in the cantilever width direction), ∆W will be zero and ∆L small (corresponding to the applied strain). We are then left with

Analogously, the relative variation of the capacitance can be written as,

with β and d, respectively, the permittivity and thickness of the dielectric. Since permittivity and thickness do not change upon strain, and geometrical variations ∆W and ∆L can be neglected as explained before, we end up with

Finally, combining Eqs. A3-4-A3-6, it emerges how the main variation of drift current upon application of strain, reflects a variation of mobility,

A:36

: STRAIN-MOBILITY TRENDS IN CRYSTALLINE ORGANIC SEMICONDUCTORS

Relative variation of mobility as a function of transfer integrals is obtained using [(J i( ) /J i ) -1] × 100 (reported as percentage variation), where J i( ) is the transfer integrals obtained when mechanical strain is applied along crystallographic axis i (a, b, c) and J i is the corresponding transfer integral when no strain is applied. Similarly, relative variation of mobility as a function of effective mass is obtained using [(m * i /m * i( ) ) -1] × 100. Similar to the main chapter (Chapter 7), relative variation with respect to transfer integrals would be, henceforth, referred to as ∆µ J and that with respect to effective mass would be referred to ∆µ m A:37 Intra-chain hydrogen bonds computed for (N H) 10 decamer between amide nitrogen atoms and the top (0 = S = 0 group close to phenyl unit) and bottom (0 = S = 0 group close to CF 3 unit) neighboring oxygen atoms are reported in Fig. A5-1. The number of hydrogen bonds with the top oxygen atoms are marginally higher than the bottom ones. This is expected as the top 0 = S = 0 unit is close to the electron donating phenyl group, whereas the bottom 0 = S = 0 unit is close to the electron withdrawing CF 3 group, leading to the variation of partial charge on oxygen atoms.

Hydrogen bonds of amide nitrogen and oxygen atoms of (N H) 10 decamer with surrounding water molecules (henceforth, represented as ODW) are reported in Fig. A5-2 along with the all the possible hydrogen bonds of the decamer with water molecules, i.e., inclusive of the hydrogen bonding of water molecules with all polar groups in the decamer.

A:43 Hydrogen bonds (Dynamic) Hydrogen bonds (Dynamic) (N H -N ) 5 decamer

Intra-chain hydrogen bonds computed for (N H -N ) 5 decamer between amide nitrogen atoms and the top and bottom neighboring oxygen atoms are reported in Fig. A5-3. Similar to the previous case, the number of hydrogen bonds with the top oxygen atoms are marginally higher than the bottom ones.

Hydrogen bonds of amide nitrogen and oxygen atoms of (N H -N ) 5 decamer with surrounding water molecules are reported in Fig. A5-4 along with the all the possible hydrogen bonds of the decamer with water molecules. Compared to amide nitrogen atoms, the number of hydrogen bonds between oxygen atoms of the decamer and the surrounding water molecules is much higher. Further, it can be observed that these hydrogen bonding are the dominant type with a major contribution to all the possible A:44 Hydrogen bonds (Dynamic) Hydrogen bonds (Dynamic) A:45

N 10 decamer

Hydrogen bonds of de-protonated amide nitrogen and oxygen atoms of N 10 decamer with surrounding water molecules are reported in Fig. A5-5 along with the all the possible hydrogen bonds (from polar groups) of the decamer with water molecules. Compared to amide nitrogen atoms, the number of hydrogen bonds between oxygen atoms of the decamer and the surrounding water molecules is much higher, similar to (N H -N ) 5 decamer. Hydrogen bonds (Dynamic) From the hydrogen bond analysis, it evolves that, when de-protonated amide nitrogen atoms are present in the decamers (N 10 and (N H -N ) 10 ), hydrogen bonds between oxygen atoms of the decamer with surrounding water molecules are higher compared to those formed by amide nitrogens, be it intra-chain or with surrounding water molecules. However, when only protonated amide nitrogens are present in the decamer, intra-chain hydrogen bonds between amide nitrogen and the neighboring oxygen atoms are higher when compared to the hydrogen bondings of the polar groups with water molecules.

Conformational free energy

Conformational free energy of interaction is computed using the Adaptive Biasing Force (ABF) method and MD simulations for (N H) 3 and N 3 trimers. Simulation conditions are same as that reported in Chapter 8.

Free energy, reported in Fig. A5-6, is computed as a function of distance between interacting nitrogen atom of the central STFSI units and the top (green curve) sulfur atoms and the bottom (cyan curve) sulfur atoms of the neighboring STFSI units, A:46 independently, with these atoms defining the collective variables to obtain the free energy of interaction. Conformation free energy for (N H) 3 trimers with the interaction distance of ≈ 4 Å between the amide nitrogen and the top sulfur atoms indicate that the π stacked like or eclipsed conformation of neighboring phenyl units is the low energy conformation. The free energy distance of ≈ 8 Å between the amide nitrogen and the bottom sulfur atoms indicate that the CF 3 -SO 2 group is slightly drifted with respect to the P henyl -SO 2 group, thereby reducing the number of hydrogen bonds between the protonated amide nitrogen and the bottom oxygen atoms of the neighboring STFSI unit. Further, the conformation free energy of N 3 trimers is not surprising and indicate that the neighboring STFSI units are drifted apart and no π stacked like conformation is possible. Energy (kcal/mol) Energy (kcal/mol) Conformation free energy between the amide nitrogen and the top sulfur atoms is represented by green curve, whereas that with the bottom sulfur atoms is represented by cyan curve.

A:47

APPENDIX B

: MOLECULAR ORBITAL THEORY Charge carrier occupying any given molecule occupy a specific molecular level, whose shape, extension and energy will effect the charge transport of the charge carrier, substantially. Further, the energies and the wavefunctions of the molecular orbitals also influences the molecular conformation in the bulk. Quantum mechanical description takes into account all these effects. The Schrödinger equation can provide all the information, however, obtaining an explicit analytical solution of the Schrödinger equation is almost impossible in majority of the cases and the development of several approximations maintaining the validity of the equations, becomes important. Advancement of modern computers led to the popularity and rapid development of Computational Methods of Quantum Chemistry, which allows the quantum treatment of either a single molecule or groups of molecules, with a full quantum description possible. A brief review of popular quantum chemistry methods is provided in this chapter, following as closely as possible the approach of Jensen [1] Born-Oppenheimer Approximation

The Schrödinger equation for a molecule can be written, without considering the magnetic interactions, as

where, r i and R i takes the set of coordinates of the electrons r i and nuclei R i , and.

B:1

The kinetic energy of electrons, where i runs over all the electrons

The kinetic energy of the nuclei, when α runs over all the nuclei.

The electron nucleus Coulombic attraction, where r i and R α denote the position vectors of electrons and nuclei, respectively and Z, e represent the atomic number and elementary charge, respectively.

Electronic Coulombic repulsion between electrons i and j.

Nuclear Coulombic repulsion between nuclei α and β.

The Born-Oppenheimer approximation relies on the mass disparity (and thus the velocity) of nuclei and electrons, which leads to an imbalance in time between the motion of nuclei (slow process) and electrons (fast process), such that electrons can settle in molecular orbitals (MO) at a given nuclear configuration. This entails that, the electrons adiabatically follow the nuclei in a succession of states of equilibrium ψ e ( r, R), solutions of stationary Schrödinger equation with the fixed nuclei. Therefore, the electronic wavefunctions are considered to depend on nuclei positions, parametrically. The wavefunction of the system can then be factorized as

where φ n depends only on the nuclear positions and ψ e describes the electronic states in a fixed nuclei configuration. The electronic Schrödinger equation can then be written as

with the R coordinates constant/ H e ( R) = T e + V en + V ee is the Hamiltonian of the electronic system and E e ( R) is the eigenvalue of the wavefunction ψ e ( r, R). Substituting Eqs. B-2 and B-3 in the general equation (Eq. B-1) and having the operator T n acting only on the nuclear wavefunction, results in, B:2

where,

is the Hamiltonian of the nuclear system, with the eigenvalues E M and the eigenfunction φ n ( R). The nuclei are subject to the effective potential

, sum of the repulsive potential V nn and the attractive one given by the electronic contribution E e ( R).

The original Schrödinger equation for the molecule (Eq. B-1) has been decoupled into two equations, one for electrons and the other for nuclei. To achieve this result, the operator T n can act only on the nuclear wavefunction φ n ( R), and in reality the action of T n on Eq. B-3 is given by,

so that, in summary Born-Oppenheimer approximation consists in neglecting the last two terms in Eq. B-5. For many systems of interest, this is a reasonable approximation. For some other systems, this is not the case as the transitions between two Born-Oppenheimer surfaces happen through the coupling of vibrational levels.

Symmetrization and Anit-Symmetrization Operators

If the state of the single particle u a is a ket in the Hilbert space H, the system composed of N identical particles live in the Cartesian product H N . A simple composition of such one-particle state can be denoted as,

A transposition of any two particle i and j (with i > j) is performed by the operator T ij given by,

Eq. B-7 is hermitian and unitary as can be easily verified. The operator P σ performing a permutation σ of the set of identical particles, B:3

Thus, a permutation P σ is unitary, but in general not hermitian because the transpositions {T ij } does not always commute with each other. Symmetrization and anti-symmetrization operators, which generate the states of bosons and fermions respectively, are introduced as follows,

Using the properties from Eqs. B-11 and B-12,

we obtain, Eqs. B-13, B-14

i.e., operators S and A are idempotent. As, S and A are unitary and idempotent, it follows that they are projection operators; they identify H S and H A and the two subspaces of H N . The subspace H S contains all the symmetric state of H N , while subspace H A contains all the anti-symmetric ones. Further, given that,

these subspaces are also orthogonal, i.e, H S ⊥ H A B:4

Hartree-Fock

The electronic time dependent Schrödinger equation obtain invoking the Born-Oppenheimer approximation can be exactly solved only in a few, very simple, cases. Employing variational principle approximate solutions for the ground state can be progressively refined, such that, any approximate wavefunction has an energy above or equal to the exact energy; and this energy is equal only if the approximate function is the exact wavefunction. A trial wavefunction ψ γ containing certain number of parameters γ can be constructed and the best wavefunction of this form ψ * γ can be obtained my minimizing the energy, as a function of the parameters γ,

Hartree-Fock theory employs the ansatz that the electronic wavefunction has the form of a single Slater determinant,

The electronic wavefunction can be written in terms of zero, one and two electron terms,

where h 0 is the additive constant given by the repulsion of the nuclei at a particular position, and h i is the electron-nuclei attraction and g ij is the describes the interaction between two electrons, given by,

Using the concepts from B, the energy for Eq. B-17 can be obtained as,

where ĥ and ĝ are the one electrons and two electron operators, given by, ĥ

Introducing, the Coulomb operator (J i ) and exchange operator (K i ), defined as,

and taking the total and Coulomb and exchange operators, we have,

Eq. B-21, can be rewritten as,

It is now, to define a set of MOs that makes the energy a minimum, while the variation of such orbitals is such that the MOs remain orthogonal and normalized. This is a constrained optimization and can be solved by means of Lagrange multipliers, i.e. finding the stationary point of the Lagrangian functions,

Using Eq. B-28, the variation of the energy with respect to the functional variation of the MOs can be written as, B:6

where, F is the Fork operator, given by,

The operator F can be considered as an effective one-electron energy operator, including the attractive potential of the nuclei as well as the average of the repulsive Coulombic and exchange interactions of the other electrons (with j = i). In particular the exchange interaction enter into the picture as a purely quantum mechanical effect coming from Slater determinant form of the wavefunction for identical fermions. The variation of the Lagrangian (Eq. B-29) can be written as,

As the variations |δϕ i and |δϕ i * can be chosen independently, Eq. B-32, can be split into two equations, given by,

Taking the complex conjugate of Eq. B-34 and subtracting it from from B-33, we have,

Eq. B-35 entails that the matrix of the Lagrange multipliers λ ij is Hermitian. From Eq. B-33, the Hartree-Fock (HF) equations can be written as,

but, as λ ij is Hermitian, a transformation of the MOs exist such that,

For the numerical solution of the HF equations, the unknown MOs are expressed in terms of set of unknown functionals, the basis sets. The basis set is to be chooses carefully, in a way that functions of the basis are somehow able to capture the underlying physics of the problem. Employing a basis set allows to pass from an infinite dimension problem (in Hilbert space, L 2 ) to a finite dimensional algebraic problem.

The behavior of such functions should be able to capture the physics of the problem, allowing to approximate better the real solution. Every molecular orbital is expanded as a linear combination of atomic orbitals (LCAO).

where {χ a } are forming the basis of dimension M . This leads to the Roothan-Hall (RH) equations, a generalized eigenvalue problem for a M × M matrix, and the forth order tensor, G abcd involving the two electron integrals.

B:8

The solution of RH equations (Eq. B-39) are M eigenvalues i with their associated set of coefficients {C i1 , C i2 , • • •, C iM }, i 1 to M . Each index i represents an MO giving the solution of M molecular orbitals; and if N is the number of electrons, then only N/2 of them are occupied (in a closed-shell system). Once the solution og RH equations is obtained, a new FOrk matrix can be constructed. This process is iterative in nature, repeated in self-consistent field (SCF) calculations till the convergence is achieved. It is important to notice that, the occupied orbitals are variationally optimized, because they do contribute to the total energy and then to the Fork operator at a particular SCF step of the optimization cycle.

In the solution of the HF method, the electrons are sitting on the effective orbitals coming from an average interaction. If the basis set is sufficiently large, the HF solution accounts to ≈ 99% of the total energy, but unfortunately, the remaining ≈ 1% are often very important to describe the chemical phenomena. This missing piece in HF is the electronic correlation. In an attempt to account for the correlation, an ansatz can be employed such that a linear combination of Slater determinants corresponds to the exited state configurations. This gives the origin of the so called Hartree-Fock methods, such as configuration interaction (CI), coupled-cluster (CC) and Møller-Plesset perturbation theory (MP). These approaches allows for the systematic increase of accuracy, but dramatically increasing the computational cost, with the number of exited levels considered.

Semi empirical Methods

The cost of performing an HF calculation scales a the forth power of the number of basis functions (the cost of calculating the two electron integrals needed to calculate the matrix F ab , at every step). The first step in reducing the computational cost is to consider only the valance electrons, explicitly, while accounting the core electrons by reducing the nuclear charge (effective screening). A minimal basis set can be used for the valance electronic states. While for hydrogen atom 1s function is considered, other atoms in the second and third row of the periodic table are described with four functions of s-and p-orbitals. Most semi-empirical methods use only s-and pfunctions, and Slater type orbitals as the basis function [1] The central assumption of the semi-empirical methods is the Zero-Differential Overlap (ZDO) approximation, which neglects all the products of the basis functions located on different atoms. The integrals that are different from zero are considered as parameters and are determined based on the theoretical calculations or experimental data. This is also called as parametrization of a particular semi-empirical method.

B:9

NDDO

Employing ZDO approximation leads to the Neglect of Diatomic Differential Overlap (NDDO) method. The three-and fore-centered integrals vanish and the overlap integral matrix is reduced to

where, µ and ν are two atomic orbitals. The one-electron integrals are given as,

where V A is the atomic potential (nucleus and the core electrons) of atom A and mu A and ν A are the atomic orbitals of the atom A. Due to orthogonality of the atomic orbitals, µ A |∆ 2 /2 + V A |ν A = 0, unless µ = ν/ The two electron integrals are approximated as,

INDO

The Intermediate Neglect of Differential Overlap (INDO) [2,3,4] approximation neglects all the two-centered integrals that are not of the Coulomb type, in addition to those neglected by the NDDO approximations. In order to preserve the rotational invariance, i.e. the total energy should be independent of rotation of the coordinate system, the integrals of the form µ A |V A |ν A and µ A ν B |ĝ|λ A σ B must be taken independent of the orbital types (an integral involving the p orbital must be the same as with the s orbital). This leads to the INDO method, involving the following approximations in addition to those made by NDDO method. One-electron integral then becomes,

Two-electron integrals are given by,

The non-zero two-electron integrals are usually denoted by γ, given by,

ZINDO ZINDO [5] is a parametrization of INDO that covers a wide range of elements in the periodic table, and is widely used to compute transfer integrals for organic electronic materials.

Density Functional Theory

The Density Functional Theory (DFT) allows us to calculate the ground state properties of a system without dealing (at least in principle) with may electron wavefunctions |ψ . In DFT the main quantities of interest is the electron density,

where x i is the representation of all degrees of freedom of particle i, position and spin x i = ( r i , σ i ). In the Schrödinger equation for the electrons in Born-Oppenheimer approximation, the Coulomb potential arising from the nuclei is treated as a static potential, given, by, As F is the same for all N -electron systems, the Hamiltonian and hence the groundstate |ψ 0 , are completely determined by N and Vext . The average energy is then given by,

where the energy of the external potential, as in the Thomas-Fermi theory, can be expressed as,

Density Functional Theory is based on two fundamental theorems by Hohenberg and Kohn [6,7].

Theorem1: The external potential v( r) is uniquely determined, up to a constant, by the ground state electron density n 0 ( r).

As a consequence, all the properties of the system are determined by the electronic ground state density. Indeed, since n 0 ( r) determined both the external potential v( r) and the number of electrons N = d rn 0 ( r), it also characterizes unambiguously, the Hamiltonian H and thus implicitly all the properties that can be derived from H through the solution of time-dependent or time-independent Schrödinger equation, like the wavefunctions |ψ . Further, there exists a functional E[n], that expresses the energy in terms of the electron density, for any particular v( r), given by,

where, F HK is a universal functional, completely independent of the external potential v( r) and it thus te same for any set of N electrons.

Kohn-Sham Scheme

The idea is to replace the system of interacting particles in the external potential v( r) with a system of non-interacting particles in some other external potential v s ( r), such that, (-h2 2m ∆ 2 + v s ( r))ϕ( r) = i ϕ i ( r) (B-58)

These are Khon-Sham equations [8], i is the orbital energy of the corresponding Kohn-Sham orbital ϕ i , and the density of the N particle system is given by,

The effective potential v s ( r) will be choose such that the ground state density of the auxiliary system n s ( r) is equal to the original ground state density n( r). The key point of Kohn-Sham theory is to calculate the kinetic energy under the assumption of non-interacting electrons,

where the notation implies that T depends on ϕ i which in turn depend on the ground state density n. In reality, the electrons are interacting and Eq. B-60 does not provide the total kinetic energy, but usually, the most of it. The remaining kinetic energy is absorbed into an exchange-correlation term, and in general the expression can be written as, B:13

Equating Eq. B-61 with Eq. B-56, E xc can be obtained as,

i.e., the kinetic correlation energy plus the potential correlation and exchange energy.

The exchange part can also be written following the HF theory as,

leaving the main problem to treat correlation. For the second theorem, the solution is the minimum of the functional (Eq. B-61) but as T s is written as the dependent form of effective orbitals, a direct minimization is not possible. In the Kohn-Sham scheme we can write for interacting systems, can only be calculated explicitly once an approximate E xc [n] has been chosen and since it depends on n( r), which is unknown, the solution of Kohn-Sham equations has to be found self-consistently.

The improvements brought by the Kohn-Sham theory come at the price of reintroducing the molecular orbitals (3N variables), analogous to the Hartree-Fock theory, but in a method that is still quite simple, for being capable to capture many-particle correlations (once the correct xc has been chosen).

APPENDIX C

: LIST OF PUBLICATIONS Publications related to this thesis