Origin of Ferroelectricity in BiFeO3-Based Solid Solutions

We investigate the origin of ferroelectricity in the BiFeO3–LaFeO3 system in rhombohedral R3c and tetragonal P4mm symmetries by ab initio density functional theory calculations and compare their electronic features with paraelectric orthorhombic Pnma symmetry. We show that a coherent accommodation of stereo-active lone pair electrons of Bi is the detrimental factor of ferroelectricity. A Bloch function arising from an indirect Bi_6p–Fe_3d hybridization mediated through O_2p is the primary origin of spontaneous polarization (Ps) in the rhombohedral system. In the orthorhombic system, a similar Bloch function was found, whereas a staggered accommodation of stereo-active lone pair electrons of Bi exclusively results in paraelectricity. A giant Ps reported in the tetragonal system originates from an orbital hybridization of Bi_6p and O_2p, where Fe-3d plays a minor role. The Ps in the rhombohedral system decreases with increasing La content, while that in the tetragonal system displays a discontinuous drop at a certain La content. We discuss the electronic factors affecting the Ps evolutions with La content.

Karimi et al. [28] reported that a phase boundary between the ferroelectric R3c and the paraelectric orthorhombic Pnma phases exists at x~0.23 at room temperature. Rusakov et al. reported that single-phase materials with R3c symmetry can be prepared after annealing for composition 0 ≤ x ≤ 0.1, and the Pnma phase is stable at 0.50 ≤ x ≤ 1; these results were verified by synchrotron radiation X-ray diffraction, electron diffraction, and high-resolution transmission electron microscopy [38]. They also found that an incommensurate phase in orthorhombic Imma symmetry is formed at 0.19 ≤ x ≤ 0.30. Karpinsky et al. [29] proposed a temperature-composition phase diagram, in which the ferroelectric R3c (x < 0.15) and the paraelectric Pnma (x > 0.4) phases are mediated through a bridging anti-polar phase in orthorhombic Pbam symmetry.
One possible reason for the wide variety of experimental reports on the phase diagrams is incomplete solubility of La on the A site [39]. Moreover, the abovementioned phases are energetically competing with each other and some of them are likely to be energetically degenerate [31]. Therefore, ab initio studies based on density functional theory (DFT) are expected to provide clues for uncovering the ground-state crystal structure and the phase stability.
Lee et al. investigated the effect of the La doping on the variation of the off-center distortion and the orbital mixing in BiFeO 3 by experiments in conjunction with DFT calculations [30]. They reported that both an Fe-O bond anisotropy and off-center cation displacements are suppressed by the La doping. As a result, the degree of Fe 3d-4p orbital mixing decreases in the solid solution samples. An impact of the La content on the polarization and the electronic band structure was also reported by You et al. [40]. They reported that the La doping induces a chemically driven rotational instability. It modifies the local crystal field along with the electronic structure, which gives rise to a direct-to-indirect transition of the bandgap and provides an enhancement in ferroelectric photovoltaic effect. In contrast, Tan et al. reported that the La doping has little influence on P s in tetragonal BiFeO 3 [34]. In spite of extensive research by DFT studies [30][31][32][33][34][35], the P s 's evolution with the La content and its electronic origin still remain unclear.
The purpose of this study is to elucidate the origin of ferroelectricity in rhombohedral R3c, tetragonal P4mm in the BiFeO 3 -LaFeO 3 system (Bi 1−x La x FeO 3 ). The electronic feature and structural distortion in the paraelectric orthorhombic Pnma symmetry are also investigated because Bi 1−x La x FeO 3 with x ≥ 0.5 is of the orthorhombic phase [29]. We show DFT energy evolutions with x but focus on the relation between the orbital hybridizations and the ferroelectric (paraelectric) distortions. We show that the Bloch function arising from a Bi_6p-Fe_3d hybridization mediated through O_2p is the primary origin of P s in the rhombohedral system. In the orthorhombic system, a similar Bloch function and the resultant structural distortion are constructed, whereas a staggered accommodation of stereo-active lone pair electrons of Bi never allow the presence of P s . We discuss a large P s and its dependence on x in the tetragonal system.

