Domain walls in a perovskite oxide with two primary structural order parameters: first-principles study of BiFeO$_3$

We present a first-principles study of ferroelectric domain walls (FE-DWs) in multiferroic BiFeO$_3$ (BFO), a material in which the FE order parameter coexists with anti-ferrodistortive (AFD) modes involving rotations of the O$_6$ octahedra. We find that the energetics of the DWs are dominated by the capability of the domains to match their O$_6$ octahedra rotation patterns at the plane of the wall, so that the distortion of the oxygen groups is minimized. Our results thus indicate that, in essence, it is the discontinuity in the AFD order parameter, and not the change in the electric polarization, what decides which crystallographic planes are most likely to host BFO's FE-DWs. Such a result clearly suggests that the O$_6$ rotational patterns play a primary role in the FE phase of this compound, in contrast with the usual (implicit) assumption that they are subordinated to the FE order parameter. Our calculations show that, for the most favorable cases in BFO, the DW energy amounts to a several tens of mJ/m$^2$, which is higher than what was computed for other ferroelectric perovskites with no O$_6$ rotations. Interestingly, we find that the structure of BFO at the most stable DWs resembles the atomic arrangements that are characteristic of low-lying (meta)stable phases of the material. Further, we argue that our results for the DWs of bulk BFO are related with the nanoscale-twinned structures that Prosandeev et al. [Adv. Funct. Mats. (2012), doi: 10.1002/adfm.201201467] have recently predicted to occur in this compound, and suggest that BFO can be viewed as a polytypic material. Our work thus contributes to shape a coherent picture of the structural variants that BFO can present and the way in which they are related.


