Surface reconstructions and related local properties of a BiFeO3 thin film

Coupling between lattice and order parameters, such as polarization in ferroelectrics and/or polarity in polar structures, has a strong impact on surface relaxation and reconstruction. However, up to now, surface structures that involve the termination of both matrix polarization and polar atomic planes have received little attention, particularly on the atomic scale. Here, we study surface structures on a BiFeO3 thin film using atomic-resolution scanning transmission electron microscopy and spectroscopy. Two types of surface structure are found, depending on the polarization of the underlying ferroelectric domain. On domains that have an upward polarization component, a layer with an Aurivillius-Bi2O2-like structural unit is observed. Dramatic changes in local properties are measured directly below the surface layer. On domains that have a downward polarization component, no reconstructions are visible. Calculations based on ab initio density functional theory reproduce the results and are used to interpret the formation of the surface structures.

The surfaces of perovskite oxides undergo surface relaxation 1,2 and reconstruction [3][4][5][6][7][8][9] owing to the breaking of translational symmetry. The resulting surface structures can have a strong influence on the functional properties of the materials when they possess large surface-to-volume ratio 4,10 . Surface structures also govern the performance of devices that rely on interfacial coupling or interactions [11][12][13][14] . Depending on the charges of the crystallographic termination planes, perovskite surfaces can be distinguished into two groups: polar surfaces with a net charge and non-polar surfaces without a net charge. The atomic, electronic and magnetic properties of such surfaces can differ significantly from one group to the other 1,6-9,15-18 . In ferroelectric oxides, surfaces can also be charged or uncharged depending on the orientation of the surface plane with respect to the spontaneous polarization P s , resulting in a strong influence on surface structure and properties [19][20][21] . For instance, surface chemistry and surface adsorption/desorption behavior in important ferroelectrics such as BaTiO 3 (BTO) and Pb(Zr,Ti)O 3 (PZT) are highly dependent on P s 19,21 . Although BTO and PZT are polarized along the tetragonal [001] axis at room temperature, it should be noted that the (001) surfaces of these structures, terminated by either a BaO or a TiO 2 plane, are charge neutral.
There is another category of ferroelectrics that possesses ferroelectric polarization and polar atomic planes simultaneously but is little studied. A prototypical material is BiFeO 3 (BFO), which exhibits both ferroelectric and G-type antiferromagnetic order at room temperature 22 . Although the multiferroic phase of BFO has a rhombohedrally-distorted R3c structure 22 , for convenience the crystallographic notation for a pseudocubic unit cell is used throughout the text unless specifically defined otherwise. As illustrated schematically in Fig. 1(a), the (BiO) + and (FeO 2 ) − layers in BFO stack alternately along the < 001> axis, leading to the formation of a polar {001} surface. In addition, off-center displacements of the Fe and O atoms with respect to the Bi sub-lattice result in a large P s of approximately 0.9-1.0 C/m 2 22 along the direction of the tensile-distorted [111] body diagonal. The polarization contains an < 001> component that interacts with the charges on the terminating surface, as illustrated in Fig. 1. Depending on the polarization direction and the termination of the atomic planes, four different surface configurations can be obtained. Figure 1(b) shows all of these surface configurations, which are defined as type I to type IV, respectively. For type I and type IV surfaces, the charges that are caused by the [001] Scientific RepoRts | 7:39698 | DOI: 10.1038/srep39698 component of P s (i.e., P z ) and that of the surface polarity accumulate additively, while for type II and type III surfaces the charges are subtractive and can compensate. Interaction between the bulk polarization and the plane polarity is expected to have a significant influence on surface relaxation and reconstruction. The study of such surface phenomena is of great importance to gain a basic understanding of physical interactions between lattices, charges and polarizations at surfaces or interfaces in BFO, which may involve electrical screening and ferroelectric domain ordering 23 . These interactions will affect the local electronic, ferroelectric and magnetic properties of the material, with implications for potential applications such as exchange bias 24 in interface-controlled devices for future nanoelectronics.
Here, we present an atomic-scale study of the surfaces of a BFO thin film grown on a DyScO 3 (DSO) (110) o substrate, where the subscript o refers to an orthorhombic lattice. By performing atomic-resolution scanning transmission electron microscopy (STEM) and spectroscopy, combined with first-principles density functional theory (DFT) calculations, we reveal two distinct surface structures, which depend on the polarization direction of domains in the film. By relying on an excellent match between the results of ab initio DFT calculations and our experiments, we determine the local ferroelectric and magnetic properties on the basis of atomic positions in the DFT-calculated model. The formation of the two types of surface structure is also discussed on the basis of polarization and the polarity of the surface plane. Figure 2(a) shows a high-angle annular dark-field (HAADF) STEM 25 image of a BFO (001) film on a DSO substrate, recorded along a < 100> direction. In this image, the bright dots (i.e., peaks in intensity) correspond to heavy Bi atomic columns, while the less bright dots correspond to FeO columns. A signal from the pure O columns cannot be distinguished from the background. Along the film normal, the (BiO) + and (FeO 2 ) − planes stack alternately, as in the perovskite structure ( Fig. 1(a)). The chosen sample area contains two domains. A domain wall (DW), as marked by yellow dashed lines, can be traced by following a reversal in the shifts of FeO columns with respect to Bi columns. Figure 2(b) illustrates this reversal more clearly in the form of intensity profiles of the Bi and FeO columns, following the red lines in the left domain and the blue lines in the right domain in Fig. 2(a). The polarization vectors P s in the two domains are denoted by arrows in Fig. 2(a).

