First-principles Landau-like potential for BiFeO$_3$ and related materials

In this work we introduce the simplest, lowest-order Landau-like potential for BiFeO$_3$ and La-doped BiFeO$_3$ as an expansion around the paraelectric cubic phase in powers of polarization, FeO$_6$ octahedral rotations and strains. We present an analytical approach for computing the model parameters from density functional theory. We illustrate our approach by computing the potentials for BiFeO$_3$ and La$_{0.25}$Bi$_{0.75}$FeO$_3$ and show that, overall, we are able to capture the first-principles results accurately. The computed models allow us to identify and explain the main interactions controlling the relative stability of the competing low-energy phases of these compounds.


I. INTRODUCTION
Magnetoelectric multiferroics, materials that simultaneously show magnetic and electric orders, are of significant interest, since the coexistence and coupling of these orders hold great potential for development of multifunctional devices 1,2 . BiFeO 3 is among the most exciting and extensively studied representatives of this family because it displays both orders at room temperature 3 .
Ferroelectricity appears in BiFeO 3 at T C ∼ 1100 K 4,5 . Below T C , it has a rhombohedrally distorted perovskite structure (space group R3c, #161) 6,7 , which differs from the perfect cubic phase by the presence of two distortions: (i) polar displacements of Bi 3+ and Fe 3+ cations with respect to O 2− anions (Bi 3+ dominates due to its stereochemically active 6s lone pairs 8 ) giving rise to a spontaneous polarization P of up to 100 µC/cm 2 along a pseudocubic 111 direction 9,10 ; and (ii) antiphase rotations R of the FeO 6 octahedra about the same pseudocubic 111 direction as the polarization (a − a − a − in Glazer's notation 11 ) 12,13 . (In the following, all directions are in the pseudocubic setting.) Below T N ∼ 640 K, BiFeO 3 also shows G-type antiferromagnetic (G-AFM) order with the nearest-neighboring Fe spins antialigned 4,14 . The canting of the Fe spins driven by Dzyaloshinskii-Moriya (DM) interaction 15,16 can give rise to a weak magnetization in this material. The DM interaction relies on the symmetry breaking caused by the FeO 6 octahedral tilts of BiFeO 3 12 ; indeed, the phase of the octahedral rotations defines the sign of the DM vector and, in turn, that of the weak magnetization. In bulk BiFeO 3 an incommensurate cycloidal spiral is superimposed on the G-AFM order, yielding a zero net magnetization 17 . This cycloid, however, can be suppressed by doping in bulk systems 18 and by epitaxial constraints in BiFeO 3 films [19][20][21][22] . Therefore, ferroelectricity can coexist with weak ferromagnetism in BiFeO 3 at ambient conditions. Additionally, a 180 • deterministic switching of the DM vector and weak magnetization by an electric field has been reported from a combined experimental and theoretical study of BiFeO 3 films grown on DyScO 3 substrates 22 . It is proposed that the magnetoelectric switching is the result of a peculiar polarization reversal that is found to occur in two steps, a 109 • rotation followed by a 71 • rotation (or vice versa); further, the FeO 6 octahedral tilts are believed to reverse together with the polarization, resulting in the observed reversal of the weak magnetic moment. Note that octahedral tilts will typically not follow polarization in a single-step 180 • reversal and, therefore, a two-step switching path is crucial for controlling the weak magnetization in BiFeO 3 by an electric field. These observations make BiFeO 3 a promising candidate for applications in magnetoelectric memory elements. However, to be technologically relevant, switching characteristics have to be optimized such that coercive voltages are below 100 mV and switching times fall in the range of 10-1000 ps 23,24 . Hence, the current challenge is to optimize the ferroelectric switching in BiFeO 3 while retaining the two-step path and magnetoelectric control.
One of the efficient strategies for optimizing polarization switching in BiFeO 3 is doping by La. Indeed, since polarization in this compound largely originates from the 6s lone pairs of the Bi 3+ cations, their substitution by isovalent, lone-pair-free cations leads to a reduction of the polar distortion 3,25,26 . For example, it has been experimentally demonstrated that 15-20% La-doped BiFeO 3 films show a polarization which is up to 60% smaller than that of pure BiFeO 3 films 24,27 . Further, first-principles calculations have predicted that subsitution of Bi by La cations reduces the energy barrier between polar states by up to 50% for 25% doping. This, in turn, leads to a reduction of coercive voltages (down to 0.8 V for a 100 nm film), enabling low-power switching 24 . Additionally, a significant reduction of switching times has been demonstrated for La 0.15 Bi 0.75 FeO 3 films compared to pure BiFeO 3 in a wide range of applied electric fields 28 . Nevertheless, further improvement requires understanding the origin of the two-step polarization switching in BiFeO 3 and related materials, as well as search for other strategies for manipulating the switching energy landscape. For that purpose, dynamical simulations of polar- ization switching based on phenomenological models of the free energy can be very helpful.
Landau free-energy potentials [29][30][31][32] , together with the Landau-Khalatnikov time-evolution equation 33 , offer a practical scheme to investigate switching in ferroelectrics. In this approach, one expands the energy of the compound around the reference paraelectric phase in powers of the relevant order parameters, keeping only terms compatible with the crystal symmetry 34 . It is important to note that the reliability of such simulations depends on the choice of the free energy expansion's coefficients, which can be obtained either by fitting to experimental data or from first-principles calculations 34 . In compounds as complex as BiFeO 3 , which feature multiple primary order parameters, deriving a suitable Landau potential from experimental information is all but impossible; hence, there is a clear need for the development of first-principles approaches.
In this work we introduce the simplest, lowest-order Landau-like potential able to reproduce the energies and structures of the low-energy polymorphs of BiFeO 3 and related materials. We present an analytical approach to compute the model parameters from density functional theory (DFT) and apply it to BiFeO 3 and La 0.25 Bi 0.75 FeO 3 . We demonstrate the overall accuracy of the obtained potentials, and discuss an effective way to treat intermediate compositions. Finally, we discuss the physics captured by the model, namely, the interaction between polarization and FeO 6 octahedral tilts, how it affects the energetics of different BiFeO 3 polymorphs, as well as the effects of La doping.

II. COMPUTATIONAL DETAILS
All calculations are performed using the DFT 35,36 implementation in the Vienna Ab initio Simulation package (VASP) 37 . For the exchange-correlation potential, we use the generalized gradient approximation optimized for solids 38 , with a Hubbard U correction (within Dudarev's scheme 39 and U = 4 eV) for a better treatment of iron's 3d electrons. We treat the interaction between core and valence electrons by the projector-augmented plane wave method 37 (5d 10 6s 2 6p 3 ), 9 of La (5p 6 6s 2 5d 1 ), 14 of Fe (3p 6 3d 7 4s 1 ), and 6 of O (2s 2 2p 4 ). We use a plane-wave basis set with a cutoff energy of 500 eV. We use a 3 × 3 × 3 Γ-centered Monkhorst-Pack k-point grid for reciprocal space integrals in the Brillouin zone corresponding to a 40-atom cell that is a 2×2×2 multiple of the 5-atom perovskite unit (see Fig. 1(a)). We ensure that these choices provide a good level of convergence for the quantities of interest. All simulations are performed with the G-type antiferromagnetic order of Fe magnetic moments imposed. In the lattice optimizations, the structures are considered to be relaxed when the forces acting on the atoms are below 0.01 eV/Å. We calculate elastic constants by finite differences using the strain-stress relationship 41 .

A. Landau-like potential
In this section we introduce the potential for BiFeO 3 and La-doped BiFeO 3 as an expansion around the reference paraelectric cubic phase in powers of the following order parameters: (i) the three-dimensional electric polarization P = (P x , P y , P z ); (ii) the antiphase rotations of FeO 6 octahedra R = (R x , R y , R z ); (iii) the strain η = (η xx , η yy , η zz , η yz , η xz , η xy ), where η xx = ǫ xx , η yy = ǫ yy , η zz = ǫ zz , η yz = 2ǫ yz , η xz = 2ǫ xz and η xy = 2ǫ xy , and ǫ ij are the components of the homogeneous strain tensor. The resulting expression for the potential (per perovskite unit cell) is written as follows: Here, F 0 is the free energy of the reference cubic phase. F (P ), F (R) and F (η) are the energy contributions solely due to polarization, FeO 6 octahedral rotations and strain, respectively, which we write as follows: We truncate the expansion in P and R at the fourth order, which is the minimum required to model structural instabilities. In turn, we only consider harmonic terms for the strains. Then, F (P, R), F (P, η) and F (R, η) are the coupling terms, which we write as: F (P, η) = γ P 111 (η xx P 2 x + η yy P 2 y + η zz P 2 z )+ γ P 122 (η xx (P 2 y + P 2 z ) + η yy (P 2 z + P 2 x ) + η zz (P 2 x + P 2 y ))+ γ P 423 (η yz P y P z + η xz P z P x + η xy P x P y ); Here we restrict ourselves to the lowest-order symmetryallowed couplings between the considered order parameters. Note that, in these equations, A, B, C, C ′ and γ are the material-dependent expansion coefficients that we compute using DFT as detailed in Sec. III B. The coefficients C 11 , C 12 , and C 44 in Eq. (4), as well as γ parameters in Eqs. (5)-(7) are given in Voigt notation for compactness.

B. Computing the potential parameters
We now describe the approach to compute the expansion coefficients of the Landau-like potential introduced in Sec. III A. We mainly focus on an analytical approach, but also discuss briefly a numerical scheme for compari-son.

Training set
We first identify the states or polymorphs that we want our models to describe. We consider the ground state as well as the low-energy polymorphs of the material, including the states that might be relevant for polarization switching. We thus define a training set of first-principles results corresponding to the energies and structures of such polymorphs.
Before we continue, let us introduce a convenient notation for the polymorphs we consider: we write "P ", where the first letter, P or R, indicates whether the structure presents a polar distortion or FeO 6 octahedral tilts, respectively (if both distortions appear, we indicate it by P+R); then, [001] or [111] shows the axis along/about which the corresponding distortion is oriented; finally, "c" indicates that the cubic cell is kept fixed. Thus, for example, the polymorph P[001]c is characterized by a polar distortion along the [001] direction and its cell is fixed to that of cubic reference structure. For simplicity, we also introduce short notations for all the polymorphs of interest, such as "1c" for the state P[001]c. We summarize all the notations for the polymorphs and the corresponding order parameters in Table  I.
BiFeO 3 The starting point for constructing the training set is the already-mentioned 40-atom supercell compatible with the G-type antiferromagnetic order and the antiphase rotations of the FeO 6 octahedra. First, we run a DFT simulation to optimize the volume of the cubic phase of BiFeO 3 using this supercell. Next, we use the optimized structure to construct six polymorphs (1c to 6c in Table I) by imposing the polar distortion and/or antiphase octahedral rotation along/about either the [001] or [111] directions while keeping the volume and shape of the supercell fixed (the corresponding ionic displacement patterns are illustrated in Fig. 2). We use DFT to optimize the ionic positions in these polymorphs and calculate the energies E s of the resulting structures, where the index s labels polymorphs in the training set. Additionally, we also consider the state we call P[111]+R[111]c (7c); here, we impose a polar distortion and octahedral tilts with amplitudes typical of BiFeO 3 , but oblique to each other. This structure does not correspond to a special point of the energy landscape; therefore, we do not perform a structural optimization and only compute its energy, which is needed to obtain the coefficient C ′ P R , as we will show in Sec. III B 2.
Next, we consider the first six polymorphs mentioned above (structures 1 to 6 in Table I), but now allowing for changes in the shape and volume of the supercell (note we omit the "c" in the notation).
In all cases we extract the displacements u Bi,s (u Bi,s are in Angstrom) of the Bi cations with respect to the corresponding O anion cages. We average the values of these Table I. Polymorphs (labeled by s) included in the training set for computing the potential's coefficients and the notations for their polarizations Ps, FeO6 octahedral rotations Rs and components of the strain tensor ηs.
displacements over all Bi ions to obtainū Bi,s . Since the Bi off-centering largely determines the electric polarization in BiFeO 3 , we estimate P s for the considered polymorphs as polymorph is P 6 = P 0 (1, 1, 1) which gives the magnitude of P 6 around its experimentally determined value of 1 C/m 2 . Similarly, we compute the rotation angles of the FeO 6 octahedra R s about the pseudocubic axes, from which we obtain the amplitude of the antiphase tilt pattern,R s . Finally, in the cases where the shape and volume of the cell are allowed to relax, we extract also the components of the strain tensor η s . The obtained results constitute our training set, which is presented in Table S1 of the Supplementary Material.
La-doped BiFeO 3 Experimentally, La dopants distribute quasi-randomly in the BiFeO 3 lattice, so that the macroscopic symmetry (cubic for the paraelectric phase, rhombohedral for the ground state) is only recovered when a sufficiently large sample volume is considered. Unfortunatately, reproducing such a situation in a DFT calculation has a prohibitive computational cost; thus, here we assume that a particular highly-ordered La arrangement, where the dopants are as separated as possible from one another and which respects the cubic symmetry of the reference lattice, is a good approximation to the average experimental configuration. (For a 25 % La doping, the symmetric arrangement we use is shown in Fig. 1(b).) This approach allows us to derive Landau potentials for doped materials, with the experimentally relevant symmetry properties, from relatively inexpesive DFT calculations. Admittedly, a careful (computation-ally costly) validation of its accuracy remains for future work.
Having chosen a suitable, symmetric dopant arrangement, we optimize the cubic cell of the reference paraelectric structure using DFT. Next, we use this structure to construct the sets of polymorphs 1c to 7c and 1 to 6, in analogy to the case of pure BiFeO 3 . For the case of a 25 % La composition, the obtained values of P s , R s , η s and E s of their optimized structures are summarized in Table S2 of the Supplementary Material.
Note that we encountered difficulties in constructing the training set for La 1−x Bi x FeO 3 compositions with an intermediate content of La (0 < x < 0.25). For example, for La 0.125 Bi 0.875 FeO 3 , one can easily construct a paraelectric reference by subsituting a single Bi atom in the supercell of Fig. 1(a). However, we observed that the P[001]+R[001]c and P[001]+R[001] polymorphs relax to the lower symmetry phases displaying additional (and large) in-phase rotations of the FeO 6 octahedra. These extra distortions are secondary modes activated by the symmetry breaking associated to the combination of polar and antiphase orders together with the considered arrangement of La dopants. These distortions are not expected to occur experimentally, as the La dopants are largely disordered in real samples, and such in-phase tilts may occur locally at most. Moreover, they cannot be treated within our simple potentials (an explicit consideration of in-phase tilts would be required) and complicate the definition of the training set. Hence, here we do not compute models for such intermediate compositions.
Nevertheless, as we show in Sec. IV C, suitable potentials can be obtained by interpolation between those obtained for neighboring (well-behaved) compositions. Table II. Conditions used to derive the analytical expressions for the potential's coefficients. Es|φs,eq denotes the energy of the polymorph s corresponding to the equilibrium value of order parameter φs. λ = 1 m 4 deg 2 /C 2 is an ad hoc coefficient used to balance the units for the terms in f5c and f6c (see text).

Conditions
The approach introduced in this section allows full control of the information used to compute the parameters of the potential. To achieve that, we derive an analytical expression for each of the potential coefficients, in terms of E s , P s , R s and η s of the polymorphs in the training set (see Sec. III B 1). To obtain such formulas, we use Eqs. (1) -(7) and impose the zero-derivative condition ∂F/∂φ i,s = 0, where φ i,s is the ith component of order parameter φ evaluated for polymorph s. Let us illustrate our procedure by presenting in detail the case of the parameters A P , B P and C P of Eq. 2.
By taking the derivatives of these energies with respect to P 1c and P 2c , and setting them equal to zero, we obtain the following equations for the equilibrium values of P 1c and P 2c : and Then, by using these expressions in Eqs. (8) and (9), we obtain, respectively, E 1c and E 2c as functions of the parameters of the potential, such as and From Eq. (11) one can see that 3B P + C P = −A P /2P 2 2c . By using this in Eq. (13), one can straightforwardly obtain the analytical expression for A P : From Eq. (12), in turn, one can obtain: Finally, by combining Eqs. (13), (14) and (15), we get: Since E 1c , E 2c , P 1c and P 2c are known from our first-principles calculations described above, the coefficients A P , B P and C P can be directly computed using Eqs. (14), (15) and (16), respectively. Similarly, we can derive the analytical expressions for the remaining coefficients of our potential. The specific conditions and properties used in the derivation are summarized in Table II, and the resulting expressions are: where B 7c = (2P 2 7c,⊥ + P 2 7c, )(2R 2 7c,⊥ + R 2 7c, ), and where λ = 1 m 4 deg 2 /C 2 is an ad hoc coefficient that allows us to combine two zero-derivative conditions (for polarization and tilts, respectively) into only one. Further, we have 7c,⊥ + 2R 2 7c,⊥ R 2 7c, )); (26) and Note, that it is possible to choose other conditions, different from those in Table II includes the information about the ground state and corrects this problem. These difficulties reflect the simplicity of our low-order polynomial model, which can account (exactly) for only a small number of properties. Finally, the elastic constants C 11 , C 12 and C 44 are calculated directly from DFT.

