Abstract
This paper describes 3D particle-in-cell simulation of charge injection and transport through nanocomposite film comprised of ferroelectric ceramic nanofillers in an amorphous polymer matrix and/or semicrystalline ferroelectric polymer with varying degrees of crystallinity. The classical electrical double layer model for a monopolar core is extended to represent the nanofiller/nanocrystallite by replacing it with a dipolar core. Charge injection at the electrodes assumes metal-polymer Schottky emission at low to moderate fields and Fowler-Nordheim tunneling at high fields. Injected particles propagate via field-dependent Poole-Frenkel mobility. The simulation algorithm uses a boundary integral equation method for solution of the Poisson equation coupled with a second-order predictor-corrector scheme for robust time integration of the equations of motion. The stability criterion of the explicit algorithm conforms to the Courant-Friedrichs-Levy limit assuring robust and rapid convergence. Simulation results for BaTiO3 nanofiller in amorphous polymer matrix and semicrystalline PVDF with varying degrees of crystallinity indicate that charge transport behavior depends on nanoparticle polarization with antiparallel orientation showing the highest conduction and therefore the lowest level of charge trapping in the interaction zone. Charge attachment to nanofillers and nanocrystallites increases with vol% loading or degree of crystallinity and saturates at 30–40 vol% for the set of simulation parameters.
1. Introduction
Nanocomposite films comprising ceramic nanofillers in amorphous polymer matrices are of interest in capacitive energy storage for rapid power cycling applications. Advantages include high energy density, low cost, lightweight, and small footprint. Nanocomposite films combine the processability and high breakdown field strength of the polymer with the high dielectric constant of the fillers, engineered to an acceptable level of dielectric loss. Most of the increase in the effective dielectric constant comes from an increase in the average field in the polymer matrix with very little of the energy being stored in the high permittivity phase. Large contrast in permittivity between the two phases gives rise to highly inhomogeneous electric fields () in the “interaction zone,” defined as the interfacial region that surrounds the nanofillers and interspaces. Within this region, mobility is field-dependent and may also depend on the physicochemical bonds on the common surfaces. Morphology plays an important role in film layers with very high surface-to-volume aspect ratios. Highly inhomogeneous fields and structural inhomogeneity generally lead to a significant reduction in the effective breakdown field strength of the composite, limiting the increase in the energy storage capacity and energy density. So a careful compromise needs to be engineered at the polymer-ceramic interfaces to attain a viable solution [1]. Interface engineering attempts to improve the dispersion of filler in polymer matrix. In a well-dispersed composite, interfaces between the ceramic nanoparticle and polymer phases create effective electron scattering and transport centers, thus reducing the breakdown probability. Moreover, well-dispersed ceramic nanoparticles block degradation tree growth and thus increase the long-term breakdown strength.
The polymers currently used as matrices in the dielectric nanocomposites, including polyethylenes (PE), poly(methyl methacrylate)s (PMMA), epoxy resins, and polyimides (PI), usually possess dielectric permittivities of ~2–5. significantly lower than their inorganic counterparts. Poly(vinylidene fluoride) (PVDF) based ferroelectric polymers exhibit large spontaneous polarization and high dielectric constants (~10 at 1 kHz) because of the presence of highly electronegative fluorine on the polymer chains and the spontaneous alignment of C–F dipoles in the crystalline phases. PVDF copolymers, such as poly(vinylidene fluoride-co-trifluoroethylene), P(VDF-TrFE), and poly(vinylidene fluoride-co-hexafluoropropylene), P(VDF-HFP), have also been utilized in the formation of the dielectric nanocomposites. However, dispersion of inorganic fillers in these fluorinated polymers is always problematic because of the low surface energy of the polymers. The agglomeration of the ceramic dopants gives rise to electron conduction for a high dielectric loss and undesirable porosity for dielectric failure at much lower fields. PVDF is a semicrystalline polymer with varying degrees of crystallinity. The ferroelectric properties of crystalline PVDF stem from the net molecular dipole moment, oriented in the plane of the carbon backbone and along the polymer chain [2].
Considerable progress has been made over the past several years in the enhancement of the energy densities of the polymer nanocomposites by tuning the chemical structures of ceramic fillers and polymer matrices and engineering the polymer-ceramic interfaces. Ferroelectric ceramics, such as BaTiO3 (BT), Pb(Mg1/3Nb2/3)O3-PbTiO3 (PMN-PT), or other ferroelectrics or relaxor ferroelectrics, possess a very large dielectric permittivity. Large enhancements in the electric energy density and the electric displacement in the nanocomposites of P(VDF-TrFE-CFE) terpolymer/ZrO2 nanoparticles have been demonstrated through the interface effect with the presence of 1.6 vol% of ZrO2 nanoparticles [3]. Surface modifications include grafting, functionalization, and bonding. Interface engineering attempts to improve the dispersion of filler in polymer matrices include the grafting of polystyrene (PS) directly onto TiO2 surfaces via atom transfer radical polymerization (ATRP) to result in a dielectric constant enhancement of over three times that of bulk polystyrene. The high permittivity core lends high capacitance, while the flexible polystyrene shell endows the materials with good dispersability and film-forming properties. A second method is the functionalization or Phosphonic acid surface-modification of BT nanoparticles to yield well-dispersed P(VDF-HFP) composites [4]. Chemical functionalization of the BT nanoparticles with ethylene diamine moieties on the surface also rendered the particles as homogeneous distributions in the poly(vinylidene fluoride-co-chlorotrifluoroethylene) [P(VDF-CTFE)] and poly(vinylidene fluoride-ter-trifluoroethylene-ter-chlorotrifluoroethylene) [P(VDF-TrFE-CTFE)] matrices [5]. Another method is chain-end functionalization of the ferroelectric polymer and subsequent covalent-bonding of the polymer nanocomposites to form uniform nanoparticle dispersion. Ferroelectric polymers have also been terminated with Phosphonic acids together with subsequent utilization of the reactive end-groups of the polymer for directly coupling with oxide nanoparticles to yield the covalent-bonded nanocomposites. It was demonstrated that covalent assembly of polymer matrix and inorganic filler not only led to highly dispersed nanofillers without additional surface modification but also provided great stability and offered enhanced interfacial interactions for high polarization responses under the applied fields [6]. High-energy-density polymer nanocomposites based on surface-functionalized TiO2 nanocrystals used as dopants in a ferroelectric P(VDF-TrFE-CTFE) have been shown to possess comparable dielectric permittivities of 42 and 47, respectively, measured at room temperature and 1 kHz. High dielectric performance in the nanocomposites is realized via the large enhancement in polarization response at high electric fields and changes in polymer microstructure induced by the nanofillers. The energy density of the P(VDF-TrFE-CTFE) nanocomposites on the TiO2 concentration was reported to peak at around 10 vol% TiO2 content. This feature is likely associated with the interface effect that is proportional to the interfacial area. The thickness of the interface region is generally estimated to be ~20 nm when the volume fraction of the polymer chains residing in the interfacial area reaches a maximum at 10 vol%. The large interface area in the nanocomposites would produce the Maxwell-Wagner-Sillars (MWS) interfacial polarization at low frequencies contributing to the creation of the interaction zone with the Gouy-Chapman diffuse layer, thereby greatly affecting polarization and the dielectric response of the polymer matrix. The incorporation of the TiO2 nanoparticles into the polymer induces an improved electric displacement, which accounts for high energy densities observed in the nanocomposites [7]. Finally, surface hydroxylation treatment using hydrogen peroxide to modify the surface of BT nanofillers dispersed in a ferroelectric copolymer host has resulted in up to two orders of magnitude reduction in the leakage current of nanocomposite thin-film capacitors. This reduction is observed concurrently with the enhancement of the effective permittivity and breakdown strength of the thin-film nanocomposites [8].
Characterization studies on PVDF have involved dynamic pyroelectrical measurements and thermally stimulated depolarization (TSD) experiments, from which the permanent polarization of the dipoles in the crystalline phase, the frozen-in polarization of the dipoles in the amorphous phase, and the effect of the excess charges of the MWS interface polarization can be identified [9]. The influence of uniaxial stretching on the dielectric properties and remanent polarization of α-PVDF into α-phase oriented films were investigated by differential scanning calorimetry (DSC) and infrared spectroscopy (FTIR) and verified the increase in crystallinity by the reduction of the amorphous-crystalline interface, stable remanent polarization, and the permittivity at 1 kHz [10]. Experimental studies on the dependence of polarization fatigue on crystallinity of P(VDF-TrFE) copolymer films indicated that high crystallinity of the film correlated with slower fatigue rate of the film. A possible explanation was the space charge trapped at the boundaries of crystallites and/or captured by defects in both the amorphous and crystalline phases as the major contributor to polarization fatigue [11].
In matching nanofillers with the polymer matrix, the challenge is to understand the role of the interaction zone where the very large area to volume ratio of the interfaces in nanocomposites would have significant impact on the electrical and dielectric properties of the polymers.
2. Modeling Rationale
2.1. Nanocomposite Models
The physics of the interfaces is still not well understood. Key questions remain on how the interfaces are related to the structures of the components and how the interfacial properties change with the filler size, orientation, loading, and the applied electric field. Detailed modeling guided by theoretical insight and in-depth dielectric characterization is necessary to derive understanding before rational designs can be used to guide the synthesis of hybrid structures with optimally engineered interfaces. Creation of such a model begins with an appropriate representation for the nanofiller/nanocrystallite and its charge transport behavior.
Several models of nanoparticles are discussed in the literature [12]. The Tanaka Multicore 3-layer model consists of (1) a first layer on the nanoparticle surface only a few nm thick forming a region of chemical bonding that is ionic, covalent, hydrogen, or van der Waals’ force; (2) a second layer corresponding to a region of some ordered polymer chains, which is some 10 nm in thickness; and (3) a third layer which is dominated by the far-distance force originating from the electric double layer (EDL) which is several 100 nm thick, corresponding to the Debye shielding length. Another is the Lewis 2-layer model comprising an inner bound Stern layer and a diffuse Gouy-Chapman (Debye-Hückel) outer layer created around each nanoparticle to determine the dielectric properties. The classical EDL is similar to the Lewis model. It assumes a monopole net charge for the core particle with a Stern-Helmholtz first layer of counter-charge tightly attached to the particle surface forming a double-layer charge structure. The second Gouy-Chapman layer of diffuse charge is loosely attached by Coulomb force, electrically screening the first layer. Charge transport is enabled by the increase in nanoparticle loading within the dielectric composite eventually leading to overlap of the diffuse double layers forming a conducting path and thereby greatly affecting polarization and dielectric response of the polymer matrix. Bulk charge accumulation is reduced due to this internal conductivity, and the dielectric breakdown strength of the nanocomposite is improved.
In this paper, we propose to extend the EDL to model ferroelectric nanofillers and nanocrystallites by substituting a dipole for the core. This eEDL model would allow initial charge attachment on the opposite-signed end of the dipole to form the Stern-Helmholtz layer. Subsequent arriving charge may be repelled by the initial attached charge forming the diffuse transport layer. Charge transport will occur by meandering of the trajectories due to Coulomb repulsion and use of a field-dependent mobility, . The MWS effect due to charge build-up could lead to polarization inversion. It is well known that the polarization and dielectric properties of the PVDF based ferroelectric polymers originate from their crystalline domains. We propose to simulate semicrystalline PVDF as nanocrystallites in an amorphous polymer matrix where the ferroelectric nanoregions are treated with the eEDL model. Ceramic nanofillers such as the Perovskite BT in an amorphous polymer matrix are also represented with the eEDL model. The charge particle transport simulation will determine electrode field modification and allow the passage of charge particles, tracking of trajectories, and calculation of charge fractions.
2.2. Current Simulation Efforts
Current simulation efforts on nanocomposites are based on effective medium theory (EMT) for charge percolation to calculate complex permittivity using Maxwell-Garnett, Bruggeman, and Loyenga models [13]. For continuum calculations of energy density, effective permittivity may be estimated by any of the preceding methods, including the Lichtenecker logarithmic rule. Another uses Matlab to create phenomenological 2-layer models to estimate relative permittivity incorporating water uptake and matrix reorganization although actual model details are lacking [14]. Percolation methods [15, 16] are applied when conductive or semiconductive fillers are used in the insulator matrix to improve effective dielectric constant, . As the volume fraction of the fillers increases to the vicinity of the percolation threshold , of the composites can be described by the power law: where is the dielectric constant of the polymer matrix and is an exponent of near unity. Using this strategy, the dielectric constant of composites can be increased one to two orders of magnitude over the polymer matrix. The end result is the qualitative prediction of possible conductive paths through the nanocomposite medium. This picture of the electrical behavior of particulate nanocomposite materials is attributed to the formation of a connected network of conductive inclusions within the insulating polymer matrix. The implicit assumption is that charge is able to transfer between localized states through hopping conduction and/or electron tunneling.
Bipolar charge transport models capable of handling leakage current up to prebreakdown levels have been successfully applied to layered polymer films [17–19]. Continuum charge transport models are not suited to simulate material with morphology, especially at the nanometer length scale. Probing the role of interfaces in nanocomposites for charge transport requires quantitative details of charge migration through the interspaces, a requirement beyond the capability of continuum models. Particle simulation or specifically a class of “particle-in-cell” (PIC) methods is needed to represent the discrete charge entities and track their migration as particles through a nanocomposite rheology consisting of random nanofillers and/or nanocrystallites in an amorphous matrix. In order to study the role of interfaces, it is necessary to facilitate charge particle transport through the interspaces with a field-dependent mobility under the composite field while allowing for particle-particle and particle-nanofiller/nanocrystallite interactions. Nanoregions with ferroelectric properties may be modeled as spherical dipoles with remanent dipole moment and polarization orientation that are random, in-plane, parallel, and antiparallel. This paper will study the effects of nanofiller loading, degree of crystallinity, and polarization orientation on charge conduction, field modification, and distribution of charge fractions.
2.3. The Proposed eEDL Model
The eEDL model for ferroelectric nanofillers and nanocrystallites is shown in Figure 1 where the core is replaced by a dipole of the appropriate physical diameter and dipole moment. In Figure 1(a), incoming positive charge is repelled by the positive end and attracted towards the negative end of the dipole. Charge is allowed to attach on impact forming the bound Stern-Helmholtz layer. Subsequent incoming charge may be repelled to form the diffuse outer Gouy-Chapman transport layer. The cumulative charge build up on opposing ends of the dipole leads to MWS polarization as shown in Figure 1(b). The gradual charge deposition and formation of the diffuse layers as charge percolates through the nanocomposite film create the interaction zone which consists of the interspaces formed from overlapping interfacial regions. Trajectories for charge that make it through the film to the counter-electrode are curvilinear paths that meander through the interspaces, as shown in Figure 1(c). Charges that arrive at the counter-electrode contribute to the conduction of the film. These charges are considered as neutralized and therefore do not contribute to the field.
(a)
(b)
(c)
3. Problem Formulation
3.1. Electrostatic Fields
The composite field due to the bias, space charge, and polarization charge is given by
The electrostatic field, , is the gradient of a scalar potential, , resulting in the Poisson equation for charge conservation where is the space charge:
The potential for a point charge, , is
When the particle count, , is small to moderate, fields from the space charge are computed from superposing or summing the point source solutions for all particles considered:with as the Green function fundamental solution to a point charge in the volume . The polarization, , from sphere-shaped ferroelectric materials may be written aswhere is the ratio of the permittivity of the sphere to the polymer matrix, is the Clausius-Mossotti factor, and is the susceptibility. The dipole moment, potential, and field are given bywhere is the intrinsic dipole moment, is the distance of the evaluation, and is the nanocrystallite radius. The total field is obtained by summing both dipole contributions in (5).
3.2. Solution of Poisson Equation with the BIEM
The solution of the Poisson equation is obtained using an integral equation method derived from the integral form of the divergence theorem. Volume integrals for homogeneous space may be reduced to surface integrals thus reducing the dimensionality of the problem region. The boundary integral equation method (BIEM) field solver is based on the source distribution technique (SDT) where mobile space charge and bulk trapped charge are treated as volume sources and conducting boundaries and material interfaces are replaced with appropriate distributions of free and bound interfacial polarization and trapped charges to satisfy the specified boundary and interface conditions. Once the charge distributions are ascertained, field parameters are evaluated by superposition of the integral contributions from all sources.
The solution to the Laplacian given by the following equation:is a Fredholm integral of the first kind for potential:where is a charge distribution on surface and , the Green function, is given bywhich naturally satisfies the vanishing far-field condition. The field is related to the potential byand it can be computed by integration of analytically differentiable kernels where, for example, the flux which is the normal field on the surface is given bywhere the coefficient in the second term on the right is the Cauchy principal value (CPV) integral with a reentrant angle of . Enforcement of boundary conditions for potential and flux and interface conditions at material interfaces for continuity of tangential and normal ():for the Poisson equation results in for Dirichlet, Neumann, and dielectric interface conditions, respectively, where and are the potential and normal derivative on the surface. Kernel functions, including the Green function and its derivatives (, , and ), are integrated numerically using Gauss-Legendre quadrature. Singular kernels are accurately computed using minimum-order sampling by tying the quadrature weight function to the singularity. Details on the computation of these integrals are discussed in detail for axisymmetric and 2D geometries [19, 20]. Discretized forms of these three equations are solved simultaneously for on electrodes and on interfaces. Equations (14) invoke superposition of sources that include free charge, , on electrodes (); polarization charge, , and trapped charge, , on material and physical interfaces ; and mobile and trapped space charge, in the volume . With as electrode area and as the current density, finite numbers of positive and negative charge particles, , are injected at each time step, , from the cathode and the anode into the respective adjacent polymer according towhere the charges are injected from randomized locations on the surface of the corresponding electrodes.
Free charge at the two electrode-polymer interfaces and polarization and trapped charge at the polymer-polymer interface are related to the electric fields on sides 1 and 2 bytherefore allowing direct inference of the fields.
3.3. Matrix Algebra for Rapid Computation and Reuse
Elements of the discretized matrices on the left hand sides of (14) are calculated from spatial integrals and therefore need not be recomputed for every time step. Large portions of the matrices may be saved and reused in matrix algebra for rapid computation of the sources or independent variables, and . The forcing functions on the right hand side may be approximated as summations of contributions from point charges as shown in (5). More details are available in the literature [20].
For large numbers of particles, , and spatially varying distributions a “scatter-gather” method may be used to scatter the space charge, , to on the vertices of the inscribed volume mesh using, for example, trilinear interpolation: where are the trilinear interpolation functions and is the mesh volume. This particle-mesh scheme facilitates integration by gathering the contributions from the vertices, , of all volume meshes, :to compute the field at the position. This method requires only () calculations for a mesh of size . Following (17), the current density, , can be similarly expressed aswith () as the drift velocity. Contributions of these volume sources to the field are evaluated as the algebraic summations of the field from each volume mesh as shown in (18).
3.4. Equations of Motion
The electrodynamic simulation of space charge migration through the nanocomposite film involves the time-dependent integration of the two equations of motion:where and are the mass and the charge of the particle, is its velocity, and is the electric field at the particle location. The equations above are followed in time as they evolve in velocity and position space (phase space) to determine the trajectory of each particle. The driving terms in (20) are the Coulomb and dipole forces. Due to the small acceleration time and short mean free path, the charge is assumed to migrate with a drift velocity, , where is the mobility in the polymer.
4. Simulation Algorithm
4.1. Particle Simulation
Particle simulation is achieved within a computational cell with the solution of the Poisson equation for charge conservation and integration of the equations of motion performed at each time interval as a two-step field solve-particle push algorithm [21]. Nanofillers and nanocrystallites are randomly distributed within the cell to the prescribed vol% loading using “hard sphere” logic, that is, allowing contact but no overlap. High loading will require some ordering, like hexagonal close-packed (HCP) to achieve a volume fraction of 0.9069. The 3D PIC model simulates the dynamics of charge particle injection, particle-particle interactions, and particle attachment to nanofillers and nanocrystallites during migration through the nanocomposite layers between the top and bottom electrodes. In accordance with the proposed eEDL model, particles initially attach to nanofillers and nanocrystallites on contact to form the bound Stern-Helmholtz layer on the opposite charged end of the dipole. Subsequent particles form the diffuse Gouy-Chapman transport layer due to Coulomb and dipole force repulsion. The self-consistent field solution is obtained using the BIEM with Dirichlet boundary conditions for the electrodes, homogeneous Neumann boundary conditions for the four vertical side-walls, and includes contributions from the mobile, attached, and polarization charge. Charge particles are injected from the electrodes over an field range using Schottky emission at low-to-moderate fields and Fowler-Nordheim tunneling at high fields. Charge particles migrate with Poole-Frenkel field-dependent mobility subject to velocity saturation at high field. Recombination is based on Monte Carlo statistics. Particle tracking involves a predictor-corrector algorithm to integrate the equations of motion subject to Coulomb, dipole, and bias field forces. A self-consistent particle-particle, particle-mesh (P3M) scheme may be implemented for very large numbers of particles [22–24].
4.2. Computational Cell
A computational cell is used to cater to the nanoscale simulations of discrete, randomly distributed nanofiller and nanocrystallite particles in an amorphous polymer matrix. The computational cell is a cuboid of edge dimension 1 μm, volume of 1 μm3 bounded by 4 vertical side walls with zero flux condition (), and top (anode) and bottom (cathode) electrodes at constant potential to maintain the bias field. Discrete numbers of charge particles, commensurating with the current density and time step, are continuously injected from the electrodes as a function of the averaged electrode field. Injection locations are randomly dispersed over the electrode surface. Charge particles migrate through the nanofiller and nanocrystallite distribution under the composite field, initially attaching on contact to form the bound layer. Subsequent charge particles are repelled by the bound layer and meander through the interspaces formed by the interaction zone, outlining curvilinear trajectories. Charge particles that arrive at the counter-electrode are considered to be neutralized and therefore have no further contribution to the composite field.
To minimize local error, displacement, , is required to be smaller than the Debye length; that is, , where is the largest charge number density. The time step, , needs to be shorter than the dielectric relaxation time, , characteristic of charge fluctuations to decay. The stability criterion of the explicit algorithm is given by the Courant-Friedrichs-Levy (CFL) limit, :which represents the ratio of mobile charge velocity, , to trace velocity, . In accordance with (22), for stability, the trace velocity cannot be faster than the speed of the charge [25].
Simulation proceeds with the following steps where setup consists of tasks (1) and (2) and time-dependent solution involves iteration of tasks (3) and (4) as follows.(1)Generate geometric model of computation cell with surface meshes for numerical discretization and integration.(2)Generate random distribution of nanoregions and assigned polarization orientation to given vol% loading.(3)Solve for surface sources on electrodes and interfaces.(4)Integrate surface sources to compute field on charge particles for use in dynamic tracking through nanoregions using the PIC algorithm.
4.3. Time Integration Algorithm
The two differential equations given by (20) and (21) are integrated numerically using a predictor-corrector method of order or second-order accurate in time. The dependent variables, and , are iteratively improved in two steps at each time level. To start, we assume that, at time level , the sets () and () are known. Then, the predictor-corrector equations are as follows. Predictor: Corrector: Update the following:
In the predictor stage, the initial values for the () iteration level are improved. These values are then used in the corrector stage to obtain the updated values for the next time level. The time-dependent algorithm loops through the following steps.(1)Inject a number of particles according to Schottky or Fowler-Nordheim injection.(2)For each particle, we have the following. Predictor Compute field at . Compute and . Compute new and . Corrector Compute field at . Compute new , . Particle push is obtained.(3)Update iteration arrays.(4)Check particle location: (a) remove if deposited on counter-electrode or (b) reintroduce from opposite side-wall if exited from side-wall.(5)If no mobile charge left in polymer, exit. Else return to step (1).
The preceding tacitly assumes that the field is “frozen-in” relative to the time step of charge migration and the transient time to attain terminal velocity is much shorter than the mean free path or the time between collision events between oppositely charged particles.
5. Charge Injection and Transport
5.1. Charge Injection
Charge injection from a metal like Al into vacuum requires a work function of ~4.1 eV. For injection into a polymer such as LDPE, the work function is ~1.2 eV. Charge injection from the electrode into the polymer is further facilitated by lowering of the potential barrier through a combination of the image force and the applied field. Assuming Boltzmann statistics for charge carriers in the organic semiconductor, electrons need to migrate from the highest occupied molecular orbital (HOMO) to the lowest unoccupied molecular orbital (LUMO) through a barrier defined as the difference between the Fermi level in the injecting electrode and the center of the density of states (DOS) distribution at the LUMO-level. On the other hand, holes are injected to the HOMO through a barrier defined as the difference between the Fermi level in the injecting electrode and the center of the DOS distribution at the HOMO-level.
At moderate applied fields, electrons or holes within the electrode arrive at the metal-polymer interface with enough energy to leave the metal surface. Charge injection from a metal electrode into the LUMO band of the polymer by Schottky barrier thermionic emission are given bywhere is the Richardson constant ( A/m2·K2) and and are the energy barriers to injection in . The combined effect of the image force and the applied field results in a lowering of the barrier potential given by
At higher applied fields, the slope gets steeper and the barrier is further lowered so that the tunneling length is much shorter, increasing the probability for tunneling through the barrier. Charge injection from a metal electrode into the polymer may be treated using the Fowler-Nordheim quantum mechanical tunneling model given bywhere , , is the electron effective mass, is Planck’s constant (4.1356 × 10−15 eV·s), and is the effective potential barrier. is equal to or for positive and negative charge, respectively. This occurs when the applied field is high enough that the potential barrier seen by the emitted electrons is sufficiently thin to allow tunneling through the barrier. The transition from Schottky emission to Fowler-Nordheim tunneling occurs at about 166 V/μm [18].
5.2. Field-Dependent Mobility
At low-fields and low bulk samples, carriers are almost in equilibrium with the lattice vibrations so the low-field mobility is mainly affected by phonon and Coulomb scattering. The mobility increases until the velocity approaches the random thermal velocity. In a moderately large electric field, less thermal fluctuation is required to free charge allowing for higher conduction via the Poole-Frenkel mobility:where is a constant and is as defined in (27). At higher electric fields, mobility decreases with increasing electric field due to increased lattice scattering at higher carrier energies, and the velocity of carriers begins to saturate. The following expression is used to implement a field-dependent mobility that provides a smooth transition between low-field and high-field behavior:where is the low-field mobility at a field of , is the saturation velocity, and is commonly used [26].
5.3. Charge Attachment/Detachment
Charge trapping takes place at a hopping site that requires energy substantially greater than the average energy to release charge carriers. Trapping and detrapping of space charge in polymeric materials are related to the microstructure and morphology of the materials. Trapping mechanisms include physical defects such as dangling bonds which lead to shallow traps; “self-traps” due to field modification which alters the length of the polymer chain and their potential well; and chemical defects or impurities which result in deep traps. Detrapping mechanisms may be photon-assisted by illumination or phonon-assisted through lattice vibration, impact ionization, and tunneling, with the latter two occurring at high fields.
Charge trapping in the polymer bulk is not considered for now. In accordance with the eEDL model, for particle simulation, charge trapping in/attachment to nanocrystallites are assumed to be deterministic, that is, attachment on impact or probability initially, with subsequent taper-off due to Coulomb force and dipole field repulsion from build-up of the attached charge in accordance with the MWS effect, that is, as the limiting behavior.
First principles molecular dynamics (MD) calculations of trap depths and metal-polymer barrier potentials may be used together with AC-biased curve measurements of leakage current to extract practical transport and attachment/detachment coefficients.
The probability and coefficient for charge detachment/detrapping from nanocrystallites are given bywhere is the detachment cross-section (nm2), is the nanocrystallite density (μm−3), is the trap depth (eV), and is the charge velocity (cm/s). Monte Carlo selection [23] uses random numbers to determine the outcomes of charge detachment events and proceeds as follows.(1)Given .(2)Calculate .(3)Generate random number: .(4)If detachment event.(5)If no detachment.
Detachment is not considered for now due to a lack of published data. However, this model is capable of incorporating future parameters to be abstracted from empirical fits to MD calculations.
5.4. Recombination
Neglecting recombination between mobile charge and trapped bulk charge of the opposite polarity, the Monte Carlo collision model may also be used to describe particle-particle events between oppositely charged entities resulting in recombination. The probability of a collision of the th charge particle with a charge particle of the opposite polarity in a time step is given bywhere is the number density of the opposite polarity mobile charge, is the total collision cross-section which in general depends on the kinetic energy of the th particle, and is the velocity of this particle relative to the velocity of the opposite polarity particle. Therefore, is a measure of the number of collisions per unit length, is the mean free path between collisions, and is the collision frequency. From (32) it is clear that, for , and for , means the particle will eventually collide. For a finite , the probability of a collision is determined by comparing with , representing a uniform random number between 0 and 1. Following the same logic as in charge detachment, for , the particle is assumed to have sustained a collision within the time step resulting in recombination and neutralization of the two charge particles.
6. Results and Discussion
The model for the nanocomposite film is a random distribution of nanofillers and/or nanocrystallites in a 500 nm layer sandwiched on the top and bottom by 250 nm amorphous polymer layers. Nanofillers and/or nanocrystallites are represented as spherical inclusions with ferroelectric properties and much higher permittivity in the case of ceramic nanofillers. The linear dimension of the computational cell is 1 μm. Top and bottom electrodes form the anode and cathode, respectively, between which a bias field is imposed. Only unipolar injection from the anode is considered here. The four vertical side walls have zero flux and periodic boundary conditions where, in the latter, charge particles exiting a side-wall reenter from the opposite side-wall. Charge transport studies are performed for several nanoparticle polarization orientations with respect to the bias field including random, in-plane, parallel, and antiparallel (Figure 2). Table 1 summarizes the simulation parameters, including the dipole moment of BT nanofillers and PVDF nanocrystallites [27].
(a)
(b)
(c)
(d)
Sample trajectories are shown in Figure 3 for randomly distributed 150 nm nanofillers (large circles) at 10 vol% loading with dipole moment of 4.25 × 10−30 C·m and (a) in-plane; (b) parallel; and (c) antiparallel polarization. Figure 3(a) for in-plane polarization shows the rear view with 7 attachments and 3 trajectories passing through to the cathode, one of which suffers multiple deflections enroute. The side view in Figure 3(b) shows 8 trajectories terminating on upper half-spheres due to the bias field and the added attraction of the dipole field and 2 passing through to the cathode. Finally, Figure 3(c) contains a perspective view showing trajectories terminating on the tail-end of the spheres and 1 trajectory that is repelled by the front-end of the dipole, deflected laterally, and subsequently moved to the counter-electrode.
(a)
(b)
(c)
Figures 4(a) and 4(b) shows the charge fractions, normalized to total charge, of those attached to nanofillers and those that arrived at the counter-electrode for the 4 polarization cases. Lower attachment results in higher conduction or higher fractions of arrivals at the counter-electrode and vice versa. In all runs, the tracking time is adjusted so that the same amount of charge is injected. The random and in-plane cases have similar results with the latter case exhibiting a slightly higher attachment due to all dipoles being aligned in the plane to provide equal access to the opposite-signed ends. The parallel case presents the negative end of the dipoles directly to the positive charge and results in the highest attachment. The antiparallel case is just the opposite, making it more difficult for incoming positive charge to be initially repelled before they can be drawn to the underside of the dipoles. Depending on the proximity of the neighboring nanoparticles, a dipole polarization cancellation effect may cause the charge particle to continue migration therefore increasing the probability of propagating through to the counter-electrode. In particular, for the antiparallel case, very little attachment is observed to the nanofillers but there is >90% arrival at the counter-electrode.
(a)
(b)
Shown in Figure 5 are the fractions of positive charge captured on nanofillers, in transit, and arrived at the counter-electrode. The curves asymptote with time to 2.24%, 5.05%, and 92.70%, respectively. Table 2 summarizes the asymptotic fractions for all the polarization cases.
Figure 6 shows the top views of scatter plots for the randomly distributed 150 nm nanofillers (large circles) at 10 vol% loading with dipole moment of 4.25 × 10−30 C·m and attached charge (orange), mobile charge (magenta), and charge arrived at the counter-electrode (blue) for (a) random; (b) in-plane; (c) parallel; and (d) antiparallel polarizations. Figure 6(a) shows mobile positive charge attracted to the tail-end of the randomly polarized spheres. The upwards-directed in-plane polarization in Figure 6(b) shows charge attachment to the back end of the spheres. Figure 6(c) shows mobile positive charge attracted to the tail-end of the sphere, forming fully coated upper hemispheres. Finally, in Figure 6(d) the spherical nanofillers are rendered transparent to show the mobile positive charge attached to the underside or lower hemisphere.
(a)
(b)
(c)
(d)
The combined effects of 10 vol% BT nanofillers and varying degrees of crystallization (10 to 50 vol%) of the PVDF polymer matrix are studied where, for maximum effect, only antiparallel polarization is considered. Figure 7 shows side views of scatter plots and charge trajectories for randomly distributed 150 nm nanofillers (10 vol%, dipole moment of 4.25 × 10−30 C·m) and 70 nm nanocrystallites (10 to 40 vol% loading, dipole moment of 1.183 × 10−30 C·m). Curvilinear trajectories show charge particles migrating through the interspaces under the combined influence of bias and polarization fields and charge. Red trajectories signify attached charge, black trajectories denote unattached (mobile) charge, green circles are PVDF nanocrystallites, and solid blue circles are BT nanofillers. The incident trajectories at low degree of crystallization are almost normal to the plane of the nanocomposite film. With increasing vol% loading of the nanocrystallites, the incident trajectories bend outwards to reflect the increasing repulsion seen by the upward directed polarization. The proportion of red lines is also increased to show increasing attachment with higher degree of crystallization.
Figure 8 shows side views of scatter plots of randomly distributed, antiparallel polarized, 150 nm nanofillers (10 vol%, dipole moment of 4.25 × 10−30 C·m) and 70 nm nanocrystallites (10 to 50 vol% loading, dipole moment of 1.183 × 10−30 C·m) and attached charge (orange), mobile charge (magenta), and charge arrived at counter-electrode (blue). Charge attaches mostly to the tail (negative end) of the dipolar PVDF nanocrystallites and nanofillers. Solid blue circles denote BT nanofillers and green circles are the PVDF nanocrystallites.
Asymptotes of charge fractions over time are tabulated in Table 3 for the range of PVDF nanocrystallite loadings and cases with and without 10 vol% of BT nanofillers. It is observed that the combined effect of nanocrystallites and nanofillers contributes additively to the polarization field resulting in increased attachment up to about 40 vol% of nanocrystallites. For the current set of parameters, leakage conduction appears to asymptote to approximately 55%, which seems to be an inflection point beyond which the fraction arriving at the counter-electrode reincreases. This observation may be nonphysical as the higher predicted values may be attributed to charge leakage to the counter-electrode along the side-walls as seen in the 40 and 50 vol% cases in Figure 8. Therefore, this may be due to the periodic boundary condition that reintroduces exiting particles into the opposite side-wall thereby localizing their trajectories to the edges of the cell. From the 40 and 50 vol% trajectories in Figure 7, higher vol% results in increased repulsion for the antiparallel polarization case, as evidenced by both the deflection of the incident trajectories and the larger proportion of attachments. A clear need is a self-consistent set of empirical measurements to pin down the polarization limit and relate that to the degree of crystallization.
Reduction of the average anode field due to the antiparallel polarized nanoparticles reduces the amount of injected current following the Schottky emission model used here. field reduction also reduces transport through the nanocomposite film layer. Because antiparallel polarization leads to increased leakage current conduction, it is anticipated that there will be reduced accumulation of space charge in the interaction zone and decreased probability for breakdown. Table 4 summarizes the average field reduction on the anodes. The field reduction is approximately 10 V/μm for every 10 vol% increase in degree of crystallinity for this set of parameters.
The charge fractions for 10 vol% BT nanofiller and 30 vol% degree of crystallinity of PVDF is shown in Figure 9 where charge captured on nanocrystallites, nanofillers, in transit, and arrived at the counter-electrode asymptote with time to 34.86%, 7.90%, 4.34%, and 52.70%, respectively.
Figure 10(a) for 10 vol% BT nanofiller and 30 vol% PVDF nanocrystallites shows that average anode and cathode fields are reduced due to the antiparallel polarization. The fields take up to 1 transit time to reach steady-state, with the anode field lowered by the proximity of the injected positive charge. In Figure 10(b), injected current density at the anode stabilizes after an interval of 1 transit time.
(a)
(b)
Shown in Figure 11(a) is the average field at the anode which drops from the bias value of 100 V/μm due to the field from antiparallel poling of ferroelectric BT nanofillers and PVDF nanocrystallites. Field reduction is seen to increase linearly with increasing degree of crystallinity of PVDF. In Figure 11(b), injected current density assumes Schottky thermionic emission which drops due to decreasing anode field. As expected, current density reduction goes exponentially with field reduction.
(a)
(b)
From the preceding results, it is clear that antiparallel polarization has dramatically higher conduction, a property which has immediate implications on reducing residual charge accumulation and extending the lifetime and reliability of nanocomposite films. The mechanism acts to lower field, reduce charge injection, and minimize attachment to nanoparticles by orienting the dipole for field repulsion. Suitable ferroelectric materials include Perovskites, namely, BT, TiO2, and ZrO2 as nanofillers.
Another beneficial implication may be to increase the power conversion efficiency (PCE) of organic photovoltaic (OPV) devices by inserting an ultrathin film of a ferroelectric copolymer, P(VDF-TrFE), at the metal-organic interface to enhance the charge extraction efficiency. Highly crystalline P(VDF-TrFE) films prepared by the Langmuir-Blodgett method spontaneously polarize and has been shown to be responsible for the enhancement of PCE in ferroelectric OPV devices [28].
Too high vol% loading of nanoparticles and/or too strong antiparallel polarization may significantly reduce the electrode field and severely limit charge injection. Results are derived from parameters drawn from the literature. There is therefore a need for verified simulation parameters derived from empirical or first principles calculations to ensure computed results are comparable to laboratory data.
7. Conclusions
This paper has described a rapid 3D particle simulation algorithm which couples the BIEM with robust predictor-corrector time integration of the equations of motion to efficiently simulate charge transport through nanocomposite film comprised of ferroelectric nanofillers and/or nanocrystallites in amorphous polymer matrices. An extended EDL model was implemented which substituted a dipole for the monopole in the classical EDL model. The rationale was to allow for the initial charge particles to attach onto the nanofillers and/or nanocrystallites on impact to form the bound Stern-Helmholtz layer. Subsequent charge particles are repelled due to charge build-up and the MWS effect, leading to the formation of the diffuse Gouy-Chapman transport layer. Careful trajectory tracking revealed the curvilinear paths taken by charge particles that meander through the interspaces that form the interaction zone.
The eEDL model is used to simulate the ferroelectric properties of crystalline nanoparticles and ceramic nanofillers such as BT as statistical ensembles of randomly oriented polar nanoregions with the associated dipole moments embedded in a continuous amorphous matrix. Preliminary simulation results corroborate the notion of modeling the ferroelectric nanoparticles as remanent polarized nanoregions in an amorphous matrix. Extension of the EDL model to ferroelectric simulations shows promise as the predicted rate of charge migrations around both crystalline nanoregions and ceramic nanofillers for different polarization orientations is seen to be field-dependent.
Antiparallel polarization has shown much higher charge conduction. The explanation is the preferential attachment of charge to the tail (negative end) of the dipole due to charge polarity. In the antiparallel case, the mobile charge is repelled before it is attracted to the other end; a situation which may repeat as almost colinear dipole fields of aligned nanoparticles cancel, thus facilitating the charge to move through the interspaces. Charge captured on the oppositely charged end of the dipole may ultimately lead to polarization reversal.
One area of future work is to study the shape-dependence and effects of physical barriers on charge transport. The barrier properties of polymers can be significantly altered by inclusion of inorganic platelets with sufficient aspect ratio to alter the diffusion path of penetrant molecules. These models are generally based on random, parallel filler platelets perpendicular to the permeation direction (random in only two directions), introducing an orientation factor. At high aspect ratio in nanocomposites, significant decreases in permeability are predicted and observed in practice [29]. High aspect ratio nanolayers provide properties that are not possible for larger-scaled composites. The impermeable layers mandate a tortuous pathway to transverse the nanocomposite. The enhanced barrier characteristics and lower charge migration are benefits from the hindered diffusion pathways through the nanocomposite [30].
Another area of future work will involve derivation of transport and attachment/detachment coefficients from empirical data on time-dependent curve measurements of leakage current and matching with model predictions to extract and refine estimates for key parameters. The model will then be used to predict results that can be compared against experiment. Ultimately the validated model may be used for parameter exploration and optimization.
Yet another effort would be to seek corroboration with estimates from first principles quantum mechanical (QM) and molecular dynamics (MD) calculations of trap depths and metal-polymer barrier potentials.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
Support by the Department of the Navy, Office of Naval Research, Grant no. N000141310064, is gratefully acknowledged.