Structure and chemistry.
Corresponding to the two domains, the two types of surface structure are recognized in Fig. 2(a), which are separated by the DW. On the surface of the right domain, a double-atomic-layer (DL) is clearly visible, exhibiting image contrast that is similar to that of the BiO planes in the film matrix (see also an intensity plot in supplementary Figure S1). The two atomic layers have a lateral displacement of a < 010> /2 with respect to one another. According to the direction of P s , this surface corresponds to a type I surface, as defined in Fig. 1(b).
The chemistry of the DL was revealed by using atomic-resolution energy-dispersive x-ray spectroscopy (EDXS) elemental mapping, as shown in Fig. 3(a-c). Figure 3(a) shows a magnified HAADF STEM image of the DL, where chemical mapping was performed. Figures 3(b,c) show Bi and Fe elemental maps, respectively, from which it is evident that the DL contains Bi atoms, while no Fe atoms are detected. The Bi DL has a zigzag configuration along the < 100> direction (see dotted line in Fig. 3(a)), while it has a square pattern along the < 110> direction (see supplementary Figure S2). The spacing between the layers was measured to be ~0.26 nm, which is consistent with that of Aurivillius-type layers formed in 0.95(Na 0.5 Bi 0.5 )TiO 3 − 0.05BaTiO 3 thin films 26 .
Structurally, an Aurivillius-type (Bi 2 O 2 ) 2+ layer contains an O atomic plane sandwiched between two Bi atomic planes (see supplementary Figure S3). Annular bright-field (ABF) STEM imaging 27,28 was used to confirm its presence. As shown in Fig. 3(d), the ABF image on the left and its lateral average on the right reveal atomic columns as darker dots on a brighter background, including a signature from the expected pure O columns. A layer with relatively weak contrast (marked by a solid light blue arrow) can be attributed to the presence of an O atomic layer 29 . In this layer, an O atom (light blue) and a neighboring Bi atom (green) form a dumbbell-like configuration, as marked by dotted ellipses, which is reproduced by the simulated ABF image shown in Fig. 3(e). The simulation was calculated on the basis of an atomic O8 model obtained from first-principles DFT calculations (see Methods and supplementary Figures S4-S7). By using the same model, the HAADF image was also simulated and shown to have a good correspondence to the experimental image, as presented in the inset to Fig. 3(a). Based on these results, it is concluded that the DL structure is very close to that of the (Bi 2 O 2 ) 2+ unit in the Aurivillius phase.
The same analysis was applied to the surface structures on the left domain in Fig. 2(a). Figure 3(f) shows a magnified HAADF STEM image of the surface area. Although the intensities of the atomic columns decay slightly at the surface due to a reduction in specimen thickness, the atomic features remain well-resolved. Based on the atomic-resolution HAADF STEM image and EDXS maps shown in Fig. 3(g,h), the terminating plane on the surface is most likely to be a perovskite (BiO) + plane. The (Bi 2 O 2 ) 2+ -like structure that was observed on the surface of the domain with an upward P s component (i.e., a type I surface) is now absent. According to the polarization direction (i.e., a downward component of P s ), the surface of the domain can be classified as type II, as defined in Fig. 1 Lattice expansion. The lattice parameters in the surface regions were investigated by quantifying the atomic-resolution HAADF STEM images. The positions of intensity peaks corresponding to Bi atomic columns were determined from fast-acquisition HAADF STEM images by fitting two-dimensional Gaussian functions to the intensity peaks 30 . Based on the measured positions of the intensity peaks, the tetragonality (i.e., the c/a ratio) was calculated for each unit cell. Figure 4(a) shows the mean values of c/a and the corresponding standard deviations σ as error bars in the BFO film matrix, the sub-skin and skin layers (see the legend in the inset of Fig. 4(a)) for both the type I and the type II surface. In the BFO film matrix, for both types of domains the c/a ratio is nearly constant at a value of ~1.02. Considering the lattice mismatch of ~0.3% between DSO 31 and BFO 22 , the measured tetragonality of 1.02 indicates that the BFO film is fully strained due to the epitaxial relationship with the DSO substrate. The c/a ratio remains almost unchanged from the matrix to the BFO skin layer for the type II surface area, while in the skin layer in the type I surface region this ratio increases abruptly to a value of ~1.13, which is much larger than the value for the BFO film matrix on the DSO substrate.
Accompanying the large c/a ratio in the BFO skin layer on the type I surface, the off-center displacements of the FeO columns with respect to the centers of mass of the surrounding four Bi atomic columns also show a dramatic increase, suggesting a strong coupling between the polarization and the lattice distortion. As shown in Fig. 3(d) and the inset to Fig. 4(a), in the BFO film matrix the central FeO columns (magenta) are shifted down and right along approximately the < 011> diagonal direction, with an average displacement |D BF | of ~30 pm ( Fig. 4(b)). This value is very close to a displacement of 33 pm for Fe atoms reported for BFO grown on a TbScO 3 substrate with a lattice mismatch of < 0.14% 32 . In contrast, the off-center displacements of the FeO columns in the BFO skin layer are ~50 pm, i.e., 67% larger than those in the film matrix. In addition, the displacement vector undergoes a reorientation towards the < 001> direction. We also measure shifts of the O atoms, as marked by open arrows in Fig. 3(d,e), which do not exhibit a collinear relationship with the shifts of the FeO columns.
In the domain with the type II surface, no changes in off-center displacements of FeO columns were measured, consistent with the measured variation in c/a ratio in Fig. 4(a).
In order to understand the atomic structures and the large difference in lattice expansion in the skin layer between the two types of surface, we performed first-principles DFT calculations (see Methods), starting with free-surface relaxation of BFO (supplementary Figure S4). Although the calculation exhibits a different pathway for free-surface relaxation compared to our experimental observation, extremely large atom movements still occur in the vicinity of a type I surface. This behavior reveals a strong but natural demand to lower the system energy by ionic displacements, leading to a compensation of accumulated charges on this type of surface.
Based on the experimental results shown in Figs 2 and 3, supercells were constructed, comprising 2 × 2 × 6 cells of BFO together with a structural unit of Aurivillius or δ -Bi 2 O 3 containing 8 O and 4 Bi atoms (i.e., a Bi 2 O 2 -like unit) at the surface as a reference model for a type I surface, in which the nominal valency of the full supercell is zero (see Methods and supplementary Figure S5). Considering the need to screen the positive charge induced by the matrix P s and the surface (BiO) + layer at the type I surface (Fig. 1), the surface structure was tuned to electronegativity by deliberately introducing Bi vacancies or O adatoms to the reference model. The resulting models were named BiX or OX (where X is an integer), corresponding to the residual Bi occupancy or the numbers of O adatoms at the surface (see Methods and supplementary Figure S6). Since it is very difficult to compare the total energies of two systems that contain different numbers of atoms, the c/a ratio of the BFO skin was used as a figure of merit for determining the best matching model for the type I surface. All of the relaxed structures are shown in supplementary Figure S6, alongside corresponding HAADF STEM image simulations. Our ab initio result confirms that, following an increase in negative surface charge in the Bi double-layer, the c/a ratio of the BFO skin increases accordingly in all of the relaxed models (supplementary Figure S7), in agreement with the experimental result and in accordance with the surface charge compensation mechanism. Furthermore, based on the good correspondence between the simulated and experimental HAADF STEM images (supplementary Figure S1), it is reasonable to conclude that the presence of outer surface O adatoms improves the stability of a type I surface. A closer look at the averaged ABF STEM image in Fig. 3(d) reveals weak contrast from O atoms at the outer surface, as indicated by the black arrow, which is reproduced well in the simulated image in Fig. 3(e) based on the relaxed O8 model only. Moreover, a comparison of the c/a ratio and the displacement vector D BF further confirms the match between the DFT calculations and the experimental measurements, as shown in Fig. 4. Based on these results, the relaxed O8 model is considered to be representative of the real structure.
Properties of the skin layer. The spontaneous polarization of the surface region of the BFO film was calculated on the basis of the atomic positions in the O8 model shown in Fig. 5(a). The variation in polarization with distance from the film matrix to the surface was evaluated by using the 'sliding' unit cell method 33 Figure 5(b) shows the modulus of P s plotted as a function of the 'sliding' BFO unit cell layer. In the film matrix (i.e., from layer 1 to layer 8), |P s | remains constant at ~0.867 C/m 2 . This value is close to a prediction of 0.873 C/m 2 34 . The polarization decreases slightly in the sub-skin unit cells, followed by an abrupt increase to ~1.4 C/m 2 in the skin unit cell layer, which is very close to that reported for the tetragonal phase 35 . The slight drop in the sub-skin layers arises mainly from a decrease of the in-plane components of P s . From the model, we note that the in-plane components of polarization in the skin unit cell are comparable with that in the matrix. The increase in the out-of-plane component results in a rotation of the polarization vector in the skin unit cell layer in the < 110> plane. This phenomenon is similar to the M A phase transformation of BFO 36 , suggesting a similar relationship between the bound charge and the lattice strain.
According to the DFT calculations, the reconstruction in the surface layer also affects the local magnetic properties of BFO. The results of our calculations show that BFO still stabilizes in a G-type antiferromagnetic structure in both the type I and the type II surface. This is expected, considering the fact that the magnetic G-type to C-type order transition for BFO films occurs at a c/a ratio of ~1.19 as a result of an increase in compressive lattice strain 37 . The calculated Fe magnetic moment in the ground state decreases from ~3.8 μ B in the strained film matrix to ~3.0 μ B in the skin unit cell layer for the type I surface. In contrast, the magnetic moment remains at 3.8 μ B for all of the Fe ions in the unreconstructed type II surface, which is comparable to the value of 4.11 ± 0.15 μ B in a BFO single crystal 38 .