Numerical approach
The numerical approach that we introduce in this section allows to compute the potential coefficients using the information from all considered structural polymorphs. We focus here on the case of pure BiFeO 3 , noting that exactly the same procedure can be applied to La 0.25 Bi 0.75 FeO 3 .
We work with the BiFeO 3 polymorphs from the training set introduced in Sec. III B 1, namely, 1c to 7c and 1 to 6 of Table I. Based on the energies (E s ) and structural parameters (P s , R s , and η s ) obtained from DFT, we construct an overdetermined system of linear equations, with the potential parameters as unknowns, using the expressions for the energy and zero derivatives corresponding to all polymorphs. (For the 7c state, the zero-derivative condition does not apply. Also, we use the elastic constants C 11 , C 12 and C 44 directly obtained from DFT. ) We find that the potential obtained using this approach provides less accurate predictions for the properties of the low-energy polymorphs of BiFeO 3 as compared to the analytical approach introduced in Sec. III B 2 (see details in Sec. SII of the Supplementary material). More specifically, the numerically-determined potential does a better job at reproducing some features (e.g., the polarization of the P[001] state) that we disregard in our analytical approach. In turn, it is less accurate when it comes to capture some critical properties (e.g., the ground state energy). We conclude that, while this fitting approach might work well for more complete, higher-order potentials (for example, such as the one introduced in Ref. 42), it seems less suitable for computing the parameters of our low-order model. Therefore, in the following, we are going to discuss only the results obtained using the analytical approach.