I. INTRODUCTION
Ferroelectric (FE) crystals are insulators displaying a macroscopic polarization that can be switched by the action of an electric field. The ground state of a bulk ferroelectric, where a unit cell is periodically repeated in space, is a perfectly ordered structure where the atomic patterns that give rise to the polarization are the same in every part of the crystal. However, real ferroelectric crystals need to minimize the electrostatic energy arising from their finite size (i.e., the electric field outside a ferroelectric sample must rapidly decay to zero), and as a result the formation of regions where the polarization points in different directions will occur. Like in the case of other ferroic materials such as ferromagnets or ferroelastics, these regions are called domains, and the area in between domains is called domain wall (DW).
The presence of DWs in a material can cause important changes to its macroscopic behavior. 1 Thus, for example, recent findings about the conducting 2,3 and photovoltaic 4 properties of ferroelectric DWs (FE-DWs) have fostered the interest in using them as functional parts in devices in nanoelectronics. 5 The vast majority of the recent discoveries on FE-DWs have featured multiferroic BiFeO 3 (bismuth ferrite or BFO), a material that has concentrated many experimental efforts during the past decade because of the promise it holds for room-temperature applications based on magnetoelectric effects. 6 For all these reasons, currently there is a strong interest in improving our atomistic understanding of FE-DWs, in particular in the case of BFO.
While the width of magnetic DWs is in the range of 10 to 100 nm, FE-DWs are believed to be much thinner. The atomic structure of ferroelectric walls can be studied experimentally using techniques such as high-resolution transmission electron microscopy, which permits a spatial definition of up to around 0.5Å. However, the atomic displacements that characterize FE-DWs are of the order of 0.02Å, and therefore direct imaging and interpretation of the structure of these walls is still a challenge. 7 Firstprinciples calculations have also been used to obtain information about FE-DWs. Some studies [8][9][10][11] were carried out on ferroelectric perovskites where the polarization is the only structural order parameter, such as BaTiO 3 and PbTiO 3 . They report that FE-DWs in these materials are atomically thin. Other recent first-principles studies 2,12 have focused on multiferroic BFO. The width computed for BFO's DWs is at most 1 nm, and some of them are reported to show a chiral pattern that is abstent in former studies.
In the rhombohedral phase of BFO that is stable at ambient conditions (R3c space group), the O 6 octahedra of oxygens that surround the Fe ions are rotated in antiphase about the polarization axis. These rotations constitute a second structural order parameter (OP) in addition to the polarization; such an O 6 -rotation related OP is usually referred to as anti-ferrodistortive (AFD). AFD arXiv:1211.5116v1 [cond-mat.mtrl-sci] 21 Nov 2012 modes constitute the most usual structural instability in the family of perovskite oxides, and have been studied in a systematic and extensive way starting with the classic works of Glazer 13 and others. Interestingly, AFD and FE distortions often compete in perovskite oxides, and it is rare to see them occur simultaneously in a material. 14 In fact, as far as we know, BFO is the only perovskite oxide that presents both OPs in its equilibrium phase at ambient conditions, which makes it a unique case for study.
Whenever two OPs occur simultaneously in a given phase of a material, we tend to assign the lables primary and secondary to them. Such a classification is useful to rationalize the structural phase transition and domain variants occurring in the crystal: The primary OP is the one that drives the phase transition and it would occur even in absence of the secondary OP; the primary OP also defines the number and properties of the structural domains that can form. In contrast, the occurrence of the secondary OP relies on the presence of the primary distortion; the secondary OP may appear as a result of the symmetry breaking caused by the primary OP (its character would thus be improper) or because of a destabilizing coupling with the primary OP (the secondary order would be triggered in that case). Interestingly, as we will see, such a distinction between primary and secondary OPs cannot be made in BFO's case.
Most of the literature on BFO focuses on the investigation of the FE polarization, which is effectively treated as the single structural OP of the compound. This is a rather natural approach: The spontaneous polarization, with its associated electrostatic energy, is the leading driving force for the formation of domains; additionally, we can easily act on it by applying an external field. As regards the AFD distortions, the oxygen octahedra are known to rotate specifically about the polarization axis, and are much more difficult to manipulate and characterize experimentally; hence, they are usually treated as if they played a secondary role. However, these AFD modes do not occur as a result of a polarization-induced symmetry breaking; further, first-principles calculations show that they constitute a genuine structural instability of the ideal cubic perovskite phase of BFO, independently of the occurrence of the polarization. Hence, it is clear that the AFD distortions do not fit the definition of a secondary OP given above. Moreover, the rhombohedral phase of BFO does not occur as a consequence of a displacive phase transition characterized by the condensation of the FE soft mode; rather, BFO's R3c phase undergoes a strongly discontinuous transformation upon heating, and the resulting high-temperature structure is probably characterized by oxygen-octahedra rotations more than it is by polar distortions. Hence, from this perspective either, there are no clear arguments to distinguish between primary and secondary OPs in BFO.
Detailed studies based on density functional theory (e.g., see Table I of Ref. 15) allow us to make this discussion more quantitative. It has been shown that the FE distortion of the ideal perovskite phase (which would lower the symmetry from the cubic P m3m space group to rhombohedral R3m) reduces BFO's energy by about 750 meV/f.u. (the precise result depends on the density functional employed in the calculations 15 ). In contrast, a pure AFD distortion (leading to the R3c space group) reduces the energy of the cubic phase by about 680 meV/f.u. The combined FE and AFD distortions render the R3c phase of BFO, which lies about 925 meV/f.u. below the cubic structure. We can thus conclude that the FE and AFD instabilities compete: If they did not interact at all, the energy of the phase combining both should be about 1430 meV/f.u. below the cubic structure (with 1430 = 750+680); the significantly smaller energy difference obtained from the simulations (925 meV/f.u.) clearly denotes a competition.
Hence, the simultaneous occurrence of FE and AFD distortions in BFO's R3c phase has to be attributed to the fact that, individually considered, both constitute very strong instabilities of the cubic phase of the material. Such instabilities clearly compete, not cooperate; yet, their simultaneous occurrence renders a net reduction of BFO's energy, thus overcoming the adverse effect of the competition. Let us note that, because the polarization and rotation axes coincide, one might be led to believe that the two OPs cooperate in some way, but such an interpretation would be wrong; rather, the coincidence of the two axes defines the conditions for which the FE-AFD competition is weakest.
In view of these facts, one should probably say that BFO presents two primary order parameters that need to be considered on equal footing in a theoretical discussion of this material. Bearing this in mind, we have conducted a first-principles study to reexamine the FE-DWs occurring in the R3c phase of bulk BFO. Interestingly, we find that the DW energetics is essentially determined by the type of discontinuity affecting the AFD patterns, and not so much by the change in the FE polarization accross the wall. Our observations have a rather natural explanation in the context of recent discoveries highlighting what we may call BFO's polymorphism, as we are able to relate the lowest-energy DWs found with the lowlying meta-stable phases that this compound is known to present. More precisely, the atomic arrangement at the most stable FE-DWs resembles the structure of a different BFO polymorph. Thus, our results suggest that multi-domain configurations in BFO display a peculiar form of polytypism. The connection of our work with the novel nanoscale-twinned structures 16 predicted to occur in BFO under various contidions (high temperature, high pressure, appropriate chemical substitution) are also discussed, offering a coherent picture of structural diversity and energetics in this material.
The paper is organized as follows. In Section II we discuss the details of the configurations we studied, introducing a rigourous classification of the FE/AFD DWs that can occur in BFO; we also describe the firstprinciples methodology used. In Section III we present and discuss our results, revealing the structural details that determine the DW energetics. Other aspects of the DWs, as e.g. their electronic structure, are also addressed, and the connection of our results with existing experimental information is discussed. In Section IV we summarize the work and present our conclusions.