Discussion and Conclusions
Our results reveal the detailed atomic structure in the surface region of a BFO film on a DSO substrate. Two types of surface structure exist, depending on the polarization direction of the underlying ferroelectric domains and the termination of the polar atomic planes at the surface. A structure similar to the Bi 2 O 2 unit of the Aurivillius phase is found on the surface of film domains that have an upward component of ferroelectric polarization (type I surface). The same structure may be present in an image shown in Ref. 23 for a BFO film grown on a different substrate with a different thickness by molecular beam epitaxy, although it was not mentioned in that paper. Our observations are also numerically consistent with a BFO surface structure viewed along the < 110> direction grown on a 5 nm La 0.7 Sr 0.3 MnO 3 /SrTiO 3 substrate 39 , although the surface contraction layer was then attributed to the formation of a BFO dead layer instead of Bi 2 O 2±δ . These findings indicate a general tendency of BFO film surfaces to reconstruct into this kind of structure.
The fact that the appearance of the two surface structures strongly depends on the polarization of the film domains implies a coupling between the structure and surface charges triggered by ferroelectric polarization. For the type I surface, the structural unit of the Aurivillius phase acts as a reservoir for negative charges from the excess of O in the outer Bi plane, which stabilizes the surface structure. The negative charges contribute to the compensation of positive charges associated with the upward component of polarization and the (BiO) + polar surface termination. In contrast, the negative surface bound charges induced by the downward component of P s on the type II surface can be self-compensated by the positive charges of the (BiO) + termination layer and therefore do not require an additional Bi 2 O 2±δ structure. Indeed, free surface relaxation of the R6 structure shows that only small changes in O column splitting from 52 to 76 pm in the {100} projected plane (marked by orange arrows in supplementary Figure S4) are needed to stabilize the structure of a type II surface. This change is presently below the detection limit of the applied experimental techniques.
Recently, Gao et al. 21 studied the surface of a ferroelectric PZT film. In their work, a charged 180° DW and suppressed polarization were found at negatively polarized surfaces (i.e., with P s pointing opposite to the film surface), while no reconstruction was detected at positively polarized surfaces (i.e., with P s pointing towards the film surface). These observations were interpreted in terms of surface screening, possibly due to the presence of O vacancies on the former surface and Pb vacancies on the latter surface, controlled by ferroelectric P s alone. In the case of a BFO film, due to differences in crystal structure, chemistry and stacking polar atomic planes, we observe a different reconstruction that depends in an opposite way on polarization. A prominent reconstruction then occurs at a positively polarized type I surface, while no reconstruction is detected for a negatively polarized type II surface. Our results provide clear evidence for the combined contributions of ferroelectric polarization and polar atomic plane termination. These contributions can be considered in terms of a capability to alter the surface atomic structure in a different way from that for a PZT film, which is controlled by P s only. Based on our results and discussions, it is reasonable to expect an analogous behavior for surface reconstructions in other ferroelectric oxides with odd valences of cations, such as LiNbO 3 .
A large c/a ratio (~1.25) and an enhanced value of |P s | (~1.5 C/m 2 ) in a tetragonal-like phase of BFO have already been reported 35 . This phase was found to be stabilized together with a rhombohedral-like phase by large compressive strain. For the rhombohedral-like phase, different c/a ratios were reported from 1.05 to 1.07 35,36 . In addition, Beekman et al. reported an S′ phase with a c/a value of 1.09 40 . Our study shows that the BFO skin unit cell layer has a comparable c/a ratio and magnetic ordering to the S′ phase, but possesses a very large value of P s , as in the tetragonal-like phase. Our results provide an important structural basis for tuning the local multiferroic properties of BFO by controlling electrical boundary conditions.