A. BiFeO3
In this section, we analyze how accurately the potential introduced in Sec. III A predicts the properties of BiFeO 3 polymorphs. We begin by computing the coefficients of the potential following the analytical approach described in Sec. III B 2. The resulting values are presented in Table III. Next, we use the computed potential to calculate the equilibrium properties (P s , R s , η s , and E s ) of the polymorphs 1c to 6c and 1 to 6. Since in these polymorphs the form of P s is either (0, 0, P s ) or (P s , P s , P s ) (the same holds for R s ), in the following we will discuss single components of P s and R s (P s and R s , respectively). We plot the values of P s , R s and E s predicted using our potential versus their DFT counterparts as shown in Fig. 3 (all these values, as well as the components of η s are also presented in Table S1 of the Supplementary Material). We note that, if the model prediction and DFT value match exactly, the corresponding point lays on the black dashed line.
First, we discuss the BiFeO 3 polymorphs with the fixed cubic cell (no strain relaxation). From Figs. 3(a) and 3(b) one can see that, for these polymorphs, our model predicts P s and R s in nearly perfect agreement with DFT. As shown in Fig. 3(c), it also reproduces accurately their energies and, therefore, their relative stability. Next, we consider the BiFeO 3 polymorphs with allowed strain relaxation (Figs. 3(d)-(f)). In this case, our potential also provides accurate predictions for all considered quantities for the most of the considered polymorphs; in particular, it yields the correct ground state of BiFeO 3 (P[111]+R [111]). There is only one polyrmorph for which the model is less accurate, namely, P[001]. In this case, the DFT optimized structure has a large distortion along the z axis (the c/a ratio is approximately 1.27), accompanied by a large P z ; this is usually called supertetragonal phase 43,44 . This behavior is not well captured by our potential, as it underestimates the polarization and strains components compared to the DFT values (P z = 1.039 versus 1.624 C/m 2 ; η xx = −0.012 versus −0.044; η zz = 0.065 versus 0.216). This issue is also reflected in the energy predicted for this phase. From the DFT results one can see that, among the polymorphs with only polar distortion, the strain relaxation stabilizes the supertetragonal phase over the rhombohedral one (P[001] is lower than P[111] by 0.023 eV/f.u). Our model does predict the energy lowering of P[001] state due to the strain relaxation (the negative γ P 111 coupling results in large P z and η zz ). However, since it underestimates P z and η zz , this energy reduction is not enough to stabilize P[001] over P [111]. Note that these deficiencies were to be expected, as we decided to use a minimal amount of DFT information on the supertetragonal phase when deriving the parameters of our model (see Table II), because this state is not relevant for our ultimate purpose of studying polarization switching in the rhombohedral phase of BiFeO 3 . Moreover, we checked that, if we try to capture the supertetragonal c/a, this makes it difficult to obtain a correct prediction for the ground state, as the P[001] state tends to become dominant.