A. Definition of the combined FE/AFD-DWs
BFO crystallizes in rhombohedral space group R3c, with a primitive unit cell that contains 10 atoms. Its rhombohedral lattice parameters are a rh = 5.6343Å and α rh = 59.348 • . 17 It is easier to visualize its structure using the pseudocubic 40-atom unit cell shown in Fig. 1. The Bi and Fe cations are displaced from their high symmetry positions along the [111] pc pseudocubic axis, and the material develops a polarization in this direction. At the same time, the O 6 octahedra surrounding Fe cations rotate about this [111] pc axis in antiphase, which corresponds to the so-called a − a − a − pattern in Glazer's notation. 13 Now, the cubic perovskite structure presents 8 equivalent directions of the 111 pc type. In a given domain, the electric polarization P can point along any of these 8 directions; further, once a direction for P is chosen, there are two equivalent ways (related by a phase shift) in which the O 6 octahedra can rotate. Hence, in total we can have 16 different FE+AFD domains in the rhombohedral phase of BFO. We now introduce a notation that will help us to describe the relationships between these domains. Let ω i be the three-dimensional pseudovector quantifying the rotations of the O 6 octahedron at cell i; the components ω iα , with α = x, y, and z, pertain to individual rotations about each of the pseudocubic principal axes of the perovskite structure, which are sketched in Fig. 1.
(From now on, all vectors and planes will be referred to the pseudocubic axes unless we mention otherwise. 18 ) As mentioned earlier, BFO's O 6 rotations around 111 are modulated in antiphase when we move from one cell to any of its first-nearest-neighboring cells. Such a modulation corresponds to the R q-point of the Brillouin zone of the ideal cubic 5-atom perovskite cell; for the particular case of rotations around the [111] axis, the pattern can be mathematically expressed as ω i = ω 0 (−1) nix+niy+niz (1, 1, 1), where the n iα variables are integers defining the location of cell i in the lattice.
Let us consider how a FE-DW can affect this ideal pattern. To do that, let us choose a cell i that is well inside the first domain; let us assume that ω i = ω 0 (1, 1, 1) and that the polarization in this first domain is also along [111], so that P I = P 0 (1, 1, 1). Then, let us pick a cell i that is well inside the second domain, and which can be reached from cell i by advancing an even number of cells along the pseudo-cubic direction α. Thus, we have n i β − n iβ = 2nδ αβ , n being an integer. Obviously, the two domains separated by the DW are related by symmetry, which implies that ω i and ω i must be symmetryrelated as well. We can have the following possibilities: (1) ω i = ω i , which coincides with what we would see in the monodomain situation, i.e., the pattern of O 6 rotations is essentially unaltered by the FE-DW. We will denote this case as a "AFD-DW of type 0". Note that, since in BFO the O 6 -rotation axis coincides with the direction of the spontaneous polarization, we can conclude that the polarization of the second domain must be either P II = P 0 (1, 1, 1) or P II = P 0 (−1, −1, −1); the former possibility implies there is no FE-DW ("FE-DW of type 0") and the latter corresponds to a 180 • FE-DW ("FE-DW of type 3", as the three polarization components change sign). We have thus identified the first two possible types of FE/AFD-DWs that can occur in BFO, and which we will denote, respectively, by "0/0" (i.e., FE-DW of type 0 and AFD-DW of type 0, which means that there is no discontinuity at all) and "3/0" (i.e., FE-DW of type 3 and AFD-DW of type 0, which means that the discontinuity only affects the polarization, which rotates by 180 • when moving from domain I to domain II).
(2) One rotation component changes sign with respect to the ideal situation. For the sake of the discussion, let us assume that the affected direction is z, so that ω i = ω 0 (1, 1, −1). This case, which we will denote as "AFD-DW of type 1", does imply a discontinuity of the a − a − a − rotational pattern. Indeed, we must have some sort of phase boundary for the O 6 rotations about the z axis as we cross the DW. As we will see later, our simulations indicate that such a phase boundary is the most trivial one that we could imagine: The z-rotations of neighboring cells accross the DW occur in-phase, thus breaking the anti-phase modulation sequence of the ideal pattern. Here again, we can deduce that the polarization of the second domain must be either P II = P 0 (1, 1, −1) or P II = P 0 (−1, −1, 1) implying that the FE-DW must be either of "type 1" (the 71 • FE-DW in which only one polarization component changes sign) or of "type 2" (the 109 • FE-DW with two components changing sign), respectively. We have thus identified two new types of FE/AFD-DWs, which we denote by "1/1" and "2/1", respectively.