Methods
Sample preparation. An epitaxial BiFeO 3 thin film was grown on an epi-polished single crystalline DyScO 3 (110) o substrate (CrysTec GmbH, Berlin) by using pure-oxygen pressure radio frequency magnetron sputtering. This low-energy technique allows thin films with very high atomic structural perfection to be produced 41 . Prior to deposition, the chamber was pumped down to a background pressure below 1 × 10 -5 mbar. The substrate was then heated to the deposition temperature of 650 °C in an oxygen atmosphere at 2 mbar, which was maintained throughout film growth. Owing to the reduced loss of volatile Bi in the high oxygen condition, a stoichiometric BiFeO 3 target was used. The target was pre-sputtered for several minutes to remove possible contamination on the surface. Layer deposition lasted for 3 hours. After growth, the sample was cooled naturally to room temperature.
Cross-sectional lamellar specimens for STEM investigations were prepared by focused ion beam (FIB) milling using an FEI Helios Nanolab 400 s dual-beam system 42 . The lamellae were cut from two directions, < 100> and < 110> , to investigate the atomic details systematically. The as-prepared lamellae were carefully thinned using 2.5 kV Ar ion milling in a Bal-Tec Res-120 system, followed by final cleaning with a 500 eV Ar ion beam (Fischione Nanomill, Model 1040) to remove possible damaged layers introduced during the previous milling procedures.