DFT Calculations
Density functional theory (DFT) calculations were conducted using the generalized gradient approximation [41] with a plane-wave basis set. We used the projector-augmented wave method [42] as implemented in the Vienna ab initio simulation package (VASP) [43]. We employed the Perdew-Burke-Ernzerhof gradient-corrected exchange-correlation functional revised for solids (PBEsol) [44], a plane-wave cut-off energy of 520 eV, an electronic iterations convergence of 1 × 10 −6 eV, and a criterion for ionic relaxations of 1 meV/nm. The Γ-centered k-point mesh was set to 3 × 3 × 3 for the structural optimizations and 5 × 5 × 5 for density of states (DOS) and band structure calculations.
Within the simplified generalized gradient approximation (GGA) + U approach [45], we added on-site Coulomb interaction parameters of U-J of 6 eV to Fe-3d throughout the calculations. The on-site Coulomb interaction parameters of U-J for Fe-3d has been employed in the range of 2-6 eV for BiFeO 3 [34,[45][46][47]. The bandgap value is enlarged when U-J is increased, while the essential feature, such as P s , and the valence-band electronic structure remain unchanged. One main reason why we adopted U-J = 6 eV is as follows: the bandgap becomes narrow for a specific Bi-La arrangement on the A site and eventually vanishes when the arrangement of Bi and La is an ordered configuration along the polar c axis, as will be described later. In order to maintain the bandgap above ca. 2 eV, we set U-J to 6 eV for Fe-3d throughout the calculations.
Considering the spin configuration in BiFeO 3 can be approximated as the G-type antiferromagnet [46], we set the spin arrangement in which the adjacent Fe ions have an antiparallel spin configuration as much as possible irrespective of the La content (x) on the A site.
We employed three symmetries shown in Figure 1 and compared their total energies (see Figure 2): ferroelectric rhombohedral R3c, ferroelectric tetragonal P4mm, and nonpolar orthorhombic Pnma. The antiparallel spin configuration of adjacent Fe ions enforces a change in space group from R3c to R3 (or P3) for the rhombohedral cells, from Pnma to P2 1 /m for the orthorhombic ones, while the tetragonal P4mm remains unchanged, the details of which are summarized in Tables S1-S3. Although the orthorhombic cells were optimized in monoclinic P2 1 /m symmetry, we restricted the monoclinic angle of β to 90.0 degree, and then the crystal system is regarded as orthorhombic.  Because our calculations were performed in the anti-ferromagnetic spin configuration, the total DOS is identical both in the majority and minority spin bands. In a similar manner, the band structure in the majority spin band is the same with that in the minority spin band. We display the DOS in the majority spin band in the right panel and that in the minority spin band in the left one, e.g., see Figure 3c. For simplicity, we show the band structure only in the majority spin band. As for wavefunctions (partial charge densities), we visualize the majority spin component at an iso-surface level of 1 × 10 −4 , e.g., see Figure 4. The Fermi energy is set to zero throughout our manuscript.  For building a BiFeO 3 -LaFeO 3 solid solution cell, there exists several choices of the arrangement of Bi and La. We adopted the rock-salt arrangement of Bi and La, especially along the polar c axis, as much as possible, as can be seen in Figures S1 and S2. This is because the electronic structure in the rhombohedral cell with a layered Bi and La ordering along the polar c axis has a non-realistic metallic feature ( Figure S10), which is not consistent with the experimental fact of an insulating nature of BiFeO 3 -LaFeO 3 solid solutions [47,48]. We therefore avoid such a Bi-La ordering along the specific crystallographic axis and adopt the rock-salt-like orderings of Bi and La. The details of the structure parameters are listed in Tables S4-S10 for the rhombohedral cells, Tables S11-S18 for the tetragonal cells, and Tables S19-S21 for the orthorhombic cells.
From the structure parameters of the optimized cell (Tables S4-S18), we obtained the atomic displacements (∆z) from the corresponding positions in the hypothetical non-polar paraelectric lattice. We also calculated the Born effective charges (Z*) [49] by densityfunctional perturbation theory. We estimated spontaneous polarization (P s ) using the following equation where m i denotes the site multiplicity of the constituent atom i and ∆z i ·Z * i its dipole moment. The summation in Equation (1) is taken over the cell with cell volume (V). . For comparing E total in different symmetries, x should be the same. Indeed, x has a discrete value depending on the total number of the A-site atoms, i.e., six in the rhombohedral cells, eight in the tetragonal ones, and four in the orthorhombic ones (Tables S1-S3). We, therefore, set the energies obtained by fitting the DFT energies of the rhombohedral cells by a quadratic function to zero of E total . The E total of the rhombohedral cells are lower than those of the tetragonal and orthorhombic cells irrespective of x. These results are different from those reported in the literature [31,50,51]. The orthorhombic cells are slightly higher in E total by~0.1 eV while the tetragonal cells have a higher E total by~0.3 eV. In reality, the tetragonal BiFeO 3 is stabilized under a compressive strain in epitaxial films [26]. Figure 2b shows the P s of the rhombohedral cells as a function of x. The P s of the rhombohedral BiFeO 3 is 83.2 µC/cm 2 , which is in good agreement with the previous DFT calculations [26,31,52] and experiments [53]. With increasing x, the rhombohedral P s exhibits a monotonic decrease, which is consistent with other calculations [31]. The tetragonal BiFeO 3 (Figure 2c) possesses a giant P s of 133.4 µC/cm 2 , which is close to the values obtained for epitaxially strained films and from DFT calculations [26]. Zhang et al. [26] have reported that this large P s is associated with a high tetragonality c/a of 1.24 and also with a coherent displacement of Fe (∆ Fe ') by 0.033 nm with respect to the Bi sublattice. These values accord with our calculations of c/a = 1.25 and ∆ Fe ' = 0.032. The P s and c/a decrease with increasing x followed by discontinuous drops at 2/8 < x < 3/8. At x = 3/8, the c/a is almost in unity whereas the apparent P s of 52.0 µC/cm 2 exists. With a further increase in x, the P s again shows a monotonic decrease while the c/a remains constant (ca.~1.0).  Figure 4. The off-center displacements of Fe and Bi are stabilized by the following two orbital hybridizations, respectively: Fe_3d-O_2p (Figure 5a) and Bi_6p-O_2p (Figure 5b). Because Fe atoms have either a positive (↑) or a negative (↓) magnetic moment, the majority spin (↑) and minority spin (↓) bands have to be taken into account in a distinct manner.

Electronic Structures
Here, we consider the orbital interaction in the ↑ band that leads to DOS components in the valence band in the range of −7 to 0 eV (in the right panel in Figure 3c-g). Although O1 and O2 have a slightly different magnetic moment, the DOS characters are almost identical both in the ↑ and ↓ bands. When the magnetic moments are ignored, O1 and O2 have the same site symmetry, and therefore, we do not distinguish O1 and O2. Due to the same reason, we regard Bi1 and Bi2 as identical Bi atoms. Next, we consider the Fe_3d-O_2p interaction leading to the DOS component in the ↑ band. The Fe1 atom with a negative magnetic moment has the electron configuration of 3d with ↓ 5 ↑ 0 . These five ↓ electrons of Fe1 are present at deep levels in the ↓ band [54]. We focus our attention on the interaction between the empty Fe1_3d (↑) and the occupied O_2p states (Figure 5a). The hybridization of these orbitals delivers an occupied bonding state and an empty antibonding state. Therefore, Fe1 has not only the major DOS in the band but also the minor DOS in the ↑ band, as displayed in Figure 3d.
The bands in the range of −7 to −5 eV (Figure 3h) are derived primarily from the bonding states of Fe2_3d (↑)-O_2p. The Bi_6p-O_2p interaction yields a low-lying bonding state at around −4.5 eV at the Γ point (marked with a blue circle in Figure 3i), whose wavefunction is visualized in Figure 4. The distinct Fe1_3d-O_2p and Bi_6p-O_2p interactions are seen in the respective local regions of Fe1O 6 octahedra and BiO polyhedrons. We note that the interaction between the bonding states of Fe1_3d-O_2p and Bi_6p-O_2p ( Figure 5c) forms a coherent wavefunction that is spread throughout the crystal (Figure 4). Namely, the interaction between the Fe1_3d-O_2p-and the Bi_6p-O_2p-derived bonding states yields the low-lying bonding state, termed Bloch function, arising from the -Fe1-O-Bi-O-network ( Figure 4). It results also in the high-lying antibonding state in the valence band, as shown in Figure 5c. We conclude that the Bloch function stemming from the indirect Bi_6p-Fe_3d hybridization mediated through O_2p is the primary origin of P s in the rhombohedral system. Figure S3 shows (a) the crystal structure of the rhombohedral cell (x = 2/6) along with (b-f) the total and partial DOS results and (g,h) the electronic band structures in the valence band. The wavefunction of the band shown in the blue circle in h is depicted in i. The orbital interactions displayed in Figure 4 are seen also in this cell. The Bloch function arising from the Fe_3d-O_2p-Bi_6p hybridization appears at −4.5 eV (at the Γ point), while its connection along the c axis is relatively weak compared with that in the BiFeO 3 cell. The similar electronic feature was found for the rhombohedral cell (x = 4/6) in Figure S4. It is interesting to note that the Bloch function is formed through the Bi-O-Fe bond avoiding La.

Ferroelectric Tetragonal System
The large P s of the tetragonal BiFeO 3 is derived from the cooperative off-center displacements of Bi and Fe along the polar c axis with respect to the oxygen sublattice ( Figure S5 for the BiFeO 3 cell). This off-center feature is maintained in the tetragonal cell with x = 2/8. Figure S6 shows (a) the crystal structure of the tetragonal cell (x = 2/8) along with (b-f) the total and partial DOS results and (g, h) the electronic band structures in the valence band. The wavefunction of the band shown in the blue circle in h is depicted in i. The Bloch function where a Bi_6p-O2_p interaction plays a central role ( Figure S6i) is present at ca. −4.2 eV at the Γ point and spreads also along the c axis through the Bi_6p z orbital. Note that the Fe-3d state indeed does not participate in the Bloch function, which is in contrast to the rhombohedral system (Figure 4, Figures S3 and S4). Additionally, in the BiFeO 3 cell (Figure S5), the similar Bloch function with a small contribution of Fe-3d appears. Figure 6 shows (a) the crystal structure of the tetragonal cell (x = 3/8) along with (b-f) the total and partial DOS results and (g,h) the electronic band structures in the valence band. The wavefunction of the band shown in the blue circle in h is depicted in (i). The Bi_6p-O2_p hybridization results in a bonding state at ca. −4.6 eV at the Γ point. The Bloch function has a robust connection in the a-a plane whereas that is markedly weakened along the polar c axis. The discontinuous drops in P s and c/a at 2/8 < x < 3/8 stems from the in-plane feature of the Bloch function formed by the Bi_6p-O2_p hybridization.

Paraelectric Orthorhombic System
We think that the variations of P s and c/a can be qualitatively understood from a decrease in the number of the Bi pillar along the c axis. Figure S7 shows the arrangement of Bi and La in the tetragonal cells at (a) x = 2/8 and (b) x = 3/8. The tetragonal BiFeO 3 (x = 0) has a large P s , which is ascribed to the full set of the Bi pillars. The substitution of La on the A site decreases the number of the Bi pillar and two Bi pillars are maintained until x ≤ 2/8. These Bi pillars contribute to the formation of the Bloch function spreading along the c axis through the Bi_6p z orbital. With increasing x above 3/8, the number of the pillar becomes only one, and thereby the Bloch function has an in-plane feature, which is accompanied by a marked decrease in P s . Figure 7 shows (a) the crystal structure of the orthorhombic cell (x = 1/2) along with (b-e) the total and partial DOS results and (f,g) the electronic band structures in the valence band. The wavefunction of the band shown in the blue circle (g) is depicted in (h). The magnetic moment of Fe1 is negative, and then the orbital interaction shown in Figure 6a is expected for Fe1 and the adjacent oxygen atoms. The Fe1_3d states have a marked DOS component in the ↑ band, which stems from the orbital hybridization of Fe1_3d (↑)-O3_2p. An apparent DOS of Bi1_6p arising from the mixing with O3_2p appears in the valence band. The wavefunction of the band at the Γ point (Figure 7h) clearly shows that the Bloch function originating from a coherent interaction between the bonding states of Fe1_3d (↑)-O3_2p and Bi1_6p-O3_2p spreads throughout the crystal. Moreover, this Bloch function is present at −5.4 eV, which is lower by 0.9 eV than that of the rhombohedral cell (x = 2/6) ( Figure 3).

Factors Affecting Ferroelectricity
As described above, the Bloch function arising from the indirect Bi_6p-Fe_3d hybridization via O_2p is the origin of P s in the rhombohedral system. The hybridized orbital is accompanied by the formation of the covalent bonds not only in the direction of P s but also in the opposite direction of P s , e.g., see Figure 4. The Bloch function is closely related to the bond lengths of Bi-Fe (see Tables S1-S3); this is the first (structural) factor leading to a robust P s . In the rhombohedral system ( Figure S9a), the first shortest Bi-Fe bond is parallel to P s , while the second shortest Bi-Fe bond is almost normal to P s . For the tetragonal system ( Figure S9b), the first and second shortest Bi-Fe bonds are aligned along pseudo-cubic <111> direction. In the rhombohedral cells at x ≤ 4/6, the first shortest Bi-Fe length is 0.31 nm, and the second shortest one is~0.33 nm. In the tetragonal cells, the first shortest Bi-Fe length is~0.33 nm, which is comparable to that in the rhombohedral one. We note that the second shortest length is as long as~0.37 nm even at x = 0 (BiFeO 3 ). This long Bi-Fe length does not allow the Fe_3d orbital to participate in the Bi_6p-O_2p hybridization, and thereby, the Bloch function is indeed formed by the Bi_6p-O_2p interaction in the tetragonal system, as shown in Figure 6.
In the paraelectric orthorhombic cells at x ≤ 1/2, the first and second shortest Bi-Fe lengths are~0.32 nm and~0.33 nm, respectively, which are comparable to those in the rhombohedral cells. Actually, the Bloch function resulting from the Fe_3d-O_2p-Bi_6p mixture is extended throughout the crystal at x = 1/2 (see Figure 7). The most significant difference between the ferroelectric and paraelectric cells appears in an accommodation mode of stereo-active lone pair electrons of Bi [26,55,56] (derived from the mixture of Bi_6s-6p), which is the second factor that dominates the presence or absence of P s . For the ferroelectrics in the rhombohedral (Figure 4a) and tetragonal systems (Figure 6i), the lone pair electrons of Bi are directed coherently along the c axis, which is the detrimental factor of ferroelectricity. The alignment of the lone pair electrons contributes to an enhancement of off-center displacements of cations, as in ferroelectric PbTiO 3 [57]. For the orthorhombic cell (Figure 7h), the lone pair electrons of Bi1 (y = 0.75) are directed opposite to those of Bi1 (y = 0.25) along the c axis (which is the symmetry constraint). This staggered accommodation of the lone pair electrons of Bi exclusively provides paraelectricity.
The Bi-Fe bond lengths and the resultant Bloch functions can be qualitatively understood by the tilt and rotation of FeO 6 octahedra expressed by the Glazer notation [58]. The ferroelectric rhombohedral cells have an out-of-phase octahedral tilt expressed by a − a − a − , which are accompanied by the short Bi-Fe bonds along with relatively large unit-cell densities of 8.595 g cm −3 at x = 0 and 7.779 g cm −3 at x = 0.5. In contrast, the ferroelectric tetragonal cells do not have any tilt or rotation of FeO 6 octahedra at x = 0 and 0.5, leading to the long (second shortest) Bi-Fe lengths and low densities of 8.027 g cm −3 at x = 0 and 7.669 g cm −3 at x = 0.5. The paraelectric orthorhombic cells have a tilt system of a − b + a − , resulting in the short Bi-Fe lengths and high densities of 8.836 g cm −3 at x = 0 and 7.807 g cm −3 at x = 0.5.

Conclusions
We have investigated the origin of ferroelectricity of Bi 1−x La x FeO 3 in rhombohedral R3c and tetragonal P4mm symmetries by DFT calculations. In the rhombohedral system, a Bloch function arising from an indirect Bi_6p-Fe_3d hybridization via O_2p is the primary origin of P s . In contrast, the P s of the tetragonal phase stems from a Bloch function arising from a Bi_6p-O_2p mixing with a weak contribution of Fe-3d. The detrimental factor of the presence/absence of P s is an accommodation of stereo-active lone pair electrons of Bi. The paraelectric orthorhombic Pnma phase has a staggered accommodation of lone pair electrons of Bi, while the ferroelectric R3c and P4mm systems exhibit a coherent alignment of lone pair electrons of Bi. The rhombohedral system shows a monotonic decrease in P s with increasing x, which is directly associated with a weakening of the Fe_3d-O_2p-Bi_6p hybridization. In contrast, the tetragonal system displays a discontinuous drop of P s at ca. x = 0.3, which is ascribed to a transition from a 3D extension to an in-plane feature of the Bi_6p-O_2p mixed orbital.
Supplementary Materials: The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/nano12234163/s1, Figure S1: Bi and La arrangement on the A site in the rhombohedral system with Bi 1−x La x FeO 3 ; Figure S2: Bi and La arrangement on the A site in the tetragonal system with Bi 1−x La x FeO 3 ; Figure S3: Crystal structures (a), electronic density of states (DOS) (b-f), and band structures (of the majority spin band) in the valence band (g,h) of the rhombohedral cell with x = 2/6. The wavefunction of the band shown in the blue circle in h is displayed in i; Figure Table S1: Empirical formula, space group, numbers of Bi and La, and bond lengths of Bi-Fe in the unit cells of the rhombohedral cells; Table S2: Empirical formula, space group, numbers of Bi and La, and bond lengths of Bi-Fe in the unit cells of the tetragonal cells; Table S3: Empirical formula, space group, numbers of Bi and La, and bond lengths of Bi-Fe in the unit cells of the orthorhombic cells; Table S4: Structural parameters of the rhombohedral BiFeO 3 cell (space group R3) with an antiferromagnetic spin configuration; Table S5: Structural parameters of the rhombohedral (Bi 5 La)Fe 6 O 18 cell (space group P3) with an antiferromagnetic spin configuration; Table S6: Structural parameters of the rhombohedral (Bi 2 La)Fe 3 O 9 cell (space group P3) with an antiferromagnetic spin configuration; Table S7: Structural parameters of the rhombohedral (BiLa)Fe 2 O 6 cell (space group R3) with an antiferromagnetic spin configurations; Table S8: Structural parameters of the rhombohedral (BiLa 2 )Fe 3 O 9 cell (space group P3) with an antiferromagnetic spin configuration; Table S9: Structural parameters of the rhombohedral (BiLa 5 )Fe 6 O 18 cell (space group P3) with an antiferromagnetic spin configuration; Table S10: Structural parameters of the rhombohedral LaFeO 3 cell (space group R3) with an antiferromagnetic spin configuration; Table S11: Structural parameters of the tetragonal BiFeO 3 cell (space group P4mm) with an antiferromagnetic spin configuration; Table S12: Structural parameters of the tetragonal (Bi 7 La)Fe 8 O 24 cell (space group P4mm) with an antiferromagnetic spin configuration; Table S13: Structural parameters of the tetragonal (Bi 3 La)Fe 4 O 12 cell (space group P4mm) with an antiferromagnetic spin configuration; Table S14: Structural parameters of the tetragonal (Bi 5 La 3 )Fe 8 O 24 cell (space group P4mm) with an antiferromagnetic spin configuration; Table S15: Structural parameters of the tetragonal (BiLa)Fe 2 O 6 cell (space group P4mm) with an antiferromagnetic spin configuration; Table S16: Structural parameters of the tetragonal (Bi 3 La 5 )Fe 8 O 24 cell (space group P4mm) with an antiferromagnetic spin configuration; Table S17: Structural parameters of the tetragonal (BiLa 7 )Fe 8 O 24 cell (space group P4mm) with an antiferromagnetic spin configuration; Table S18: Structural parameters of the tetragonal LaFeO 3 cell (space group P4mm) with an antiferromagnetic spin configuration; Table S19: Structural parameters of the orthorhombic BiFeO 3 cell (space group P2 1 /m) with an antiferromagnetic spin configuration; Table S20: Structural parameters of the orthorhombic (BiLa)Fe 2 O 6 cell (space group P2 1 /m) with an antiferromagnetic spin configuration; Table S21: Structural parameters of the orthorhombic LaFeO 3 cell (space group P2 1 /m) with an antiferromagnetic spin configuration.

Data Availability Statement:
The data that support the findings of this study are available upon reasonable request from the corresponding author.