B. La0.25Bi0.75FeO3
Now we discuss the case of La 0.25 Bi 0.75 FeO 3 . We first compute the parameters of the potential using the ana-    Table III. Next, we compare the model predictions and DFT values for our considered polymorphs in Fig. 4 (these results are also summarized in Table S2 of the Supplementary Material, together with the corresponding strains). Let us first consider the polymorphs with fixed cubic cell ( Fig. 4(a)-(c)). Our potential provides very accurate predictions for polarizations and tilts, similarly to the case of pure BiFeO 3 . The energy and relative stability of these polymorphs is also well captured by the model. Indeed, for polar-only structures both the model and DFT predict the rhombohedral P[111]c state to be lower in energy than the tetragonal P[001]c phase. The same holds for the polymorphs having only FeO 6 rotations: R[111]c is lower in energy than R[001]c. Note that the energy difference between the structures with tetragonal and rhombohedral phases are reduced compared to the case of pure BiFeO 3 . The lowest-energy phase is P[111]+R[111]c, where polarization and tilts coexist.
For the polymorphs in which shape and volume of the cell are allowed to relax, we observe the following. First, the model predicts accurate values of the polarization in all cases except for P[001]. Indeed, for the supertetragonal state, the predicted P s and η zz are underestimated relative to the DFT values. This issue is also reflected in the energy of the polymorph: our potential predicts P[001] to be the highest-energy state, while DFT shows that this phase is the second-lowest in energy, right above the P[111]+R[111] ground state. Additionally, we find that the tilts are accurately predicted by our model for all polymorphs except R[001]; in that case, the tilt amplitude and the strains are exaggerated compared to the DFT values. Note that, as it was the case for pure BiFeO 3 , these deficiencies are the result of the limited amount of DFT information on states P[001] and R[001] that was used to derive the parameters of our model.