Scanning transmission electron microscopy and image simulations. HAADF STEM imaging 25
was carried out in an FEI Titan G3 60-300 "PICO" microscope equipped with a high-brightness field emission gun (XFEG), a monochromator unit, a probe spherical aberration (C s ) corrector and a combined image C s -C c (chromatic aberration) corrector 43 . ABF STEM imaging 27 was performed on a JEOL JEM-ARM200F probe-C s -corrected microscope with a built-in ABF attachment. Atomic-resolution EDXS mapping was performed on an FEI Titan G2 80-200 ChemiSTEM microscope equipped with an XFEG, a probe C s corrector and a super-X EDXS system 44 . All of the microscopes were operated at 200 kV. The convergence semi-angle for STEM imaging was approximately 22 mrad, while the collection semi-angle was 12-24 mrad for ABF imaging and 70-176 or 200 mrad for HAADF imaging. All high-resolution ABF/HAADF images shown in this work are shown as raw data without any post-filtering.
Multislice ABF and HAADF STEM image simulations were performed using Dr. Probe software 45 on the basis of atomic models relaxed by first-principles DFT calculations and using experimental imaging parameters. Considering the unique site occupation in the Bi75, Bi25 and O4 models (caused by the size of the supercell), surface Bi disordering was added manually to obtain a correct surface intensity. However, this operation did not change the simulated details in the BFO skin layer. A series of sample thicknesses (up to 60 nm in steps of approximately 0.8 nm, i.e., two BFO pseudocubic unit cells) was tested. A thickness of 38-40 nm was finally utilized by matching the image contrast in the experimental and simulated images. The sample misalignment angle and the aberrations of the microscope were set to zero for the simulations. Structural models were visualized using VESTA software 46 . Ab initio simulations. First-principles DFT calculations within the local spin-density approximation (LSDA) 47,48 were performed using the Vienna ab initio Simulation Package (VASP) 49,50 . All results were obtained using the projector-augmented plane-wave method 51,52 by explicitly treating the valence electrons as follows: 15 for Bi (5d 10 6s 2 6p 3 ), 14 for Fe (3p 6 3d 6 4s 2 ) and 6 for O (2s 2 2p 4 ). The calculations do not include spin-orbit corrections. Brillouin zone integrations were performed with Gaussian broadening 53  Bulk calculations were performed with a 5 × 5 × 5 Monkhorst-Pack k-point mesh 54 centered at Γ and a 550 eV plane-wave cutoff, both of which result in good convergence of the computed ground-state properties. Structural optimizations were achieved by allowing the atoms in the unit cell to relax until all of the forces on each atomic site were below 5 meV/Å and simultaneously achieving a total energy convergence of 10 −6 eV. In order to correct for the metallic behavior observed in the LDA band structure, we applied the LDA + U approach described by Dudarev et al. 55 , in which only an effective Hubbard parameter U eff = U − J enters the Hamiltonian. The magnitude of U eff was varied between 0 and 7 eV for the Fe d states (the standard LSDA result corresponds to U eff = 0 eV). The effective Hubbard parameter U eff = 2 eV was used in all of our calculations. Higher values of the U eff parameter also yielded similar results for the structural properties.
The ground state of bulk BFO is R3c and exhibits G-type antiferromagnetism 22 . However, strain from the substrate is a key factor that can be used to tune its phase, as pointed out in recent studies 56 . In our case, the use of a DSO substrate lowers the symmetry of the BFO. We therefore used an orthogonalized coordinate system with P1 symmetry. We first transformed the experimental BFO coordinates (ICSD No. 109370, from Ref. 57) to the orthogonalized coordinate system by using CellMuncher.exe included in the Dr. Probe software package 45 . A 2 × 2 × 2 supercell containing 40 atoms of BFO was chosen, in order to accommodate all possible antiferrodistortive rotations of the FeO 6 octahedra and to enforce the experimental triclinic and monoclinic geometries for the different phases. We obtained a relaxed BFO bulk cell of 40 atoms with a = 0.7791 nm, b = 0.7785, c = 0.7784 nm and α = β = γ = 90°. A schematic diagram is shown in Fig. 1 in the main text. It should also be noted that the optimized lattice constants are underestimated with respect to experimental values as a consequence of using the LSDA, with a scaling factor of ~1.02 that is consistent with Ref. 58.
In order to simulate lattice strain imposed by the DSO substrate, we compressed the lattice constant along the x-and y-axes and fixed them to the calculated lattice constant of DSO (110) o . We then scanned along the z-axis and fully relaxed all of the atoms to obtain a local energy minimum. This approach created a pseudocubic phase with a = b = 0.773 nm and c = 0.7877 nm. The resulting c/a ratio is 1.02, which is in good agreement with our experimental observations. It also gives a ground state Fe magnetic moment of 3.8 μ B , which is slightly lower than the value of 4.0 μ B in the R3c phase. Such a pseudocubic unit cell was used to build the surface structure.
In order to construct a Bi-terminated surface, a 6-layer 2 × 2 supercell with BiO-terminated surfaces (denoted R6) was constructed as a reference system (supplementary Figure S4). An additional 1.5 nm layer of vacuum was used to exclude possible interactions between the two surfaces and a dipole interaction correction was applied perpendicular to the surface in all surface calculations. One layer of eight O atoms and one layer of four Bi atoms in an Aurivillius type configuration (supplementary Figures S5(c1) and Bi100/O0 model in supplementary Figure S6) were added on top of the R6 structure to simulate the atomic behavior observed experimentally. This model acts as a starting point, because the nominal valency of the system is 0. Switching the 8 O atoms on the sites between the Bi-Bi interlayer (Aurivillius-like configuration in supplementary Figure S5(c1)) and the outermost layer (δ -Bi 2 O 3 -like occupancy in supplementary Figure S5(c2,c3)) showed little influence on the c/a ratio, with a system energy for the latter case at least 0.75 eV higher than for the former case. Therefore, site selections were not considered in the following calculations.
Charge compensation was considered by removing Bi or adding O atoms on the basis of the Bi100/O0 model. As shown in detail in supplementary Figure S6, Bi75, Bi50 and Bi25 models were constructed by removing 1, 2 or 3 surface Bi atoms, while O1 to O8 models were constructed by adding 1 to 8 O atoms on the outer Bi surface on the basis of the Bi100/O0 model. Relaxation of the structures was carried out by fixing the atoms in the bottom 3 layers of the BFO and fully relaxing the other atoms. Convergence was obtained using a 3 × 3 × 1 Monkhorst-Pack special k-point grid with a cutoff energy of 550 eV, in order to terminate the plane-wave expansion. The criterion for relaxation convergence was the same as for the bulk calculations.