B. Domain walls investigated
In this article we will consider DWs occurring at planes of low Miller indices, namely, (100) and (110), which are the ones that were discussed by previous authors. 2,12 Our initial task is to reduce the number of possible combinations of domains that we have to study. First, for a given DW many of these configurations will be equivalent by symmetry. Second, we limit ourselves to neutral DWs, i.e., those in which there is no discontinuity in the polarization component across the wall. Note that such a discontinuity would give raise to bound charges at the DW, which would lead to a large unfavorable electrostatic energy. 19,20 Hence, the electrostatically allowed cases are those in which (P I − P II ) · (1, 0, 0) = 0 and (P I − P II ) · (1, 1, 0) = 0, respectively, for the (100) and (110) planes.
Lastly, we could also limit ourselves to domains that fulfill the condition of mechanical compatibility. 21,22 To understand the origin of this condition, note that the ground state phase of ferroelectric perovskites involves a distortion of the cubic cell, i.e., an homogeneous strain.
In the presence of a DW, the distortion of unit cells at opposite sides of the DW may or may not coincide. If such strains coincide, we have mechanical compatibility; if not, an elastic energy penalty will occur. In the case of BFO, a possible mechanical mismatch at the DW comes Once the criteria described in the previous paragraphs are applied, it turns out that there are 6 possible inequivalent configurations for (100) DWs, all of them with P I = P 0 (1, 1, 1). Concerning (110) domains, there are two inequivalent ways to fix domain I: using P I = P 0 (1, 1, 1), or using P I = P 0 (1, −1, 1), which implies a domain where the polarization is parallel to the DW. For these two alternatives we find, respectively, four and eight inequivalent DW configurations. In total, we have 18 different DWs that we considered in this study. Table I lists them, indicating their FE/AFD-DW type and whether the mechanical matching condition is fulfilled or not.