C. Intermediate compositions
In this section we demonstrate how our potentials can be used to study La 1−x Bi x FeO 3 with intermediate La content, 0 < x < 0.25. We focus on the case of La 0.125 Bi 0.875 FeO 3 and check whether the properties of the polymorphs from the training set can be predicted by linear interpolation between BiFeO 3 and La 0.25 Bi 0.75 FeO 3 .
We consider two types of interpolation. First, we construct a model for x = 0.125 with coefficients obtained from interpolation of the corresponding values for the x = 0 and x = 0.25 cases. Using this model, we can easily predict the properties (P s , R s and E s ) of all the polymorhps in the training set. Second, we derive the very same properties by direct interpolation of the values obtained at x = 0 and x = 0.25. In Fig. 5 we compare the quantities thus obtained, and also include the corresponding DFT values for the polymorphs for which the information is available (see figure caption). We find that both interpolation approaches yield very similar preditions. Further, the agreement with DFT is good except for the supertetragonal P[001] phase, where our predictions suffer from the issues discussed above. Hence, we conclude that our models give us a way to treat compounds with intermediate compositions.

V. DISCUSSION
Let us now discuss the physical insights that our models provide.

A. P-R coupling
As it is well known from both experiments and computations, and correctly captured by our models, the ground state of BiFeO 3 has rhombohedral symmetry with P [111] and R [111]. It is interesting to note, though, that the DFT energy of the polar-only BiFeO 3 polymorph P[001] is lower than that of P[111]. By contrast, among the polymorphs having only FeO 6 octahedral tilts, R[111] is the lowest-energy structure. These observations yield one important conclusion: that the rhombohedral symmetry of the BiFeO 3 's ground state critically depends on the presence of the octahedral tilts, as in their absence the material would be tetragonal.
In order to understand how the rhombohedral ground state of BiFeO 3 comes about, we consider the F (P, R) part of our potential (Eq. 5) describing the coupling between polarization and octahedral rotations. The second term in Eq. (5), with C P R < 0 (see Table III), favors states where P and R are along/about any 111 direction, as for example, P [111] and R [111]. In turn, the third term, with C ′ P R < 0, favors phases where the FeO 6 tilts are about the axis defined by the polarization. Overall, these couplings lead P and R to appear together and aligned along/about the same 111 direction. Thus, these are the interactions driving the stabilization of the ground state phase of BiFeO 3 , rhombohedral and with co-existing polarization and tilts.
Does this mean, however, that P and R cooperate in BiFeO 3 ? We address this question by considering the energy diagram presented in Fig. 6 Table IV). We also show the energy of a virtual state in which P [111] coexists with R [111] but where these order parameters are not coupled. The energy of this virtual state is simply given by E 2c + E 4c , taking the cubic phase as the zero of energy. Clearly, the P[111]+R[111]c polymorph is higher in energy than the non-interacting virtual state and their energy difference arises from the coupling between P and R in P[111]+R[111]c structure. To understand this, we again consider the term F (P, R) (Eq. 5) with the corresponding coefficients B P R , C P R and C ′ P R presented in Table  III for BiFeO 3 . For the polymorph P[111]+R[111]c, we have P x = P y = P z = P 6c and R x = R y = R z = R 6c ; therefore, F 6c (P, R) = 3(3B P R + C P R + C ′ P R )P 2 6c R 2 6c for this state. In this expression, 3B P R > 0 dominates over C P R + C ′ P R < 0 and leads to a ground state energy that is higher than that of the virtual non-interacting state.
Thus, we find that, overall, the P and R order parameters compete in BiFeO 3 (B P R > 0 dominates). Nevertheless, the polar and tilt instabilities are so strong that this repulsive interaction is not enough to prevent them from occurring simultaneously. Further, the P-R competition is minimized when the order parameters are oriented along/about the same 111 axis (C P R , C ′ P R < 0), which yields the rhombohedral ground state phase of BiFeO 3 .

