Connecting experimental synthetic variables with the microstructure and electronic properties of doped ferroelectric oxides using high-throughput frameworks
Doping remains as the most used technique to photosensitize ferroelectric oxides for solar cell applications. However, optimizing these materials is still a challenge. First, many variables should be considered, for instance dopant nature and concentration, synthesis method or temperature. Second, all these variables should be connected with the microstructure of the solid solution and its optoelectronic properties. Here, a computational high-throughput framework that combines Boltzmann statistics with DFT calculations is presented as a solution to accelerate the optimization of theses materials for solar cells applications. This approach has two main advantages: i) the automatic and systematic exploration of the configurational space and ii) the connection between the changes in the microstructure of the material and its electronic properties. One of the most studied doped-ferroelectric systems,x, is used as a study case. Our results not only agree with previous theoretical and experimental reports, but also explain the effect of some of the variables to consider when this material is synthesized.
I. INTRODUCTION<!>II. METHODOLOGY A. Framework<!>Unique cells Degeneracy<!>Individual porperties<!>B. Computational details<!>One of the first systems that attracted the attention as a candidate for use as ferroelectric-based photovoltaic devices was [KNbO<!>CONFLICTS OF INTEREST
<p>Since 2009, hybrid perovskites have revolutionized solar cell research, obtaining power conversion efficiencies, PCE, up to 22% in less than 10 years [1][2][3]. In contrast to conventional silicon solar cells which are based on p − n junctions, the efficiency of this new generation of solar cells is not controlled by the Shockley-Queisser limit [4]. Ferroelectric materials exhibit a spontaneous electric polarization, facilitating the excited carrier separation, which is known as bulk photovoltaic effect, BPVE [5][6][7]. However, there are two aspects that hamper their fast implementation and deployment in the market: i) stability or durability and ii) environmental impact. The most well-studied example is CH 3 NH 3 PbI 3 or MAPbI 3 , which is not very stable in the presence of ambient moisture and also contains Pb [8][9][10]. As an alternative to hybrid perovskite thin-films, free-Pb inorganic ferroelectric perovskites combine a strong spontaneous electric polarization to separate excited carriers, with a wide range of thermal, chemical and mechanical stabilities. In this group, ferroelectric oxide perovskites, ABO 3 , are considered the most promising candidates because of the availability of low-cost methods for their synthesis [11,12].</p><p>While ferroelectric oxides perovskites excel at separating effectively excited charges, they are far from being good sunlight harvesters. Different strategies have been proposed to tailor the electronic properties of these oxides, reducing the band gap from 2 eV-3 eV to 1.5 eV [13]. The epitaxial growth of strained thin-film oxides [14][15][16][17] and promoting the rhombohedral-to-tetragonal phase transition [18] are both effective approaches to reduce the band gap of these oxides. However, doping of both A and B sites is the most common route for modifying their electronic properties [19][20][21][22][23][24][25][26]. Ferroelectric perovskites such as KNbO 3 [27][28][29], BiFeO 3 [30][31][32][33][34], or PbTiO 3 [35][36][37][38][39][40] have been doped by substitutions of their A and B sites to reduce the band gap and increase their photocurrent and PCE.</p><p>Theoretical works, most of them based on density functional theory, DFT, calculations, have helped to understand doping effects on the electronic and structural properties of ferroelectric oxides [41,42]. For instance, A-site substitutions modify the valence band through their interaction with the oxygen atoms [43]. Although these substitutions reduce the polarization of the solid [19], they can produce a direct-to-indirect transition of the band gap, reducing the probability of photoexcited carrier recombinations [43]. Contrary to A-sites, B-site substitutions alter the topology of the conduction band or create midgap states that remarkably reduce the band gap [13]. The information that can be extracted from these calculations is extremely valuable in order to guide the synthesis and characterization of new doped ferroelectric oxides, but the number of variables to be considered requieres the use of systematic approaches.</p><p>Chemical composition, dopant content, stoichiometry, defects or microstructure are some of the variables to consider when ferroelectric doped perovskites are designed. High-throughput frameworks are powerful tools to explore high-dimensional spaces [17]. Very recently, these tools have been applied to predict the electronic and thermochemical properties of doped perovskites [44][45][46]. However, these studies are focused on exploring a wide range of metallic dopants and, most of them, considering a random distribution of dopants and defects. While selecting the right elements as dopants is important for the electronic properties of the material, microstructure, defects and stoichiometry are critical in designing the final properties of the material. There are many examples in which doped perovskites exhibit ordered superstructures or, at least, there are preferential substitutions or configurations [47][48][49]. For instance, Ni atoms preferentially substitute Nb atoms that are close to oxygen vacancies [20].</p><p>Predicting off-stoichiometric material properties is a complex and usually very expensive task. Methodologies such as the virtual crystal approximation [50,51], coherent potential approximation [52] and special quasirandom structures [53] have been developed to tackle this challenge. Alternatively, some authors have proposed the combination of symmetry analysis with large-ordered supercells to perform a rigorous statistical treatment of substitutional disorder [54][55][56].</p><p>Here, we present a high-throughput framework to predict the key properties of ferroelectric doped perovskites in order to accelerate the design of new photovoltaic materials. Properties such as band gap, radial distribution function or polarization are predicted based on a statistical model, using a supercell approach and taking advantage of symmetry. As proof of concept, doped KNbO 3 -based materials were studied obtaining valuable information to optimize their performance.</p><!><p>An automatic framework was developed to explore systematically doped perovskites. Three distinct stages can be differentiated in the framework: i) supercell construction, ii) density functional theory calculations and iii) properties prediction (Figure 1). SOD package [54] was integrated in the process, creating all the inequivalent doped supercells and calculating their degeneracy using as starting point the pure KNbO 3 unit cell. This package has been successfully used before for the study of different off-stoichiometric materials [57][58][59]. SOD package requires the previous analysis of the supercell symmetry including the space group and their symmetry operators. This data was obtained using the spglib symmetry libraries [60] and PHONOPY package [61]. When more than one type of dopant is included in the supercell, an automatic iterative process was designed to use the SOD and PHONOPY packages as many times as required. All the generated supercells are fully relaxed using first principles and then single point calculations are performed to predict their structural, electronic and optical properties. These structures will be available in 2020 in the ioChem-BD platform [62] and NOMAD repository [63]. The average properties of the off-stoichiometric material are calculated using the contribution to the partition function of each supercell, P i ,</p><p>where g i is the degeneracy of the i th supercell, ∆E i is the DFT energy difference between the i th supercell and the ground state, k B is the Boltzmann constant and T is the temperature. This approach should be taken only as an approximation because, i) ∆E i should be substituted by ∆G i which would require the vibrational properties of each supercell and ii) T is only a parameter that describes how much disorder has been explored during the synthesis. Once P i (T ) is calculated for each supercell, any material property, N , for the ensemble can be calculated as,</p><p>Primitive Cell < l a t e x i t s h a 1 _ b a s e 6 4 = " a y X 9 T q r 7 N 0 H G R J J q E H R x q J 9 y r G M 8 C w T 3 m A S q + c Q Q Q q W J g G I 6 J J J Q k 4 q y T Q j e 8 s u r p F k p e 5 f l q / t K q X q b x 1 F E p + g M X S A P X a M q u k N 1 1 E A U P a J n 9 I r e r C f r x X q 3 P h a</p><!><p>< l a t e x i t s h a 1 _ b a s e 6 4 = " s i c 8</p><!><p>< l a t e x i t s h a 1 _ b a s e 6 4 = " 4 V 9 s r 7</p><!><p>All structures were fully relaxed with the VASP code [64][65][66] using the projector-augmented wave, PAW, method [67,68]. Energies were obtained using the generalized gradient approximation, GGA, proposed by Perdew et al. [69] and a plane-wave basis set with a cutoff of 500 eV. DFT+U calculations were performed according to Dudarev formalism, using a single U eff parameter (U eff = U − J) [70]. Pseudopotentials and the Hubbardlike term for each atom were selected following the standards proposed by Calderon et al. [71]. The Brillouin zone was sampled using a k-points mesh of 2500 and 8500 k-points per reciprocal atom for the geometry optimization and single point calculations, respectively. Hellmann-Feyman theorem was used to calculate the forces over the ions, including the Harris-Foulkes corrections [72]. Wavefunction was converged when the energy difference between two consecutive SCF steps was smaller than 10 −6 eV and 10 −8 eV for geometry optimizations and density of states, DOS, single-point calculations, respectively. Geometries were considered optimized when forces acting over each atom were below 10 −3 eV/ A. Radial distribution functions, g(r), and their running integration numbers, are calculated for each fully optimized structure. Polarization, P α , was calculated in the KNbO 3 film as [73][74][75],</p><p>where Ω is the volume of the supercell, N is the number of cations in the supercell and δα i is the relative displacement in the α direction between the cations and the anions, O, located in the same β − γ plane. Z * i,α is the component of the Born effective charges perpendicular to the β−γ plane computed using density-functional perturbation theory [76,77]. This approach only considers the polarization due to atom displacements with respect to their centrosymmetric positions, but it is accurate enough to predict the KNbO 3 bulk P , 35 µC cm −2 , compared to experimental values, 37-41 µC cm −2 [78][79][80]. [11,81]. This solid solution has been widely studied because it fulfils the main requirements to be used in photovoltaic devices, exhibiting ferroelectricity and a customizable band gap in the visible region [11,23] However, the properties of this solid solution are highly dependent on dopant concentration [23,27,82], synthesis method [83][84][85][86] and temperature [25,87,88]. In order to understand how these variables modify the doped system, a model containing 40 atoms (2×2×2 KNbO 3 supercell) was built. This supercell size has been previously used to explore the effect of different type of dopants in KNbO 3 [19]. Two Nb atoms were substituted by Ni atoms and 4 K atoms were substituted by Ba atoms, which corresponds to x = 0.5 for the</p><!><p>x and it is in the upper limit of the dopant content that has been experimentally synthesized for this system [11,23,25]. The model size and dopant concentration were chosen by balancing the computational cost with two important requirements: i) the presence of at least one O vacancy, O V , for charge compensation and ii) the inclusion of at least two atoms of each dopant to understand the dopant distribution. Tetragonal KNbO 3 was used as the primitive cell for building the initial supercell model because it is the non-centrosymmetric polymorph with a range of temperature stability that is reached in most of the synthesis approaches [11,27,83]. Around 1700 possible supercells were found, and almost 970 unique structures were fullyoptimized and studied using first principles. Relative energies, ∆E, are grouped in three clusters (Figure ): one between 0-1 eV, a second one between 1-2.5 eV and the last one between 2.5-4 eV. KNbO 3 undergoes important structural changes when it is doped. Li et al. reported a very detailed phase diagram, describing the most stable KBNNO phases as function of composition and temperature [25]. However, some discrepancies can be found when literature is reviewed. For instance, the maximum content of dopants, x, that stabilizes the cubic phase, varies from 0.1 to 0.5 [23,28,83]. This could be due to the almost identical a and c parameters for the tetragonal phase when x is increased [87], making the typical splitting of signals in XRD difficult to detect [84,89]. A similar trend was found when the supercells were fully optimized. The a/c ratio for the optimized tetragonal primitive cell is 0.95, however in the doped structures, an average a/c ratio of 0.98 was found, confirming the cubic phase stabilization when x is increased. Microstructure or short-range order is also modified because of the dopants. The addition of Ni 2+ cations, with a larger ionic radius (0.69 Å) than Nb 5+ cations (0.64 Å) produces stress, which is partially reduced by the substitution of K + cations (1.38 Å) by Ba 2+ (1.35 Å) cations and the formation of oxygen vacancies [27]. Oxygen vacancies, O V also play an important role in balancing the charge [20]. As it was mentioned before, our models include one O V , which is enough to balance the charge of the system. However, additional calculations were performed to estimate the formation energy of a second O V . Formation energies for this second vacancy were always found to be higher than 3.5 eV, which makes them thermodynamically unstable. Dopant and vacancy distribution have been linked to the electronic and optical properties of the solid solution, hence radial distribution function, g(r), has been explored for each element (Figure 3). For instance, Ni presents a first coordination shell with O atoms at two different distances (Figure 3a). It is also possible to extract the number of atoms that correspond to each peak if combining g(r) with the running coordination number, n(r). The first coordination shell for Ni is constituted of 5 O atoms, one of these O atoms is at a shorter distance and the other 4 O atoms exhibit, in principle, similar distances at low temperatures (Figure 3c). Here again, we would like to highlight that the term temperature is used to describe the level of disorder that has been explored during the synthesis. At low temperatures, only the most stable structure (ground state) or supercells with similar stability contribute to g(r). However, at high temperatures, more structures contribute in the description of the properties of the ensemble. The coordination number for the first shell of O atoms around Ni slightly increases with the temperature. This fact corresponds to a distribution in which the O V is not any more coordinated to both Ni atoms, but one of them is partially substituted by a Nb atom. Moreover, Ni-O distances are enlarged and they are not as homogeneous as they are when observed at low temperatures, which is in agreement with experimental observations in which Ni 2+ ions occupy the center of a distorted octahedra [86]. The opposite trend can be observed for Nb atoms (Figure 3b) which present a coordination number of 6 at low and mid temperatures, which slightly decreases with the temperature. The coordination number for Nb decreases slower than the coordination number for Ni (Figure 3c) because there are three times more Nb atoms than Ni, and only one Nb is close to the vacancy. The shape of the first peak (Figure 3b) also shows that the Nb atoms presented a distorted octahedra coordination, including a shorter Nb-O bond at high temperatures, when the Nb that is located next to the vacancy. The same analysis has been applied to Ba and K atoms (Figure 3c). At low temperature, all 4 K atoms are strongly preferred around the vacancy. As mentioned, Ni 2+ and K + present larger ionic radii than Nb 5+ and Ba 2+ , so we can expect that they will occupy the A and B sites around the vacancy to minimize stress. However, the relative size difference between K + and Ba 2+ is smaller than the Ni 2+ -Nb 5+ pair so the exchange between K + and Ba 2+ around the vacancy is observed as soon as higher temperatures are explored.</p><p>Structural changes in KNbO 3 are paired with important modifications on the optoelectronic properties of the solid solution. These changes not only depend on the dopant concentration [11,23,90] but also on the synthesis. Both sol-gel methods [83] and pulsed laser deposition techniques [84] seem to obtain solid solutions with larger band gaps than samples obtained through traditional solid state methods in which calcination and sintering are common steps [11,27] Interestingly, different electronic and optical properties are reported when using the same synthesis method, but changing the temperature at which the sample is sintered [85,86]. When electronic properties are considered, the most important change is the reduction of the band gap, for any dopant concen-</p><p>f F e X c + l q 0 b T j 5 z B n / g f P 4 A x 3 2 M 7 g = = < / l a t e x i t > c < l a t e x i t s h a 1 _ b a s e 6 4 = " F Q 2 q c 3 l p 4 A + S i F w L y X P p X I r B I s tration or synthesis method. This band gap reduction seems to be related to the substitution of Nb atoms by Ni atoms, which creates occupied mid-gap states above the valence band of the KNbO 3 . In order to analyze the effect of doping and defects on the electronic structure of KNbO 3 , a band unfolding technique to recover an effective primitive cell picture has been applied [91,92]. Using this approach, it is possible to make a direct comparison between the band structure of KNbO 3 (Figure 4a) and 4b). PBE+U predicts an indirect band gap, E g , around 1.75 eV (Figure 4a), however, experimental values are reported between 3.08 eV and 3.30 eV. [93,94] Despite of the differ-ence, PBE+U functional represents a good trade-off between accuracy and computational resources and allow us to examine the trend in the electronic structure when the effect of dopants is considered. The band gap value is reduced to 1.60 eV (Figure 4b) for the most stable doped structure, including important modifications in the band structure. For instance, the first unoccupied band, mainly projected on the Nb 3d xy orbital at the tetragonal phase [17], is pushed to higher energies (Figure 4b) . Simultaneously, a new occupied flat band appears just above the valence band of KNbO 3 , which corresponds to Ni 3d orbitals. The position of this band is the factor that determines the band gap of the material. Different theoretical works have demonstrated that lower band gaps are obtained when Ni atoms are coordinated to 6 O atoms at low dopant concentrations [20,23]. The crystal field splitting of the Ni orbitals is modified when the vacancy is included in its first coordination shell, stabilizing the energy of the occupied orbitals and, at the same time, increasing the band gap [20]. However, at high dopant concentrations, there is a strong preference for locating both Ni atoms around the O V , even at high temperatures (Figure 3c). This strong preference explains how band gap values for</p><p>x first decrease at low dopant concentration, but then are higher when x is increased [11,23,27]. As it was mentioned above, sintering temperature also modifies the optoelectronic properties of the solid solution for a given dopant concentration [85,86] The effect of temperature can be analyzed with the statistical framework proposed in this work.</p><p>First, the density of states, DOS, of the ground state structure is compared with a DOS of the ensemble at 500 K (Figure 4c). In both cases, the mid-gap Ni states can be observed just above the valence band. However, the position of this Ni states is slightly higher in energy than in the ground state, so the band gap is reduced. The influence of the temperature on the band gap is depicted in Figure 4d. Band gap values first tend to decrease below 1.5 eV, and then the band gap for the ground state is recovered at higher temperatures. In order to get a more direct comparison with the experiments, DOS were also calculated with the HSE functional [95,96] for 20 supercells which have been selected randomly. There is a direct linear correlation between the band gap values obtained with PBE+U and HSE functionals (Figure S1) This correlation has been used to extrapolate the band gap values obtained for the ensemble with the PBE+U functional to HSE band gap values at different temperatures. Using HSE, band gap values are located between 2.5 eV and 2.7 eV. Experimental band gap range for x = 0.5 (1.7-1.9 eV [11,23,25]) is between the values obtained for PBE+U and HSE functionals.</p><p>At this point, most of the attention has been focused on Ni atoms, however, Ba atoms could play an important role, too. Previous studies reported that dopants at A sites modify the valence band and can change the directto-indirect transition of the band gap nature of the band gap. Recent experimental reports, in which KNbO 3 has been doped by Ni but not with Ba, have obtained band gap values higher than 3.0 eV [90] which could mean that Ba atoms play also an active role reducing the band gap. Structurally, we already found that Ba and K ions can be exchanged more easily around O V than Ni and Nb when temperature is increased. In order to analyze the effect of Ba atom positions around the O V , different supercells were built in which the positions of the Ni and Nb atoms around the vacancy remain as in the ground state, but the number of Ba and K around the vacancy are modified (Figure 4e). Band gap values clearly decrease when the number of Ba ions around the vacancy increases, so Ni 3d orbital energies are not only modified by the presence of a vacancy, but also by the environment around that vacancy. This fact could explain the band gap reduction of the ensamble at the mid-temperature range. At these temperatures, the substitution of one or two K ions by Ba ions around the vacancy is observed in the radial distribution function and running coordination numbers.</p><p>Preserving polarization is as important as reducing the band gap of KNbO 3 in order to keep the ferroelectric properties of the material. Doping usually reduces the polarization of the system [19], so prediciting the polarization of the solid solution is essential to evaluate the performance of these materials. Here, we present the calculated polarization for the ensemble for different temperatures (Figure 4f). Two main conclusions can be extracted from Figure 4f. First, polarization substantially decrease if we compare these values with KNbO 3 bulk. Second, the polarization is reduced when temperature is increased. Both facts are in agreement with experimental measurements, obtaining values in the range of 10-20 µC cm −2 , which are reduced when temperature is increased [11] IV. CONCLUSIONS Solids solutions are relevant systems for a wide variety of technologies. However, they are complex materials whose properties depend on many variables such as dopant concentration, synthesis method or temperature. This complexity is even greater when more than one property has to be monitored. This is the case of doped ferroelectric materials where band gap and polarization need to be optimized. Here, we have presented a framework to study doped ferroelectric materials combining Boltzmann statistics and DFT calculations. This approach presents two main advantages: i) the automatic and systematic exploration of the configurational space and ii) the connection between the changes in the microstructure of the material and its photovoltaic performance.</p><p>One of the most studied doped-ferroelectric systems, [KNbO 3 ] 1−x [BaNi 1/2 Nb 1/2 O 3−δ ] x , has been used as a proof of concept. More than 900 inequivalent supercells were optimized and their electronic properties calculated. Radial distribution function and running coordination number shows a strong preference of Ni atoms to occupy the positions around the oxygen vacancy, even at high temperatures. However, K atoms are exchanged for Ba atoms around the vacancy when temperature is increased. These micro-structural changes are connected with the modifications observed on the electronic structure. Ni 3d occupied orbitals generate mid-gap states that reduce the band gap. The strong preference of Ni atoms for positions next to the O mathrmV at high dopant concentrations stabilizes the mid-gap states, increasing the band gap with respect to intermediate dopant concentrations. The number of Ba atoms around the vacancy also modifies the position of the mid-gap states, modifying the band gap. Doping reduces the band gap, but also reduces the polarization of the solid solution. The predicted values are in the same range than previous experimental reports.</p><p>The presented framework has demonstrated that it can bridge the microstructure of solids solutions with their optoelectronic properties and, most importantly, connecting them with some of the variables that can be modified during the synthesis process. We believe that this framework could be used as an efficient tool to strengthen the connection between theoretical and experimental results, thus accelerating the discovery, design, and optimization of new doped photosensitized ferroelectric materials.</p><!><p>There are no conflicts to declare.</p>