C. First-principles methods
Our calculations are based on density-functional theory (DFT). 23 Most of them were done using the localdensity approximation (LDA) [24][25][26] for the exchange and correlation functional, although for comparison we repeated some calculations using two flavors of the generalized-gradient approximation: the functional of Perdew, Burke, and Ernzerhof 27 (PBE), and its adaptation to solids 28 (PBEsol). To obtain a better description of iron's 3d orbitals, we added a Hubbard-U term to the energy of the system, following the prescription of Dudarev and coworkers. 29 We used the implementation of this formalism in the Siesta code. 30 Siesta uses a basis set of localized numerical orbitals, which makes it computationally very efficient. In our case, we included in our basis set s (DZ), p (SZ), and d (DZ) orbitals centered at Fe; s (DZ), p (DZ), d (SZ), and f (SZ) orbitals centered at Bi; and s (DZ), p (DZ), and d (SZ) orbitals centered at O. The SZ and DZ terminology means that we used either one (singleζ) or two (double-ζ) radial functions for each of the angular forms, respectively. 30 We used norm-conserving Troullier-Martins pseudopotentials 31 for a efficient treatment of the interaction between ion cores and valence electrons; the pseudopotentials were generated in the following configurations: Fe's 4s 2 3d 6 , Bi's 6s 2 6p 3 , and O's 2s 2 2p 4 ; partial core corrections were included in all cases. We used a real-space grid to represent the electronic density, the Hartree and exchange-correlation potentials, and the matrix elements between basis orbitals; this grid is equivalent to one generated by a plane-wave basis cut off at 300 Ry. The integrations in reciprocal space were done using a Monkhorst-Pack grid equivalent to a 4×4×4 one for a 5-atom perovskite unit cell.
We also used the Vasp code 32 for our calculations.
Vasp uses plane-waves as basis sets, and therefore allows for a more systematic numerical treatment of the equations to be solved at the cost of more computational time. We used the projector-augmented wave method 33 to represent the interaction of electrons with ionic cores solving for the following electrons: Fe's 3s, 3p, 3d, and 4s; Bi's 5d, 6s, and 6p; and O's 2s and 2p. We generated a plane-wave basis set using a kinetic energy cutoff of 500 eV. The integrations in reciprocal space were done using the same k-point grids as in Siesta.
When using Vasp we set the Hubbard-U equal to 4 eV, a value determined by requesting the computed magnetic interactions to be in quantitative agreement with those obtained from calculations with hybrid functionals. 34 U ≈ 4 eV has become a frequent choice in first-principles studies of BFO with this same kind of formalism, as it leads to qualitatively and semi-quantitatively correct results for all the properties investigated so far. Due mainly to the differences in the treatment of the ion cores, with Siesta we need to use U = 2 eV to reproduce the Vasp results for the ground state of BiFeO 3 ; thus, we used this value for the Siesta calculations.
Unless mentioned otherwise, the results in our tables and figures are the ones obtained with Vasp. As a first step in characterizing our DWs we computed the energy needed to create them. In order to do so we designed an 80-atom periodic simulation cell for each configuration listed in Table I, containing two domains of 40 atoms and two DWs. We then wanted to find the minimum energy configuration for each case without making any prior assumption on the location of the DW in terms of parallel ion planes, nor on the specific atomic arrangement occurring at the DW. Thus, for each DW type considered, we started with some reasonable atomic configuration corresponding to the two domains involved, and performed a short molecular dynamics simulation, starting from random velocities, using Siesta. Our tests indicated that this procedure leads to an unrestricted and efficient exploration of the energy surface, allowing the system to find the lowest-energy DW configuration in a reliable and reproducible way. We then relaxed the lattice vectors and atomic positions until the stresses and forces were small. The configurations thus obtained were further relaxed using Vasp until the forces were below 0.01 eV/Å and the stresses were below 0.1 kbar. At all times the magnetic moments of the Fe ions kept an antiferromagnetic arrangement of type G (closest neighbors have antiparallel spins); this is the configuration that BFO shows experimentally at a local scale at ambient conditions. Table II includes all the DW energies that resulted after this process, computed as where E is the energy of the DW configuration, E 0 is the energy of bulk BFO (computed for the same supercell for better numerical accuracy), S is the surface area of the cell face parallel to the DW, and we divide by two because each cell contains two DWs. Before discussing our results, let us note that the E DW values that we obtain are significantly lower that those reported in the previous studies of Seidel et al. 2 and Lubk et al. 12 A comparison is shown in Table III, where we include our results for the lowest-energy DWs that involve polarization rotations by 71 • , 109 • , and 180 • , respectively. For consistency with the previous studies, all calculations reported in this Table were done with cells that contained 120 atoms, and we have restricted ourselves to cases satisfying the mechanical compatibility conditions. (Our DW energies obtained using 120-atom cells are very similar to the ones we calculated with 80atom cells, implying that the results of Table II are well converged in this regard.) There are some minor differences between the Vasp calculations done by the authors of Refs. 2 and 12 and ours; yet, we have tested that even if we modify our numerical approximations to make them as similar as possible to theirs, our results change very little. Moreover, Table III also includes the results we obtaind with Siesta, which are very similar to the values found with VASP. In addition, let us note that similarly low DW energies have been obtained by Wei et al. 35 in an investigation of FE-DWs in BFO thin films under various epitaxial strain conditions. What can explain the discrepancies with the previous published literature? As some of us have recently shown, 15 the energy surface of BFO presents a very large number of local minima, which seems to be mainly a consequence of the many bonding complexes allowed by Bi's chemical versatility. We reached exactly the same conclusions in the present investigation of BFO's DWs, as we found that many possible DW configurations can render a local minimum of the energy surface. Naturally, the physical relevance of high-lying local energy minima is questionable, and we found ourselves forced to implement a careful search for the global energy minimum for each DW type. Indeed, by thermally agitating and then quenching our configurations during the optimization process, we have had access to minima with (much) lower energy than the ones previously reported. Hence, our results stand a better chance of being physically relevant.
We found that the [111](100)[111] 2/1 DW configuration has the lowest DW energy. It occurs when two polarization components and one rotation component change sign at the DW (this is a 109 • DW from the perspective of the electric polarization). The relaxed atomic struc- ture that we obtained is shown in Fig. 2(a); here one can already appreciate that the DW discontinuity is very sharp, and the O 6 octahedra are not strongly distorted at or near the DW plane. Figures 3(a) and 3(b) contain the values of the displacements of Bi and Fe cations along the pseudocubic axes, relative to the center of mass of the corresponding surrounding O 12 dodecahedron and O 6 octahedron, respectively. Such distortions are the origin of the local polar dipoles giving raise to BFO's macroscopic polarization, and thus allow us to monitor the Pchange at the DW. The Fe/Bi displacements are plotted as a function of the distance between the corresponding cations and the center of the DW, which is assumed to be located in between the two planes of Fe ions whose environment is most different from the one they see in bulk BFO. We can thus see that these polar displacements are almost identical to the ones in bulk, except for small distortions in a very narrow area around the DW center. The O 6 rotations are quantified in Figs. 3(c) and 3(d); they are the three components of (−1) nix+niy+niz ω i , in the notation of Section II.A, with the two graphs containing information for the two inequivalent lines of neighboring O 6 octahedra that run towards the DW. The O 6 rotation patterns are also very similar to the bulk ones, except for the discontinuity of one AFD component at the DW.
This is zero for a regular octahedron, but takes a value of about 0.1 in bulk BFO, as the polar displacements cause variations in the lengths of the O 6 edges. We see that this measure shows very small differences in O 6 octahedra distortions with respect to the bulk ones. In all, the result of the detailed analysis of the atomic structure points to a narrow DW whose thickness is around 1 nm. Our finding that the most stable DW configuration is precisely [111](100)[111] 2/1 might seem somewhat arbitrary; on the contrary, this result leads to very suggestive conclusions. In this DW, the AFD component that changes sign is the one perpendicular to the wall; thus, if we focus on the two planes of octahedra next to the DW, we see that they display a a + b − b − Glazer rotation pattern, i.e., the O 6 tilts occur in phase (+ superscript) about the [100] direction perpendicular to the DW and in antiphase (− superscript) about the [010] and [001] directions within the DW plane. Further, the discontinuity in the FE polarization -which is mainly captured by the Bi displacements shown in Fig. 3(a) occurs in the plane of the DW, which results in an antipolar pattern of Bi displacements around the DW center. Interestingly, such structural features are exactly the ones characterizing the P nma phase of bulk BFO, which is experimentally known to occur at high temperatures 36 and under hydrostatic pressure, 37,38 and which is relatively close in energy to the R3c ground state according to previous first-principles calculations. 15 In fact, BFO's P nma phase would match the structure of our DW almost identically if we added to it a polar distortion along the direction of the in-phase oxygen-octahedra rotations.
Interestingly, first-principles theory suggests that such a ferroelectrically-distorted P nma phase, which would present the polar space group P na2 1 , may be a lowlying meta-stable polymorph of BFO, 15,39 although such a structure has not been observed experimentally yet. Hence, our results suggest that the lowest-energy DW found in BFO owes its stability to the fact that it can adopt an atomic arrangement that mimicks a low-energy polymorph of the bulk material.
Interestingly, our results for this DW are also reminiscent of the so-called nanoscale-twinned phases that have been predicted to occur in BFO under various conditions (e.g., high temperature, hydrostatic pressure, chemical substitution of Bi by rare-earth cations). 16 In fact, the novel phases described by Prosandeev Fig. 2(b). This is a 180 • FE-DW where the cation displacements go in opposite directions at opposite sides of the wall, as shown in Figs. 4(a) and 4(b). As regards the network of O 6 octahedra, there is no AFD discontinuity involved in this boundary; yet, as we can see in Figs. 4(c) and 4(d), because the lattice needs to accommodate the large change in the polar cation displacements, the O 6 rotation angles at the DW deviate noticeably from the bulk-like values. From Fig. 4(e) we notice that there are some Fe-O bond lenght variations accompanying these distortions. However, the O-O bond length dispersion is very similar to the bulk one, as illustrated in Fig. 4(f). Let us also note that Lubk et al. computed a very high energy of 829 mJ/m 2 for their most stable 180 • FE-DW which, as in our case, involves no AFD phase boundary; they also reported that the polarization rotates in two steps from one domain to the other, forming a sort of chiral pattern. Our finding of a sharp 180 • FE-DW with an energy of 82 mJ/m 2 and around 1 nm in thickness indicates that the configuration described by Lubk et al. is very unlikely to occur in