B. Effects of La doping
Let us first consider how La doping affects the electric polarization. From Figs. 4(a) and 4(d) one can see that, for all considered polar polymorphs in the training set, a 25% La doping leads to reduction of P . Indeed, for the polymorphs with fixed cubic cell ( Fig. 4(a)) we obtain a reduction of P by 11 − 19% for P[001]c, P[111]c and P[111]+R[111]c, and an even larger reduction for P[001]+R[001]c (≈ 59%). When we allow the cell to relax ( Fig. 4(d)), the obtained P reduction is in the range of 5 − 20%.
Next, let us turn to the effect of La doping on the FeO 6 octahedral tilts. As one can see from Figs. 4(b) and 4(e), the presence of 25% La has a relatively small effect (reduction) in the amplitude of tilts. More precisely, we find that the R Our models allow us to rationalize the most important results described above. Let us start by noting that the P-R couplings (B P R , C P R and C ′ P R in F (P, R)) are not significantly affected by the doping. Hence, they do not play a significant role to explain the La-induced effects. Table IV. Energies Es of BiFeO3 polymorphs calculated using DFT and predicted by the Landau-like potential introduced in this work (the coefficents of the potential are obtained using the analytical approach described in Sec. III B 2). Energy values are relative to the energy of the reference cubic phase and given in eV per 5-atom unit cell.  doping, indicating a weaker ferroelectric instability of the cubic phase. Additionally, both B P and C P increase and the relevant combination, 3B P + C P > 0, becomes larger; hence, the quartic couplings have a stronger effect on the energy landscape compared to pure BiFeO 3 . All these changes cooperate to yield shallower ferroelectric energy wells associated to F (P ) for La-doped BiFeO 3 , with smaller equlibrium polarization and lower energy barrier between states of opposite P. Note that this is consistent with previous studies on the effect of Ladoping on the switching characteristics of BiFeO 3 24,28 . As regards the tilt energy given by F (R), Table III shows that the presence of La weakens the cubic-phase instability (A R becomes less negative); by contrast, the quartic term (3B R + C R > 0) gets reduced upon doping, thus favoring larger tilts. These changes oppose each other, and result in the generally observed moderate reduction in the amplitude of the FeO 6 rotations.

