Research Article  Open Access
Meng H. Lean, WeiPing L. Chu, "Model for Charge Transport in Ferroelectric Nanocomposite Film", Journal of Polymers, vol. 2015, Article ID 745056, 17 pages, 2015. https://doi.org/10.1155/2015/745056
Model for Charge Transport in Ferroelectric Nanocomposite Film
Abstract
This paper describes 3D particleincell 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 metalpolymer Schottky emission at low to moderate fields and FowlerNordheim tunneling at high fields. Injected particles propagate via fielddependent PooleFrenkel mobility. The simulation algorithm uses a boundary integral equation method for solution of the Poisson equation coupled with a secondorder predictorcorrector scheme for robust time integration of the equations of motion. The stability criterion of the explicit algorithm conforms to the CourantFriedrichsLevy limit assuring robust and rapid convergence. Simulation results for BaTiO_{3} 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 fielddependent and may also depend on the physicochemical bonds on the common surfaces. Morphology plays an important role in film layers with very high surfacetovolume 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 polymerceramic interfaces to attain a viable solution [1]. Interface engineering attempts to improve the dispersion of filler in polymer matrix. In a welldispersed composite, interfaces between the ceramic nanoparticle and polymer phases create effective electron scattering and transport centers, thus reducing the breakdown probability. Moreover, welldispersed ceramic nanoparticles block degradation tree growth and thus increase the longterm 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 fluoridecotrifluoroethylene), P(VDFTrFE), and poly(vinylidene fluoridecohexafluoropropylene), P(VDFHFP), 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 polymerceramic interfaces. Ferroelectric ceramics, such as BaTiO_{3} (BT), Pb(Mg_{1/3}Nb_{2/3})O_{3}PbTiO_{3} (PMNPT), 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(VDFTrFECFE) terpolymer/ZrO_{2} nanoparticles have been demonstrated through the interface effect with the presence of 1.6 vol% of ZrO_{2} 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 TiO_{2} 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 filmforming properties. A second method is the functionalization or Phosphonic acid surfacemodification of BT nanoparticles to yield welldispersed P(VDFHFP) 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 fluoridecochlorotrifluoroethylene) [P(VDFCTFE)] and poly(vinylidene fluoridetertrifluoroethyleneterchlorotrifluoroethylene) [P(VDFTrFECTFE)] matrices [5]. Another method is chainend functionalization of the ferroelectric polymer and subsequent covalentbonding 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 endgroups of the polymer for directly coupling with oxide nanoparticles to yield the covalentbonded 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]. Highenergydensity polymer nanocomposites based on surfacefunctionalized TiO_{2} nanocrystals used as dopants in a ferroelectric P(VDFTrFECTFE) 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(VDFTrFECTFE) nanocomposites on the TiO_{2} concentration was reported to peak at around 10 vol% TiO_{2} 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 MaxwellWagnerSillars (MWS) interfacial polarization at low frequencies contributing to the creation of the interaction zone with the GouyChapman diffuse layer, thereby greatly affecting polarization and the dielectric response of the polymer matrix. The incorporation of the TiO_{2} 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 thinfilm capacitors. This reduction is observed concurrently with the enhancement of the effective permittivity and breakdown strength of the thinfilm 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 frozenin 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 amorphouscrystalline interface, stable remanent polarization, and the permittivity at 1 kHz [10]. Experimental studies on the dependence of polarization fatigue on crystallinity of P(VDFTrFE) 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 indepth 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 3layer 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 fardistance force originating from the electric double layer (EDL) which is several 100 nm thick, corresponding to the Debye shielding length. Another is the Lewis 2layer model comprising an inner bound Stern layer and a diffuse GouyChapman (DebyeHü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 SternHelmholtz first layer of countercharge tightly attached to the particle surface forming a doublelayer charge structure. The second GouyChapman 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 oppositesigned end of the dipole to form the SternHelmholtz 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 fielddependent mobility, . The MWS effect due to charge buildup 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 MaxwellGarnett, 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 2layer 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 “particleincell” (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 fielddependent mobility under the composite field while allowing for particleparticle and particlenanofiller/nanocrystallite interactions. Nanoregions with ferroelectric properties may be modeled as spherical dipoles with remanent dipole moment and polarization orientation that are random, inplane, 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 SternHelmholtz layer. Subsequent incoming charge may be repelled to form the diffuse outer GouyChapman 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 counterelectrode are curvilinear paths that meander through the interspaces, as shown in Figure 1(c). Charges that arrive at the counterelectrode 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 sphereshaped ferroelectric materials may be written aswhere is the ratio of the permittivity of the sphere to the polymer matrix, is the ClausiusMossotti 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 farfield 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 GaussLegendre quadrature. Singular kernels are accurately computed using minimumorder 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 electrodepolymer interfaces and polarization and trapped charge at the polymerpolymer 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 “scattergather” 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 particlemesh 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 timedependent 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 twostep field solveparticle 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 closepacked (HCP) to achieve a volume fraction of 0.9069. The 3D PIC model simulates the dynamics of charge particle injection, particleparticle 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 SternHelmholtz layer on the opposite charged end of the dipole. Subsequent particles form the diffuse GouyChapman transport layer due to Coulomb and dipole force repulsion. The selfconsistent field solution is obtained using the BIEM with Dirichlet boundary conditions for the electrodes, homogeneous Neumann boundary conditions for the four vertical sidewalls, 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 lowtomoderate fields and FowlerNordheim tunneling at high fields. Charge particles migrate with PooleFrenkel fielddependent mobility subject to velocity saturation at high field. Recombination is based on Monte Carlo statistics. Particle tracking involves a predictorcorrector algorithm to integrate the equations of motion subject to Coulomb, dipole, and bias field forces. A selfconsistent particleparticle, particlemesh (P^{3}M) 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 μm^{3} 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 counterelectrode 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 CourantFriedrichsLevy (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 timedependent 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 predictorcorrector method of order or secondorder 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 predictorcorrector 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 timedependent algorithm loops through the following steps.(1)Inject a number of particles according to Schottky or FowlerNordheim 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 counterelectrode or (b) reintroduce from opposite sidewall if exited from sidewall.(5)If no mobile charge left in polymer, exit. Else return to step (1).
The preceding tacitly assumes that the field is “frozenin” 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 LUMOlevel. 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 HOMOlevel.
At moderate applied fields, electrons or holes within the electrode arrive at the metalpolymer 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/m^{2}·K^{2}) 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 FowlerNordheim 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 FowlerNordheim tunneling occurs at about 166 V/μm [18].
5.2. FieldDependent Mobility
At lowfields and low bulk samples, carriers are almost in equilibrium with the lattice vibrations so the lowfield 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 PooleFrenkel 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 fielddependent mobility that provides a smooth transition between lowfield and highfield behavior:where is the lowfield 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; “selftraps” 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 photonassisted by illumination or phononassisted 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 taperoff due to Coulomb force and dipole field repulsion from buildup 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 metalpolymer barrier potentials may be used together with ACbiased 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 crosssection (nm^{2}), 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 particleparticle 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 crosssection 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 sidewall reenter from the opposite sidewall. Charge transport studies are performed for several nanoparticle polarization orientations with respect to the bias field including random, inplane, 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) inplane; (b) parallel; and (c) antiparallel polarization. Figure 3(a) for inplane 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 halfspheres 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 tailend of the spheres and 1 trajectory that is repelled by the frontend of the dipole, deflected laterally, and subsequently moved to the counterelectrode.
(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 counterelectrode for the 4 polarization cases. Lower attachment results in higher conduction or higher fractions of arrivals at the counterelectrode and vice versa. In all runs, the tracking time is adjusted so that the same amount of charge is injected. The random and inplane 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 oppositesigned 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 counterelectrode. In particular, for the antiparallel case, very little attachment is observed to the nanofillers but there is >90% arrival at the counterelectrode.
(a)
(b)
Shown in Figure 5 are the fractions of positive charge captured on nanofillers, in transit, and arrived at the counterelectrode. 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.
 
10 v% BaTiO_{3}; dipole moment = 4.25 10^{−30} C·m; = 150 nm. 
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 counterelectrode (blue) for (a) random; (b) inplane; (c) parallel; and (d) antiparallel polarizations. Figure 6(a) shows mobile positive charge attracted to the tailend of the randomly polarized spheres. The upwardsdirected inplane 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 tailend 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 counterelectrode (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 counterelectrode reincreases. This observation may be nonphysical as the higher predicted values may be attributed to charge leakage to the counterelectrode along the sidewalls 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 sidewall 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 selfconsistent set of empirical measurements to pin down the polarization limit and relate that to the degree of crystallization.
 
Antiparallel poling of PVDF nanocrystallites and BaTiO_{3} nanofillers. 
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.
 
: field reduction due to antiparallel poling. 
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 counterelectrode 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 steadystate, 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, TiO_{2}, and ZrO_{2} 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(VDFTrFE), at the metalorganic interface to enhance the charge extraction efficiency. Highly crystalline P(VDFTrFE) films prepared by the LangmuirBlodgett 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 predictorcorrector 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 SternHelmholtz layer. Subsequent charge particles are repelled due to charge buildup and the MWS effect, leading to the formation of the diffuse GouyChapman 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 fielddependent.
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 shapedependence 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 largerscaled 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 timedependent 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 metalpolymer 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.
References
 P. Barber, S. Balasubramanian, Y. Anguchamy et al., “Polymer composite and nanocomposite dielectric materials for pulse power energy storage,” Materials, vol. 2, no. 4, pp. 1697–1733, 2009. View at: Publisher Site  Google Scholar
 M. Poulsen and S. Ducharme, “Why ferroelectric polyvinylidene fluoride is special,” IEEE Transactions on Dielectrics and Electrical Insulation, vol. 17, no. 4, pp. 1028–1035, 2010. View at: Publisher Site  Google Scholar
 Q. Wang and L. Zhu, “Polymer nanocomposites for electrical energy storage,” Journal of Polymer Science, Part B: Polymer Physics, vol. 49, no. 20, pp. 1421–1429, 2011. View at: Publisher Site  Google Scholar
 P. Kim, N. M. Doss, J. P. Tillotson et al., “High energy density nanocomposites based on surfacemodified BaTiO_{3} and a ferroelectric polymer,” ACS Nano, vol. 3, no. 9, pp. 2581–2592, 2009. View at: Publisher Site  Google Scholar
 J. Li, J. Claude, L. E. NorenaFranco, S. I. Seok, and Q. Wang, “Electrical energy storage in ferroelectric polymer nanocomposites containing surfacefunctionalized BaTiO_{3} nanoparticles,” Chemistry of Materials, vol. 20, no. 20, pp. 6304–6306, 2008. View at: Publisher Site  Google Scholar
 J. Li, P. Khanchaitit, K. Han, and Q. Wang, “New route toward highenergydensity nanocomposites based on chainend functionalized ferroelectric polymers,” Chemistry of Materials, vol. 22, no. 18, pp. 5350–5357, 2010. View at: Publisher Site  Google Scholar
 J. Li, S. I. Seok, B. Chu, F. Dogan, Q. Zhang, and Q. Wang, “Nanocomposites of ferroelectric polymers with TiO_{2} nanoparticles exhibiting significantly enhanced electrical energy density,” Advanced Materials, vol. 21, no. 2, pp. 217–221, 2009. View at: Publisher Site  Google Scholar
 M. N. Almadhoun, U. S. Bhansali, and H. N. Alshareef, “Nanocomposites of ferroelectric polymers with surfacehydroxylated BaTiO_{3} nanoparticles for energy storage applications,” Journal of Materials Chemistry, vol. 22, no. 22, pp. 11196–11200, 2012. View at: Publisher Site  Google Scholar
 D. Rollik, S. Bauer, and R. GerhardMulthaupt, “Separate contributions to the pyroelectricity in poly(vinylidene fluoride) from the amorphous and crystalline phases, as well as from their interface,” Journal of Applied Physics, vol. 85, no. 6, pp. 3282–3288, 1999. View at: Publisher Site  Google Scholar
 A. B. Silva, R. Gregorio, J. V. Esteves, and C. Wisniewski, “Influence of the amorphouscrystalline interface on the dielectric and ferroelectric polarization of αPVDF,” in Proceedings of the 11th International Conference on Advanced Materials, Rio de Janeiro, Brazil, 2009. View at: Google Scholar
 Z.G. Zeng, G.D. Zhu, L. Zhang, and X.J. Yan, “Effect of crystallinity on polarization fatigue of ferroelectric P(VDFTrFE) copolymer films,” Chinese Journal of Polymer Science (English Edition), vol. 27, no. 4, pp. 479–485, 2009. View at: Publisher Site  Google Scholar
 D. Pitsa and M. G. Danikas, “Interfaces features in polymer nanocomposites: a review of proposed models,” Nano, vol. 6, no. 6, pp. 497–508, 2011. View at: Publisher Site  Google Scholar
 R. Mekala and N. Badi, “Modeling and simulation of high permittivity coreshell ferroelectric polymers for energy storage solutions,” in Proceedings of the COMSOL Conference, Boston, Mass, USA, 2013. View at: Google Scholar
 I. A. Tsekmes, R. Kochetov, P. H. F. Morshuis, J. J. Smit, and T. Andritsch, “Modeling the dielectric response of epoxy based nanocomposites,” in Proceedings of the IEEE Electrical Insulation Conference (EIC '14), pp. 47–50, 2014. View at: Google Scholar
 Y. P. Mamunya, V. V. Davydenko, P. Pissis, and E. V. Lebedev, “Electrical and thermal conductivity of polymers filled with metal powders,” European Polymer Journal, vol. 38, no. 9, pp. 1887–1897, 2002. View at: Publisher Site  Google Scholar
 J. Hicks, A. Behnam, and A. Ural, “A computational study of tunnelingpercolation electrical transport in graphenebased nanocomposites,” Applied Physics Letters, vol. 95, no. 21, Article ID 213103, 2009. View at: Publisher Site  Google Scholar
 M. H. Lean and W.P. L. Chu, “Effect of gaseous void on bipolar charge transport in layered polymer film,” Journal of Physics D: Applied Physics, vol. 47, no. 7, Article ID 075303, 2014. View at: Publisher Site  Google Scholar
 M. H. Lean and W.P. L. Chu, “Dynamic charge mapping in layered polymer films,” IEEE Transactions on Dielectrics and Electrical Insulation, vol. 21, no. 3, pp. 1319–1329, 2014. View at: Publisher Site  Google Scholar
 M. H. Lean and W.P. L. Chu, “Simulation of charge packet formation in layered polymer film,” COMPEL: The International Journal for Computation and Mathematics in Electrical and Electronic Engineering, vol. 33, no. 4, pp. 1396–1415, 2014. View at: Publisher Site  Google Scholar  MathSciNet
 M. H. Lean and W.P. L. Chu, “Dielectric breakdown model for polymer films,” IEEE Transactions on Dielectrics and Electrical Insulation, vol. 21, no. 5, pp. 2259–2266, 2014. View at: Publisher Site  Google Scholar
 M. H. Lean, “Particle simulations of ion cloud in a magnetic field,” IEEE Transactions on Magnetics, vol. 34, no. 5, pp. 3118–3121, 1998. View at: Publisher Site  Google Scholar
 R. W. Hockney and J. W. Eastwood, Computer Simulation Using Particles, McGrawHill, New York, NY, USA, 1981.
 J. M. Dawson, “Particle simulation of plasmas,” Reviews of Modern Physics, vol. 55, no. 2, pp. 403–447, 1983. View at: Publisher Site  Google Scholar
 C. K. Birdsall, “Particleincell chargedparticle simulations, plus Monte Carlo collisions with neutral atoms, PICMCC,” IEEE Transactions on Plasma Science, vol. 19, no. 2, pp. 65–85, 1991. View at: Publisher Site  Google Scholar
 D. M. Caughey and R. E. Thomas, “Carrier mobilities in silicon empirically related to doping and field,” Proceedings of the IEEE, vol. 55, no. 12, pp. 2192–2193, 1967. View at: Publisher Site  Google Scholar
 S. Gottlieb, C.W. Shu, and E. Tadmor, “Strong stabilitypreserving highorder time discretization methods,” SIAM Review, vol. 43, no. 1, pp. 89–112, 2001. View at: Publisher Site  Google Scholar  MathSciNet
 C.G. Duan, W. N. Mei, W.G. Yin et al., “Simulations of ferroelectric polymer film polarization: the role of dipole interactions,” Physical Review B: Condensed Matter and Materials Physics, vol. 69, no. 23, Article ID 235106, 2004. View at: Publisher Site  Google Scholar
 Y. Yuan, P. Sharma, Z. Xiao et al., “Understanding the effect of ferroelectric polarization on power conversion efficiency of organic photovoltaic devices,” Energy and Environmental Science, vol. 5, no. 9, pp. 8558–8563, 2012. View at: Publisher Site  Google Scholar
 D. R. Paul and L. M. Robeson, “Polymer nanotechnology: nanocomposites,” Polymer, vol. 49, no. 15, pp. 3187–3204, 2008. View at: Publisher Site  Google Scholar
 P. C. Lebaron, Z. Wang, and T. J. Pinnavaia, “Polymerlayered silicate nanocomposites: an overview,” Applied Clay Science, vol. 15, no. 12, pp. 11–29, 1999. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2015 Meng H. Lean and WeiPing L. Chu. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.