reality.
Among all the considered DW configurations, let us discuss in some detail two additional cases, namely, the [111](110)[111] 1/1 DW that is sketched in Fig. 2(c) and whose structural features are shown in Fig. 5, and the [111](110)[111] 1/1 DW characterized by the results in Fig. 6. Both of them are 71 • FE-DWs, their respective energies being 152 mJ/m 2 and 142 mJ/m 2 for the 80-atom unit cell. Interestingly, while the former satisfies the mechanical compatibility condition mentioned above, 21,22 the latter does not. Yet, in spite of this fact, the second DW has the lower energy. Indeed, generally speaking, our results show that the mechanical mismatch does not involve a large energy penalty. The relative small importance of this mechanical mismatch should probably be attributed to the small shear strains associated with the BFO's R3c phase, which result in a small deformation of the DW plane of little energetic consequence. Hence, one should probably not consider the mechanical-matching criterion in studies of rhombohedral BFO, and should be cautious about its assumption for other materials or phases displaying small shears.
From the three cases studied in detail so far, we see that low DW energies are associated to configurations in which the O-O bonds keep lengths similar to the ones displayed in bulk BFO. To further investigate this trend, We computed this measure for each of the DW configurations of Table II, and plotted the DW energy as a function of it in Fig. 7. We see a very clear correlation between low DW energies and low variations in the O-O bond lengths. We can further understand this trend by looking at the rotations that are present in the center of the domains, independently of the exact atomic rearrangements that take place at the DW. Let us consider for a moment a cubic perovskite crystal where we remove the A and B cations, and we keep the network of regular O 6 octahedra. We now impose the formation of an AFD-DW, so that at each side of it the octahedra rotate a fixed angle, in a − a − a − fashion, about one of the 111 axes. The oxygen atoms at the DW plane are shared by two octahedra, one from each of the domains separated by the wall. Then, when the two different a − a − a − -like patterns freeze in, the rotation of each octahedra will try to send  the oxygen atoms moving in one of four possible directions that form 90 • among them. Depending on which directions are chosen by the two octahedra sharing a single oxygen, three different cases can arise, as depicted in Fig. 7: the two directions may coincide (case A), they may form 90 • (case B), or they may form 180 • (case C). In order to accommodate the two conflicting movements required from a single oxygen ion, the octahedra at the DW will distort from their regular shape in the 90 • case, and more so in the 180 • case.
In a material where AFD and FE distortions are present, the O 6 octahedra are more distorted than in the case we have just described. However, if the rotation angle is relatively large, as in BFO, most of the displacement of the O atoms out of high-simmetry positions is due to the rotation about 111 . We can therefore apply the analysis of the previous paragraph to our BFO DWs, and label them following the same notation, since the kind of mismatch is the same for every pair of octahedra on opposite sides of the wall. We have done so in Table II. As expected, we see that type A DWs are the lowest in energy, followed by type B DWs; finally, we find that type C DWs tend to be the most unfavorable ones. This is also illustrated in Fig. 7.

