We report a model describing the molecular orientation disorder in CH3NH3PbI3, solving a classical Hamiltonian parametrised with electronic structure calculations, with the nature of the motions informed by ab initio molecular dynamics. We investigate the temperature and static electric field dependence of the equilibrium ferroelectric (molecular) domain structure and resulting polarisability. A rich domain structure of twinned molecular dipoles is observed, strongly varying as a function of temperature and applied electric field. We propose that the internal electrical fields associated with microscopic polarisation domains contribute to hysteretic anomalies in the current-voltage response of hybrid organic-inorganic perovskite solar cells due to variations in electron-hole recombination in the bulk.
Solar cells based on hybrid organic-inorganic perovskites display unusual device physics, which is still in the process of being understood.1–5 One unusual aspect is the notable hysteresis in current–voltage curves, depending on rate of measurement and starting point within the curve.3 When the solar cell is held at short-circuit (in the dark or light), the photovoltaic performance decreases considerably. When the photovoltaic cell is operating at near open-circuit voltage, performance builds.
The highest apparent efficiencies are produced after tempering the solar cell in forward bias at 1.4 V.3 Similar behaviour is observed in both mesoporous and planar solar cell architectures and independent of the presence of a hole transport layer (e.g., SpiroMeOTAD). This effect is strongest in planar device architectures, as would be expected if it resulted from the electric field applied between the contacts. This behaviour occurs over time scales up to hundreds of seconds.
Perovskite (ABX3) structured materials tend to have a large dielectric constant due to the relative ease of polarising the cell structure. In particular, distortion of the edge-sharing BX6 octahedra can produce an overall electric dipole between the A and B lattice sites. Inorganic perovskites are known to exhibit a range of ferroelectric and ferroelastic phase transitions.6 The hybrid perovskite analogues are formed by replacing the element at the icosahedral A site with an isovalent molecule. In this work, we study the photovoltaic absorber methyl-ammonium lead iodide, CH3NH3PbI3 (MAPI). Here
Ferroelectric materials have been studied in the context of photovoltaic applications for over half a century.7 A number of effects have been attributed to the lack of centrosymmetry. The anomalous photovoltaic effect was first reported in samples of PbS8 in 1946 and was subsequently linked to the existence of a ferroelectric transition, resulting in large photovoltages (ca. 10 kV).9 Similar effects have been reported for ferroelectric phases of SbSI, ZnS, and CdTe.9–11 Additionally, oxide perovskites such as BaTiO3 and KNbO3 demonstrate a bulk photovoltaic effect.12 In these materials photocurrents can be generated in the absence of asymmetric electrical contacts, unlike standard photovoltaic cells.
In this Letter, we report the implementation of a classical Monte-Carlo simulation of the domain behaviour arising from molecular rotation in hybrid perovskites. The model takes parameters from density functional theory (DFT) calculations, using both static lattice and molecular dynamic (MD) simulations. In this work, we restrict ourselves to the two-dimensional case, and allow the dipoles to freely rotate (no cage strain is applied).
We propose that the MA ions in this material are highly rotationally mobile, the interaction between these ions forms ordered domains (which respond slowly to applied electric fields), which results in a structured local potential field. We speculate that the equilibrium (open-circuit) configuration is beneficial for solar cell operation by reducing charge carrier recombination through interpenetrating percolating pathways of (electric potential) valleys and ridges for holes and electrons. At short-circuit, the electric field resulting from the built-in voltage is sufficient to disrupt this structure, suppressing long-range order and resulting in more isolated domains. The existence of intricate dipole phase behaviour and the resultant structure in internal electric fields indicate that these photoferroic characteristics, atypical of standard photovoltaic materials, must be considered in device modelling.
Low-frequency dielectric behaviour: The dielectric response of CH3NH3PbI3, and related hybrid perovskites, exhibits significant temperature and frequency dependence. At low temperatures there is a discontinuity associated with the first-order phase transition between the orthorhombic and tetragonal phases (ca. 160 K); at higher temperatures the orientation of MA becomes (partially) disordered. Work by Poglitsch and Weber in 1987 measured the complex dielectric response of methyl-ammonium lead halides (iodide, chloride, bromide) as a function of temperature between 100 and 300 K.13 The effective dielectric constant at 300 K was measured to be 33 for CH3NH3PbI3 at a frequency of 90 GHz. In 1992, Onoda-Yamamuro et al. reported a value of ca. 58 at a frequency of 1 kHz.14 In contrast, the static dielectric constant, in the absence of molecular reorientation, is predicted to be 24.1 from electronic structure calculations (PBEsol + QSGW),15 which is in good agreement with the value of 23.3 determined from a fit of permittivity measurements over 100–300 K to the Kirkwood-Fröhlich equation.14 The unusual dielectric behaviour will make analysis of impedance measurements on photovoltaic cells challenging.
Ab initio molecular dynamics: An open question in these materials is the alignment and dynamics of the MA ion. Analysis of 2H and 14N NMR spectra confirmed that MA cation rotation is a rapid process at room temperature.16 X-ray diffraction has been used to characterise the low temperature orthorhombic (Pna21, C2v symmetry), room temperature tetragonal (I4/mcm, D4h symmetry), and above room temperature cubic (Pm3m, Oh symmetry) crystal structures of MAPI.13 The position of the MA molecules is usually described with partial occupancies that satisfy the space group symmetry,13,17 e.g., in the cubic phase eight identical positions can be fitted around the standard (
To provide atomistic insight, without the assumption or restriction of lattice symmetry beyond a periodic supercell (80 atom 2 × 2 × 2 expansion of the pseudo-cubic perovskite structure21), we investigate the energetics and dynamics of MA in MAPI with ab-initio MD simulations based on the PBEsol exchange-correlation functional. We employ a timestep of 0.5 fs, which is sufficient to describe even the C–H vibrations. The MD trajectory at 300 K contains large-scale fluctuation of the ions about their equilibrium positions, including rapid rotation of the methyl group and total rotation of the methyl-ammonium ion. (A video file of the trajectory is provided elsewhere.22) The hybrid perovskites are structurally soft materials; so far almost all published calculations and analysis have assumed perfect crystals, while these data indicate that such structures are not representative of the materials at room temperature. Further structural analysis is on-going.
Due to the permanent molecular dipole of the methyl-ammonium ion,5 its ensemble average position and the dynamics of its movements are of interest in explaining the dielectric response and electrical behaviour of devices made from MAPI. Custom codes were written to analyse the MD trajectories, identifying the C–N bonds across the periodic boundaries and calculating this molecular alignment (of the eight MA ions) relative to the pseudo-cubic unit cell. The distribution of spherical coordinates over the MD ensemble enables us to make statements about the average distribution of molecular direction, relative to the crystallographic unit cell. Here the histogram is in binned in (θ, ϕ) spherical angles, and hence the bins are not a constant solid angle. The Oh symmetry of the ideal cubic perovskite phase leads to a 48-fold reduction of the phase space onto its reflection domain.
When we plot the data without considering the symmetry of the MA ion environment (Figure 1(a)), there is little that we can say other than a preference for the ion to align with the faces of the cube (Δθ = π/2). The limited simulation time leads to the evident incomplete coverage of spherical phase space. Therefore, we reflect the data onto the first octant, and further exploit symmetry to reduce the internal coordinates to contain the domain between the unique [1, 0, 0] (X) faces, [1, 1, 0] (M) edges, and [1, 1, 1] (R) diagonals (Figure 1(b)) to increase the signal to noise ratio by a factor of 48 (see the supplementary material23 for more details). This reveals a high density of ensembles in a distribution around facial alignment, a lowered distribution around edge alignment, and an increased distribution in a disordered halo around diagonal alignment.
![FIG. 1. Density maps (2D histograms, in spherical coordinates) of MA alignment within the pervoskite cage structure as determined by ab initio molecular dynamics at 300 K. The data are centered on ϕ, θ = 0 being facial orientation. The symmetry folded data bounds the segment between diagonally aligned in the cube (bottom right), pointing at an edge (top right) and pointing at a face (top left). This covers a segment of the original data bounded by 0 < θ < π/4 and \documentclass[12pt]{minimal}\begin{document}$arcos(\frac{1}{\sqrt{3}}) < \phi < \pi /2$\end{document}arcos(13)<ϕ<π/2. In total 18547 data points (8 MA alignments per data frame) are hexagonally binned, with 36 bins in the x-axis. (a) Spherical coordinates, no symmetry, (b) Full 48-fold (Oh) symmetry. Face, edge and diagonal refer to the alignment of the molecule with respect to the perovskite cage, which correspond to the ⟨100⟩, ⟨110⟩ and ⟨111⟩ directions.](/view-large/figure/65136400/081506_1_f1.jpg)
Density maps (2D histograms, in spherical coordinates) of MA alignment within the pervoskite cage structure as determined by ab initio molecular dynamics at 300 K. The data are centered on ϕ, θ = 0 being facial orientation. The symmetry folded data bounds the segment between diagonally aligned in the cube (bottom right), pointing at an edge (top right) and pointing at a face (top left). This covers a segment of the original data bounded by 0 < θ < π/4 and
Density maps (2D histograms, in spherical coordinates) of MA alignment within the pervoskite cage structure as determined by ab initio molecular dynamics at 300 K. The data are centered on ϕ, θ = 0 being facial orientation. The symmetry folded data bounds the segment between diagonally aligned in the cube (bottom right), pointing at an edge (top right) and pointing at a face (top left). This covers a segment of the original data bounded by 0 < θ < π/4 and
We can further quantify these distributions by binning the ensemble of symmetry reduced alignment vectors by whether they are nearest (in angle) to the face, edge, or diagonal vectors. Doing so we find that the raw densities are 35% face, 42% edge, and 23% diagonally aligned (populations [6497, 7822, 4228]). Due to the symmetry of these orientations (with 6, 12, and 8 fold degeneracy, respectively), and the boundaries between these domains, these populations are not directly comparable (the solid angles they cover are different). Therefore, we weight these distributions with surface areas evaluated from a flat spherical distribution (calculated with a Monte-Carlo integration of 105 points, using the same codes). These densities, which would be equal if there was no orientation preference, are 42% face, 31% edge, and 26% diagonal, aligned.
While ab initio MD is a powerful approach, the two principle limitations due to computational expense are finite size effects (as the simulation is infinitely periodic on ≈nm) and short timescale (≈ps are insufficient for diffusion processes). To access the time and length scales necessary to represent realistic non-equilibrium structures, and so directly compare to experiment, we construct a classical model for molecular dipole interactions.
Polar molecules on a lattice: We start from the lattice dynamical theory of ferroelectricity (see Anderson24 and Cochran.25,26) We limit ourself to assuming that the dominant soft phonon mode in the system is the free rotation of the molecular dipole within an extended perovskite cage structure.
The treatment of polarisation as an effect of rotational Brownian motion is analytically challenging.27 Here we simulate these physics numerically by using a Monte Carlo method to calculate the equilibrium configuration of the dipoles. The macroscopic response of the material is expected to be linked to very slow rearrangements of domain walls as a result of statistically rare cooperative rearrangements of the microscopic dipoles. Therefore, we need to integrate a long way in simulation time to equilibrate the material.
We construct a model Hamiltonian for the dipoles (vectors pi) by summing the interaction energy of the applied unshielded electrostatic field (E0), near-neighbour dipole-dipole interactions, and local cage strain (K),
Here the energies are calculated with numerically efficient dot products operating on Cartesian three-vectors. The unit vector
At equilibrium, we associate an electric displacement D related to E0 and the polarisation density by
The polarisation density P can be calculated by a summation over the microscopic dipoles. εs refers to the static relative permittivity, rather than dielectric constant, as it is a function of temperature
We can reconstruct the dipole potential felt at an arbitrary lattice site by summing the potential contribution from all other lattice sites
The main simulation variables are the strength of the interactions. Considering a point dipole-dipole interaction between unshielded methyl-ammonium dipoles (2.29D)5 at the unit cell spacing of 6.29 Å, this energy is ≈25 meV, which we take in this work as exactly 1kBT (T = 300 K). From here we take K = 0, allowing the dipoles to freely rotate without frustration.
Equilibrium behaviour at room temperature: The domain structure of an equilibrated film with zero applied field is shown in Figure 2. Twinning of the MA dipoles occurs to minimise the free energy. This leads to aligned domains along the square axes (as cofacial alignment minimises the dipole-dipole distance compared to diagonal). Visualising the resulting dipole potential that comes from this alignment, we observe the presence of structured interpenetrating regions of high and low electric potential, following the features in the domain boundaries.
(Top) Representative orientation of dipoles in a 25 × 25 two-dimensional periodic slab representing CH3NH3PbI3, initial random configuration (left) and after equilibration under Monte Carlo at 300 K (right). (Bottom) Electrostatic potential as a result of alignment of dipoles, calculated by summation over all other lattice sites for each site. Note the emergence of twinned domain structure formed of linear dipole chains. In the dipole potential, these give rise to structured interpenetrating regions of positive and negative electric potential. (a) Random starting configuration, (b) After equilibration.
(Top) Representative orientation of dipoles in a 25 × 25 two-dimensional periodic slab representing CH3NH3PbI3, initial random configuration (left) and after equilibration under Monte Carlo at 300 K (right). (Bottom) Electrostatic potential as a result of alignment of dipoles, calculated by summation over all other lattice sites for each site. Note the emergence of twinned domain structure formed of linear dipole chains. In the dipole potential, these give rise to structured interpenetrating regions of positive and negative electric potential. (a) Random starting configuration, (b) After equilibration.
The temperature dependence of the domain structure is shown in Figure 3. At 0 K, a striped anti-ferroelectric phase is favoured, which becomes increasingly disordered as the temperature increases due to the role of configuration entropy. The room temperature phase could be viewed as superparaelectric, consisting of randomly oriented linear ferroelectric domains, while the domains are broken to give a paraelectric phase at 1000 K.
Methyl-ammonium dipole alignment as a function of temperature. Low temperature results in complete alignment and large domains, which disorder with increasing temperature leading to small domains, eventually resulting in completely disordered dipoles in the high temperature limit. The physical temperature of these transitions is difficult to quantify, being linear in the dipole-dipole interaction potential (here 25 meV) (a) 0 K, (b) 100 K, (c) 300 K, (d) 1000 K.
Methyl-ammonium dipole alignment as a function of temperature. Low temperature results in complete alignment and large domains, which disorder with increasing temperature leading to small domains, eventually resulting in completely disordered dipoles in the high temperature limit. The physical temperature of these transitions is difficult to quantify, being linear in the dipole-dipole interaction potential (here 25 meV) (a) 0 K, (b) 100 K, (c) 300 K, (d) 1000 K.
The effective simulation temperature varies linearly with respect to the Hamiltonian interaction energies. As our model currently has a point-dipole approximated interaction energy, and ignores energetic contributions from the cage strain, or ion inertia (freely rotating dipoles), simulation temperature cannot be expected to correspond directly to physical temperature.
Electric field dependence of polarisation: For a solar cell to operate the electrical contacts must be selective—a difference in work function must exist between the front and back contacts. This selectivity induces a built-in potential that at short-circuit (or in the dark) results in an electric field across the device, which acts to sweep out generated charge. Current collection at short-circuit is therefore generally maximised. Here we assume that the perovskite solar cells, without intentional doping, are p − i − n type; the potential drops linearly across the active material, producing a constant electric field. When a voltage is applied in forward bias, it works to counteract this built-in field. At the maximal power point for a relatively optimal solar cell material such, as MAPI, the operating voltage is close to the open-circuit voltage, which is close to the built-in potential and so only a small electrical field will apply across the device.
In typical perovskite solar cells, the layer thickness is of order hundreds of nanometres, and the built in voltage of ≈1 V. In the absence of charge equalisation effects, this results in an electric field of 1–10 MV m−1 across the hybrid perovskite. The interaction energy of the MA dipole with this field is U = −p.E, 0.48 meV for a upper limit field strength of 10 MV m−1. This is a relatively small perturbation compared to the dipole-dipole interaction of 25 meV.
The alignment of the dipoles as a result of the applied field is shown in Figure 4; the dipoles rearrange to partially counteract the applied field. The response in overall dipole alignment is 0.5% in the direction of the field (versus 0.04% background fluctuation for no-applied field. We emphasise that it took 104 Monte Carlo moves per lattice site to achieve this equilibrium structure, well beyond where total energy and Debye polarisation appeared to have approached their equilibrium asymptotes.
Analysis of the dipole alignment in a 25 × 25 CH3NH3PbI3 pervoskite film, equilibrated at open and short-circuit conditions. (Top) Alignment of the dipoles, (Middle) dipole potential, (Bottom) 2D Fourier transform of dipole potential (zero frequency component shifted to the centre), showing the change in periodicity of the structure. (a) Open-circuit (no applied field), (b) Short-circuit (field 10 MV m−1).
Analysis of the dipole alignment in a 25 × 25 CH3NH3PbI3 pervoskite film, equilibrated at open and short-circuit conditions. (Top) Alignment of the dipoles, (Middle) dipole potential, (Bottom) 2D Fourier transform of dipole potential (zero frequency component shifted to the centre), showing the change in periodicity of the structure. (a) Open-circuit (no applied field), (b) Short-circuit (field 10 MV m−1).
Though such a small perturbation from the built-in field exhibits negligible visual effect on the alignment (Figure 4 – Top), the effect on the dipole potential (Figure 4 – Middle) is strong, leading to deeper more segregated regions of positive and negative potential.
We quantitatively evaluate the change in the structure of the dipole potential by a two-dimensional Fourier transform (Figure 4 – Bottom). Here we see that the zero field (open-circuit) equilibrium structure is equally distributed in both axes, and the density along the origin indicating the presence of linear features. We interpret these short-circuit features as being the development of carrier traps by the dipole domain response to the built-in field; at open-circuit the extended ridges and valleys in potential could act as channels in which charge transport can take place. Understanding the true role of these features and quantifying the effect on device performance will require a sophisticated device model and an improved understanding of the nature of charge carriers and charge transport in this class of materials.
Particularly, the size of the carrier polarons in these materials will heavily affect the influence of such local inhomogeneiity in the electrostatic potential on charge transport. The Fröhlich coupling constant is estimated to be α = 1.2 using published band parameters for MAPI (m* = 0.12, ε0 = 24.1, ε∞ = 4.5) and a longitudinal optical phonon frequency of 9 THz. This indicates an intermediate electron-lattice coupling, phonon dressing leading to a 25% increase in the effective mass by the Feynman variational treatment.28 The electron and hole polaron radii would correspond to approximately 5 perovskite unit cells, sufficiently small to be influenced by inhomogeneity in the local electrostatic potential we predict. We therefore consider it plausible that such variations in electrostatic potential will drive both carrier polaron segregation (open-circuit dipole structure) and trapping (short-circuit dipole structure), leading to increased bulk recombination at short-circuit.
Recently, we have been made aware of inverted MAPI perovskite cell designs capped with a fullerene electron accepting layer, which exhibit reduced hysteresis.29 We interpret this as a result of gaps in the MAPI film resulting in high penetration of the fullerene, forming a heterojunction. This effectively quenches the recombination, as electron extraction into the fullerene phase out competes recombination. In comparison to fullerene, TiO2 is a poor (slow) electron acceptor.
Additional ferroelectric contributions: Beyond molecular dipole reorientations, additional ferroelectric contributions include distribution of free carriers (electrons and holes), as well as rotations and titling of the
Under short-circuit conditions the MAPI layer will be polarised due to the alignment of dipoles, as demonstrated in our simulations. If the system is then placed in open-circuit conditions, the polarisation of MAPI is removed; however, the depolarisation field consisting of charge carriers will take time to re-equilibrate. There will necessarily be feedback between the ferroelectric domain structure (slower process) and carrier distribution (faster process), therefore causing a further hysteretic contribution to current-voltage measurements. Additionally, the electric field across the absorber (under short-circuit conditions) would lead to an alignment of the cage polarisation (rotation and titling) and domain structure present in the hybrid perovskite film, which is likely to be over a longer time scale than the dipole reorientation. The realignment of these domains upon removal or reversal of the field is another possible source of hysteresis.
In this work, we have shown through ab initio molecular dynamics that the methyl-ammonium ion is rotationally mobile in hybrid perovskites at room temperature, and that the material is structurally soft. This material behaviour may be fortuitous in terms of facilitating transport across grain boundaries when combined with the calculated large polaron transport and small effective masses.5 The large site variation of the ions deserves further study in terms of its effect on material polarisation and ferroelectric response. Further investigation of the molecular dynamics will include expanding the simulation volume and analysing the trajectories further. Additional work with the Monte Carlo codes are required to extend the simulation to three dimensional perovskite volumes, introduce other move types (such as movement of ions both within the lattice cells, and as net migration through the film), and extend the simulated experiments to other relevant device physics. More sophisticated interaction terms for the molecular cations would increase the expected accuracy of these codes.
In conclusion, we have investigated the behaviour of the dipolar methyl-ammonium cation in CH3NH3PbI3 using numerical simulations, which have provided insights into the domain structure and polarisation fields, which will be important for developing quantitative models to explain the unusual device physics of hybrid perovskite solar cells.
Computational details: The set-up for density functional theory (DFT) calculations of the primitive unit cells of CH3NH3PbI3 with a range of molecular orientations, including structure optimisation and static dielectric response,15,21 as well as the lattice polarisation and barriers to rotation,5 have been previously reported. These were taken as the starting point for the MD simulations in this study.
Molecular dynamics (MD): Finite temperature Newtonian MD simulations were performed based on the atomic forces calculated at each timestep using DFT (Γ-point sampling with the PBESol functional and a 400 eV plane wave cutoff). The starting configuration was a fully relaxed 2 × 2 × 2 pseudo-cubic supercell (with a 3 × 3 × 3 k-grid) with MA ions aligned along the ⟨100⟩ direction. Spin-orbit coupling is not treated primarily due to the prohibitive computational expense; as orbital occupation is not changed for the undoped system, the effect on the atomic forces is expected to be negligible. Trajectory data were collected every 50 integration steps (25 fs). A Nosé thermostat (canonical ensemble) was used with a Nosé mass of 3. Custom codes were written for the analysis, with the help of the MDAnalysis library.33 A total of 58 ps (2319 frames) of data was used for analysis, after an equilibration run of 5 ps. This generated 18547 unique MA alignment vectors.
Monte Carlo (MC): The MC implementation uses a Mersenne Twister34 pseudo-random number generator; the code is serial (running at 106 s−1 on modest hardware), but efficient uses of computational resources is achieved by making ensemble runs with GNU Parallel.35 The modestly sized simulations presented here can be extended up to full device sized simulations, with well defined statistics over ensemble runs. The initial state is a lattice of randomised dipoles. The resulting classical Hamiltonian dipole Monte Carlo code, StarryNight, is available as a source code repository on GitHub.36 Codes to interpret the ab initio molecular dynamics used in the production of this paper are similarly available.37
We are grateful for useful discussions with Piers Barnes and Aurélien Leguy (Imperial College London), chiefly concerning the role of the electric field in these materials. The molecular dynamics in this work was instigated to interpret their (unpublished) neutron scattering data. We thank Laurie Peter and Petra Cameron (University of Bath) for useful discussions on hysteresis in perovskite solar cells. We acknowledge membership of the UK’s HPC Materials Chemistry Consortium, which is funded by EPSRC Grant No. EP/F067496. Additional computing resources were provided via the PRACE project UltraFOx. J.M.F. and K.T.B. are funded by EPSRC Grant Nos. EP/K016288/1 and EP/J017361/1, respectively. A.W. acknowledges support from the Royal Society and ERC (Grant No. 277757).