Polymorph
Finally, as shown in Fig. 4, the P[001]+R[001]c case is peculiar, as it presents the largest reduction in P (about 59 %) and is the only one displaying an increase of R (about 12 %). We can rationalize this result by noting that, for this state, the quartic part of the energy in F (P ) (resp. F (R)) is controlled by the B P (resp. B R ) coupling alone. Upon doping B P grows (B R decreases), which favors smaller polarizations (larger tilts). Further, because of the strong competition between polarization and tilts in tetragonal states (B P R > 0; C P R and C ′ P R do not contribute), the changes get particularly large in the case of P[001]+R[001]c. Note also a subtle difference between P[001]+R[001]c and P[111]+R[111]c. In the latter case, the relevant quartic parameter for the tilts is 3B R + C R , and the La-induced decrease in B R is partly compensated by the increase in C R ; as a result, the tilts do not grow at all (recall A R < 0 grows upon doping) and the decrease of the polarization is relatively small.
Note that all these observations are consistent with what we know about the atomistic origin of the polar and tilt instabilities in BiFeO 3 . The former rely on the presence of stereochemically active 6s lone pairs in the Bi 3+ cations; hence, their partial substitution by lonepair-free La cations naturally leads to smaller polarizations. The latter are mainly controlled by the ionic radius of the Bi 3+ cation; since La 3+ is similar in size, the doping leaves R largely unaffected.