B. Role of exchange-correlation functionals
In a previous article 15 some of us pointed out how the use of different approximations to the exchangecorrelation functional of DFT can affect the results of the total energy of BFO configurations. In particular, we found that the energy differences between phases (i.e., their relative stability) depends significantly on the employed energy functional. Therefore, since we may in principle expect a similarly important dependence of the calculated DW energies, we repeated our calculations for three representative DW configurations using two other functionals: PBE 27 and PBEsol. 28 Table IV shows small variations in the numbers obtained for the DW energies, while the relative ordering of the configurations is always the same. The variations are smaller than those in Ref. 15 because our DW configurations are all relatively similar, while in the cited work phases with very different structures (including, e.g., supertetragonal phases with O 5 pyramids instead of O 6 octahedra) were compared. We take these results as further confirmation that DFT predicts that the most stable DW configurations are those that keep the O 6 octahedra bond distortions to a minimum, and that their DW energies are in the range of several tens of mJ/m 2 .

C. Electronic structure
Seidel and collaborators 2 were the first to report the experimental observation of an enhanced electronic con-ductivity at room temperature in BFO DWs. In their article, they also included results about the electronic properties of the DWs obtained using first-principles calculations; they mentioned that at the center of the domains the local density of states (LDoS) resembles the one of the bulk, while as the DW is approached the bandgap is reduced by an amount of up to 0.2 eV depending on the DW orientation. They correlated this reduction with the presence of the conductivity at the walls. Figure 8 shows the LDoS that we obtained from the atoms that are closest to the DW in each of the configurations of Tables III and IV. As we can see, if there are bandgap closings in these structures they are too small to be captured by this first-principles methodology. In these conditions, our interpretation is that conductivity at the DWs is not caused by any significative band-gap closing resulting from structural changes at the atomic level near a DW of pure BFO.

D. Connection with experiment
Seidel et al. 2 reported an enhanced conductivity in FE-DWs of types 2 (109 • ) and 3 (180 • ), while they found that FE-DWs of type 1 (71 • ) did not conduct in their experiments. On the other hand, the experiments of Farokhipoor and Noheda 3 revealed an enhanced conductivity at type 1 FE-DWs as well. The first set of authors argued that the observed conductivity was consistent with changes in the BFO structure at the domain wall; their complementary first-principles calculations showed a small reduction in the bandgap of pure BFO that they correlated to some degree with the conducting behavior. On the other hand, Farokhipoor and Noheda concluded that conduction at type 1 FE-DWs was related to the presence of oxygen vacancies. In another paper, Seidel and coworkers 40 also demonstrated that the wall conductivity in La-doped BFO could be controlled through chemical doping with oxygen vacancies.
Our results support a mechanism for the conductivity that goes beyond the (small) effect that DW-specific structural distortions have in the local electronic properties. As mentioned above, we find negligible differences between the electronic structure of the DWs and that of the bulk material. Further, we also investigated the exchange couplings between iron spins at the DW, in order to check whether a possible reduction in the magnitude of the magnetic interactions could have an effect on the spin structure (which might, e.g., be associated with a reduced Nèel temperature at the DW) and thus on the electronic structure and conductivity; however, the observed variations never exceeded a 10 % of the bulk exchange couplings, suggesting that the effect in the magnetic (and electronic) structure will be negligible. Hence, our calculations strongly indicate that the observed enhanced conductivity is to be linked to extrinsic factors such as oxigen vacancies or other sources of off-stoichiometry at the DWs. In this sense, we are currently performing calculations to compute the energy of formation of various defects at several types of DWs, wanting to determine (i) their effect in the electronic structure and (ii) their preferencial occurrence in specific DWs. We believe that this is a promising route to explain the experimental observations of a relatively large DW conductivity from firstprinciples theory.
According to the results in Table II, the creation of type 2 and type 3 DWs in bulk, defect-free BFO is favored energetically over the creation of type 1 DWs. On electrical switching at high field experiments, Seidel et al. 2 were able to create the three types of domains on 100-nm-thick epitaxial BFO films grown on a (110) SrTiO 3 substrate. On the other hand, Farokhipoor and Noheda 3 mention that their pulsed-laser deposition films grown on SrTiO 3 -(001) substrates showed mostly DWs of type 1. Such a preferential occurrence of 71 • FE-DWs, which had also been observed in previous works, 41 is in obvious conflict with our preditions and points at the importance of factors that were not included in our simulations. Indeed, the experimental results are for thin films grown in particular orientations and conditions of misfit strain, and in the presence of a bottom substrate (and electrode) that plays a role as important as to induce a symmetry breaking in the possible orientations of the polarization (i.e., from the 8 111 polarization variants that can in principle occur, the films of Ref. 3 only exhibit 4 types of domains with the polarization always pointing towards the substrate). As regards the misfit strain, the recent first-principles investigation of Wei et al. 35 suggests that the DW hierarchy we obtained for the bulk case should be preserved for films grown on compressive SrTiO 3 -(001) substrates (with a misfit strain of about −1.5%). Hence, the discrepancy between our predictions and the experimental observations seems to indicate that there are other factors beyond the epitaxial strain -e.g., the presence of defects, the effect of interfaces with substrate and/or electrodes, the specific electric boundary conditions -that play a key role in determining which DWs occur spontaneously in BFO films.
Our first-principles results point to a picture of extremely thin ferroelectric DWs around 1 nm thick, in agreement with previous calculations for other perovskite oxides. Our ideal DWs lie on infinite parallel planes, and they do not interact with each other when they are more than a few nanometers apart. (As evidenced by the rapid convergence of the computed DW energies as we increase the separation between DWs in our simulation box.) Instead, real ferroelectrics show DWs in equilibrium that exhibit a characteristic roughness, 42 indicating that there are defects in the lattice that are pinning the walls. 5 In this respect, our simulations can be taken as the initial step in building a model of ferroelectric domains that would take into account the role of size effects and defects, a challenge that is currently out of reach of direct first-principles calculations.

IV. SUMMARY AND CONCLUSIONS
We have presented a theoretical study of domain walls (DWs) in perovskite oxides where two primary order parameters coexist. We have taken BiFeO 3 as a sample case in which two structural instabilities of similar strengthi.e., a ferroelectric (FE) one, and one that involves concerted rotations of the O 6 octahedra (anti-ferrodistortive or AFD) -are present in the phase that is stable at ambient conditions. We have used density-functional theory to perform atomic scale simulations of the hybrid domain boundaries involving both FE and AFD discontinuities, and have thus revealed several novel aspects that need to be taken into account in situations like this one. First, and in contrast with what is implicitly assumed in most BFO studies, we have found that the DW energetics is essentially determined by the type of discontinuity in the AFD patterns, and not so much by the change in the FE polarization accross the wall. Thus, the AFD distortions, which are usually treated as the secondary (i.e., less important) order parameter in BFO, turn out to play the leading role in determining which DWs are most stable. Second, we have found that the lowest-energy DWs present atomic structures that can be directly identified with low-lying meta-stable BFO phases. Indeed, our results suggest that multi-domain configurations in BFO are similar to polytypic structures with (meta-)stable polytypes occurring at the (very thin) domain boundaries. This finding has an appealing and natural connection with the novel nano-twinned structures 16 that have been recently proposed to explain several controversial regions in the phase diagram of BFO and BFO-based solid solutions. Finally, we have discussed the implications of our results as regards current experimental work on BFO's DWs, emphasizing the need of considering various extrinsic factors to explain some of the most interesting experimental observations.