VI. CONCLUSIONS
In summary, we have introduced the simplest, lowestorder Landau-like potential for BiFeO 3 and related compounds, as well as methods that allow to compute the potential parameters from Density Functional Theory (DFT). More precisely, we have derived analytical expressions for all the model coefficients as functions of the energies and structural features (polarization, FeO 6 octahedral tilts and strains) of a small set of relevant polymorphs. We have applied the proposed approach to BiFeO 3 and La 0.25 Bi 0.75 FeO 3 , showing its overall accuracy in reproducing the DFT data. We have also showed that our models can be used -by interpolation -to predict the properties of compounds with intermediate dopant concentrations. We note that the introduced potential, as well as the analytical scheme to obtain its coefficients from DFT, can be readily applied to study the properties of other perovskite oxides characterized by the same order parameters (polarization, antiphase oxygen-octahedral tilts, strains). This includes ferroelectrics where the tilts are not important (e.g., BaTiO 3 or PbTiO 3 ), antiferrodistortive non-polar perovskites (e.g., LaAlO 3 ), or compounds where both distortions play a relevant role (e.g., SrTiO 3 ), as well as their corresponding solid solutions. In principle, an extension of our scheme to compounds where other order parameters are relevant (e.g., in-phase tilts in orthorhombic perovskites like CaTiO 3 45 ) should be straightforward.