Oxygen-vacancy induced magnetic phase transitions in multiferroic thin films

Multiferroics in which giant ferroelectric polarization and magnetism coexist are of tremendous potential for engineering disruptive applications in information storage and energy conversion. Yet the functional properties of multiferroics are thought to be affected detrimentally by the presence of point defects, which may be abundant due to the volatile nature of some constituent atoms and high temperatures involved in materials preparation. Here, we demonstrate with theoretical methods that oxygen vacancies may enhance the functionality of multiferroics by radically changing their magnetic interactions in thin films. Specifically, oxygen vacancies may restore missing magnetic super-exchange interactions in large axial ratio phases, leading to full antiferromagnetic spin ordering, and induce the stabilization of ferrimagnetic states with a significant net magnetization of 0.5 uB per formula unit. Our theoretical study should help to clarify the origins of long-standing controversies in bismuth ferrite and improve the design of technological applications based on multiferroics.

Multiferroics in which giant ferroelectric polarization and magnetism coexist are of tremendous potential for engineering disruptive applications in information storage and energy conversion. Yet the functional properties of multiferroics are thought to be affected detrimentally by the presence of point defects, which may be abundant due to the volatile nature of some constituent atoms and high temperatures involved in materials preparation. Here, we demonstrate with theoretical methods that oxygen vacancies may enhance the functionality of multiferroics by radically changing their magnetic interactions in thin films.
Specifically, oxygen vacancies may restore missing magnetic super-exchange interactions in large axial ratio phases, leading to full antiferromagnetic spin ordering, and induce the stabilization of ferrimagnetic states with a significant net magnetization of 0.5 µ B per formula unit. Our theoretical study should help to clarify the origins of long-standing controversies in bismuth ferrite and improve the design of technological applications based on multiferroics.
Finding multiferroics in which ferroelectricity and magnetism coexist and influence each other is of great fundamental and applied interests [1,2]. Salient technological features of multiferroics include the possibility of controlling the magnetization with electric fields to design efficient logic and memory devices [3,4], and of realizing large piezomagnetic coefficients to facilitate the miniaturization of antennas and sensors [5,6]. Furthermore, competition between phases displaying distinct electric polarization and magnetic ordering offers also encouraging prospects for energy conversion applications like photovoltaics and solid-state cooling [7][8][9].
Through epitaxial strain is actually possible to create new multiferroic materials in the laboratory that ex- * Corresponding Author hibit giant electric polarization and unexpected magnetic spin ordering [19][20][21]. An illustrative example is given by BiFeO 3 (BFO), in which large spontaneous polarization and ferromagnetism (FM) have been observed under moderate compressive biaxial strains at room temperature [22,23]. The origins of the net magnetization in BFO thin films, however, are not clear yet and from a technological point of view is crucially important to understand them at the fundamental level.
Most magnetic ferroelectrics with chemical formula AM O 3 and perovskite-like structure present antiferromagnetic spin ordering along the three pseudo-Cartesian directions (AFM-G), due to the dominant role of oxygenmediated super-exchange interactions between neighbouring transition metal atoms M [24][25][26][27]. In large axial ratio structures, the covalency of M -O bonds parallel to the electric polarization is significantly reduced and consequently magnetic exchange interactions, favouring parallel magnetic spins, dominate in that direction; the coexistence of "in-plane" antiferromagnetism and "outof-plane" ferromagnetism leads anyway to null crystal magnetization (AFM-C, Fig.1a,b) when small spin canting effects are neglected. Therefore, intrinsic and robust FM in principle is not expected to occur in BFO or any other similar multiferroic [28][29][30][31].
A plausible explanation for the appearance of FM in AM O 3 perovskite oxide thin films is based on extrinsic causes like point defects [23,32]. The volatile nature of bismuth and high temperatures involved in the preparation of samples, for instance, make the presence of oxygen vacancies (V O ) almost inevitable in Bi-based multiferroics [15,17]. In fact, oxygen defects may modify significantly the structural and functional properties of perovskite thin films via changes in the M oxidation states and their coupling with the lattice strain [33,34]. However, a number of theoretical works based on firstprinciples methods have agreed in that the combined action of V O and lattice strain may not affect considerably the magnetic properties of Bi-based multiferroics [35][36][37].
Here, we present new theoretical evidence showing that the presence of V O may in fact change radically the magnetic properties of multiferroic thin films via previously overlooked electro-structural mechanisms. We select BiCoO 3 (BCO) as the model multiferroic in which to perform first-principles calculations based on density functional theory (DFT) because (i) this material already exhibits a large axial ratio in the absence of any strain, and (ii) the magnetic effects that we predict can be realized on substrates that are commonly employed for growth of epitaxial oxide perovskite thin films (Fig.1c). In particular, it is found that oxygen vacancies occupying specific lattice positions can induce the stabilization FIG. 1. Sketch of multiferroic BCO in bulk and thin film geometries. a Representation of the super-tetragonal phase (T , space group P 4mm) characteristic of bulk BCO and other multiferroics exhibiting giant electric polarization. b Magnetic structure in the T phase rendering C-type antiferromagnetism (AFM-C); green arrows represent atomic magnetic moments and their orientation. c Examples of substrates in which to grow BCO thin films displaying the effects predicted in this study; other well-known materials exhibiting super-tetragonal phases are shown for comparison [19,22]. "SE" stands for magnetic super-exchange interactions. Bi, Co, and O atoms are represented with magenta, blue, and red spheres, respectively. of full antiferromagnetic (AFM-G) super-tetragonal and ferrimagnetic (FiM) monoclinic polar phases, depending on the lattice strain. As a consequence, phase competition is enriched and magnetic functionalities further enhanced in comparison to perfectly stoichiometric thin films. We show that most of the results obtained in BCO thin films can be generalized to BFO and other Bi-based multiferroics, hence our conclusions are of broad applicability and significance to the field of functional materials.

RESULTS
V O -induced effects on phase competition and functionality. Bulk BCO presents a polar tetragonal T phase with a large axial ratio of c/a ≈ 1.3 and relatively small lattice paramater a = 3.76Å [9,29] (Fig.1a). The competing structures are a non-polar orthorhombic O phase and a polar monoclinic M phase (Supplementary Fig.1); both competing phases have cells that are slightly distorted versions of the ideal cubic perovskite structure with c/a ≈ 1. The polar phases in BCO present spontaneous polarizations along quite different crystallographic directions, namely, pseudocubic [001] pc in T and ∼ [111] pc in M. As regards magnetism, the O and M phases exhibit G-type antiferromagnetism (AFM-G) with a quite high Néel temperature, T N ≈ 500 K, whereas the T phase C-type antiferromagnetism (AFM-C) with a relatively low T N of ≈ 310 K. In stoichiometric BCO thin films, and by completely neglecting temperature effects, a multiferroic T → M phase transition involving large structural, polar, and magnetic changes occurs at in-plane parameter a in = 3.91Å [9] (Figs.2a-c).  Supplementary Fig.2). For the smallest in-plane lattice parameters, a T -C phase (magnetic spin ordering is indicated along with the structure symmetry) containing oxygen vacancies in equatorial (Eq) positions ( Fig.1a) renders the lowest energy. The electric polarization and Néel temperature in T -C(Eq) are significantly lower than in the analogous stoichiometric phase, in particular, we estimate differences of ∆P ≈ −75 µC cm −2 and ∆T N ≈ −50 K for same in-plane parameters. At a in = 3.77Å, an unusual magnetic phase transition from AFM-C to AFM-G spin ordering occurs along with the appearance of a small in-plane electric polarization (P xy ∼ 10 µC cm −2 ) and change in V O position symmetry. The Néel temperature in the T -G(Ap) phase is lower than in T -C(Eq) by approximately 50 K. Furthermore, at a in ≥ 3.91Å the system adopts a monoclinic ferrimagnetic (FiM) phase with oxygen vacancies in equatorial positions, M-FiM(Eq), and a considerable net magnetization of ≈ 0.5 µ B per formula unit (Fig.2d). The Néel temperature in the M-FiM(Eq) phase is larger than in T -G(Ap) and remains close to room temperature almost independently of a in . ( The physical mechanisms responsible for the two multiferroic phase transitions represented in Fig.2d will be explained in detail in the next subsections. Let us now comment briefly on the functionality enhancement deriving from the T -G(Ap) → M-FiM(Eq) transformation by keeping in mind that the non-stoichiometric M phase is ferrigmagnetic and polar hence responsive to both external magnetic and electric fields. First, a large change in the electric polarization orientation involving a ≈ 60 • rotation is observed during the transition FIG. 2. Effects of VO on phase competition and functionality in BCO thin films. a Zero-temperature energy of competing phases expressed as a function of in-plane lattice parameter. Metastable phases are indicated by grey curves and strain-induced phase transitions by vertical dashed lines; the oxygen vacancy positions leading to lowest energies, either apical "Ap" or equatorial "Eq", are indicated within parentheses. "G" stands for G-type antiferromagnetism, "C" for C-type antiferromagnetism, and "FiM" for ferrimagnetism. b Electric polarization of stoichiometric and non-stoichiometric ground-state phases. c Magnetic transition temperature of stoichiometric and non-stoichiometric ground-state phases. d Phase transition sequence occurring in non-stoichiometric BCO thin films under increasing ain; the green arrows represent the magnetic spin ordering in each phase. e Change in volume, ∆V , change in electric polarization orientation, ∆α, and change in magnetic moment per formula unit, Mi, associated to the multiferroic T -G → M-FiM phase transition. (Fig.2e); as a consequence, and in analogy to what has been observed in Pb(Zr 1−x Ti x )O 3 alloys [38,39] and Bi(Fe 1−x Co x )O 3 thin films [40], it should be possible to realize large piezoelectric responses under small electric bias at a in ≈ 3.91Å. Second, the sizeable changes in electric polarization and total magnetization in principle should allow for control of the polarization with magnetic fields and vice versa, which hints at the likely existence of large magnetoelectric couplings [1, 2]. And third, the out-of-plane lattice parameter shrinks by an impressive ≈ 11% (Fig.2e) hence there is the possibility of realizing giant piezomagnetic responses [5,6] and multicaloric effects [8,9] through the application of external bias. In a more speculative vein, the change in V O position symmetry from Eq to Ap could lead to novel ionic transport phenomena driven by external magnetic, rather than electric, fields [41]. As we will show later, similar magnetic phenomena are likely to occur also in FIG. 3. Electronic, structural, and magnetic properties of T BCO thin films. a Stoichiometric T phase. b Non-stoichiometric T thin films with a stable VO in Eq position and AFM-C spin ordering. The orange plane indicates the two Co ions that are reduced as a consequence of creating a neutral Eq oxygen vacancy. The red arrow in the d-orbitals occupation sketch indicates the difference with respect to the stoichiometric case. c Non-stoichiometric T thin films with a stable VO in Ap position and AFM-G spin ordering. The orange plane indicates the two Co ions that are reduced as a consequence of creating a neutral Ap oxygen vacancy. d Charge density surface plot corresponding to the stoichiometric T -C phase; the surface over which the charge density is calculated is indicated by a grey plane in the accompanying structural ball-stick representation. Isovalue paths are represented with black and coloured lines. e Charge density surface plot corresponding to the non-stoichiometric T -G(Ap) phase. Regions of interest describing Co-O bonds are indicated with yellow arrows.
V O -induced magnetic super-exchange interactions in the T phase. Figure 3 summarizes the electronic, structural, and magnetic properties of stoichiometric and non-stoichiometric T BCO thin films. In the stoichiometric T -C phase (Fig.3a), the square-pyramidal O 5 crystal field splits the electronic Co d levels into nondegenerate b 2g (d xy ), doubly degenerate e g (d xz , d yz ), and nondegenerate a 1g (d z 2 ) and b 1g (d x 2 −y 2 ). Our first-principles calculations render a high-spin Co state characterised by the electronic occupation b 2 2g e 2 g a 1 1g b 1 1g and atomic spin moment 3.1 µ B , in good agreement with the available experimental data [42]. In the nonstoichiometric T -C(Eq) phase (Fig.3b), the splitting of electronic d levels remains invariant with respect to the stoichiometric case and the occupation in the two cobalt ions nearest to the neutral V O , which become reduced and are electronically equivalent, changes slightly to b 2 2g e 2 g a 2 1g b 1 1g (depending on the choice of the tech-nical DFT parameters this electronic distribution may vary somewhat - Supplementary Fig.5 and Supplementary Methods-). Interestingly, when V O is created in an apical position and for specific a in 's the magnetic spin ordering in the T phase changes to AFM-G. Figure 3c shows the electronic density of states of the two reduced cobalt ions in the T -G(Ap) phase, which in this case turn out to be electronically inequivalent. In particular, the doubly degenerate e g (d xz , d yz ) orbitals in the cobalt ion closest to the apical V O (Co3, as labelled in Figs.3d,e) undergo a significant energy reduction and become fully populated rendering the occupation state e 4 g b 1 2g a 1 1g b 1 1g , while the other reduced metal ion (Co1, as labelled in Figs.3d,e) exhibits the more usual distribution b 2 2g e 2 g a 2 1g b 1 1g . These drastic electronic rearrangements are correlated with the appearance of a strong structural distortion in the system that pushes Co3 towards the oxygen atom underneath of it (see inversion of the corresponding O 5 square-pyramid in Figs.3c,e) and tends to restore (partially) the missing magnetic super-exchange interactions along the out-ofplane direction. This super-exchange restoration mechanism, which is accompanied by an increase in covalency of the out-of-plane Co1-O-Co3 bonds and eventually leads to the stabilization of AFM-G spin ordering, is clearly imaged by plots of the electronic density in the plane containing Co1 and Co3 and oriented perpendicular to the substrate (Figs.3d,e and yellow arrows therein).
The V O -induced AFM-C → AFM-G phase transition disclosed in large axial ratio BCO thin films may shed some light on uncomprehended experimental observations of antiferromagnetic spin ordering in other super-tetragonal multiferroic phases. For instance, in T BFO thin films several first-principles works have predicted AFM-C spin ordering [43][44][45] whereas most experimental studies indicate that AFM-G dominates [46,47]. As it will be explicitly shown later, by considering the presence of oxygen vacancies in BFO thin films those theoretical and experimental results may be reconciled.
V O -induced stabilization of a ferrimagnetic M phase. Figure 4 summarizes the electronic, structural, and magnetic properties of stoichiometric and nonstoichiometric M BCO thin films. In the stoichiometric M-G phase (Fig.4a), the octahedral O 6 crystal field splits the electronic Co d levels into doubly degenerate e g (d z 2 , d x 2 −y 2 ) and triply degenerate t 2g (d xy , d yz , d xz ); a strong Jahn-Teller distortion rendering a large Q 2 value of 0.37Å [48] lifts further the degeneracy in the e g and t 2g manifolds [49] and promotes the electronic occupation state b 2 2g e 2 g b 1 1g a 1 1g (Fig.4a). Remarkably, when specific V O 's are created in equatorial positions (see next paragraph) the lowest-energy magnetic spin ordering changes to FiM and the net magnetization per formula unit in the M-FiM(Eq) phase amounts to ≈ 0.5 µ B . The two Co ions that are reduced by the neutral vacancy present same magnetic moment orientation, same electronic oc-cupancy e 4 g b 1 2g b 1 1g a 1 1g , and sit within the [111] pc plane (Fig.4b).
How is possible that the two reduced Co ions are located along the diagonal of the pseudo-cubic unit cell rather than within the equatorial plane (that is, closest to the neutral V O , in which case the total magnetization would be null)? The ground-state M-FiM(Eq) phase appears only when highly magnetized oxygen atoms (0.1-0.2 µ B ) occupying equatorial positions in the stoichiometric M-G crystal are removed (Fig.4c). In that case, as we sketch in Fig.4d, is not possible to reduce two neighbouring Co ions sitting in the equatorial plane (Co1 and Co2) due to Pauli exclusion principle. Consequently, pairs of metal ions with same magnetic moment and orientation (Co1 and Co4 in Fig.4d) become reduced. For the couple of distant Co1 and Co4 ions to change their oxidation state and magnetic moment, however, the crystal needs to undergo sizable structural distortions involving the Bi and O atoms surrounding V O (Supplementary Tables 3-4). In particular, the non-magnetic oxygen atom in apical position just above Co2 acts as a bridge between the equatorial oxygen vacancy and Co4, by lending one of its electrons to the metal ion, hence reducing it, and receiving one electron from V O (Fig.4d). This concerted electronic hopping mechanism can be inferred from plots of the electronic density in the plane containing Co1, Co2, Co3, and Co4 and oriented perpendicular to the substrate (Figs.4e,f and yellow arrows therein). As can be appreciated therein, the covalency of the Co4-O bond is significantly reduced in the M-FiM(Eq) phase as compared to the stoichiometric case, and a vertical tilt of the Co2-O bond that brings the apical oxygen closer to the equatorial V O is also evidenced.
The discovery of FiM spin ordering in the nonstoichiometric M phase motivated us to search for similar magnetic states, even if metastable, in the other BCO thin film geometries T and O. The presence of highly magnetized O atoms was acknowledged in both stoichiometric phases ( Supplementary Fig.6), however upon removal of those oxygens we did not observe the appearance of any net magnetic moment (neglecting small spin canting effects). These results confirm the importance of the electro-structural mechanisms just described on facilitating the stabilization of FiM spin ordering. For instance, in the T phase the highly magnetized O atoms appear in apical positions ( Supplementary Fig.6) hence the nonmagnetic oxygens, which occupy equatorial positions and are somewhat clamped to the substrate, cannot act as electronic bridges between distant V O 's and Co's. Meanwhile, in the O phase, which arguably is quite similar to M in terms of structure and magnetism (Supplementary Fig.1) [9,29], the lack of polar order and high dielectric permittivity makes the enabling Bi-O structural distortions (Supplementary Table 4) to be too high in energy; consequently, FiM spin ordering is frustrated. Based on these outcomes, we propose that the following three conditions are necessary for the stabilization of ide perovskites: (1) lack of inversion symmetry leading to polar order and structural deformation ease, (2) moderate axial ratio structures with c/a ≈ 1 allowing for out-of-plane concerted electronic hoppings, and (3) the existence of highly magnetized oxygen ions. cial case of BiFeO 3 (BFO) and other Bi-based multiferroic (BiMnO 3 and BiCrO 3 ) thin films. Figure 5 encloses the energy, magnetic, and vibrational properties of nonstoichiometric BFO thin films. In order to be consistent with the notation employed heretofore, we label the usual rhombohedral-like monoclinic phase as M and the large axial ratio tetragonal-like phase as T (in spite of the fact that the space groups corresponding to those structures are different in BFO and BCO) [45,50]. Our zero-temperature calculations (Fig.5a) predict a ground-state M-G(Ap) phase at a in ≥ 4.01Å, followed by T -G(Eq) at 3.91 ≤ a in ≤ 4.01Å, and T -C(Eq) at a in ≤ 3.91Å. Indeed, when the likely existence of V O is explicitly considered in the simulations a broad a in region appears in which AFM-G spin ordering is stable in the T phase (that is missing in the corresponding stoichiometric system [43,45]). The causes of the stabilization of the T -G(Ap) phase in BFO are very similar to those explained previously for T -G(Ap) in BCO thin films ( Fig.3; we note that Ap V O 's under tensile strain are in some ways equivalent to Eq V O 's under compressive strain). As mentioned earlier, these results may shed new light on the origins of some unresolved discrepancies between theory and experiments regarding the determination of antiferromagnetic ordering in T BFO thin films [43][44][45][46][47].
By creating an apical V O in the stoichiometric M-C phase, we found a thus far neglected ferrimagnetic phase in monoclinic BFO thin films, M-FiM(Ap). In this case the net magnetization per formula unit amounts also to ≈ 0.5 µ B . Nevertheless, the M-FiM(Ap) phase is metastable at zero temperature due to a small energy difference of 25 − 30 meV per formula unit with respect to the ground-state phase (Fig.5a). The atomic structure of the metastable M-FiM(Ap) and ground-state M-G(Ap) phases are highly distorted and surprisingly very similar ( Fig.5b and Supplementary Tables 5-6). Analysis of the Γ-point phonon modes (Fig.5c), however, indicates that the M-FiM(Ap) phase is vibrationally softer than M-G(Ap). Consequently, the M-FiM(Ap) phase may be entropically stabilized over M-G(Ap) under increasing temperature since the zero-temperature energy difference between the two states is relatively small and the vibrational free energy of the former phase is more favourable [9,29]. Our theoretical results, therefore, can be interpreted as evidence showing that the observation of "ferromagnetic" behaviour in BFO thin films [22] may be caused by the presence of oxygen vacancies, just as it has been suggested by other authors [23]. Figure 5d shows the spin-up spin-down charge densities calculated in non-stoichiometric M and T BFO thin films considering in-plane and out-of-plane surfaces. The reason why we could find just one FiM solution through the generation of neutral V O 's in stoichiometric M phases is now clear [see requirements (1)-(3) listed at the end of the previous section]: all oxygen atoms in the M-G phase are non-magnetic whereas all apical O in the M-C phase are highly magnetized. The stoichiometric T -C phase also displays highly magnetized oxygen atoms in apical positions however, as we have explained before, FiM spin ordering hardly can be generated in geometries exhibiting c/a 1. Furthermore, we repeated the same type of calculations and analysis in monoclinic-like BiMnO 3 and BiCrO 3 thin films, which fulfill conditions (1) and (2) explained above, and found that the appearance of FiM spin ordering is also correlated with the presence of highly magnetized oxygen atoms (Supplementary Fig.7). In particular, a M-FiM(Ap) phase displaying a net magnetization of ≈ 0.2 µ B is found in BiCrO 3 . These results confirm that the simple rules provided in this work for prediction of FiM phases in non-stoichiometric multiferroic thin films are robust and general.
In summary, by using first-principles calculations we have disclosed a number of previously overlooked electrostructural mechanisms induced by the presence of oxygen vacancies that facilitate the stabilization of unexpected magnetic states in multiferroic thin films. In particular, AFM-G and FiM spin orderings may naturally appear in large axial ratio and monoclinic phases under certain lattice strain conditions. Our theoretical results may clarify the origins of some long-standing controversies in BFO, the paradigm of single-phase multiferroics and one of the most intensively studied functional materials. We provide general and simple rules to fundamentally understand and systematically predict FiM phases in non-stoichiometric multiferroics, thus offering new approaches for the rational engineering of bettered functional materials. The present work shows that oxygen vacancies should be considered as a design opportunity to create new funcionalities, especially as related to magnetism, in multiferroic thin films.

METHODS
Density functional theory calculations. Firstprinciples spin-polarized calculations based on density functional theory (DFT) are performed with the generalized gradient approximation proposed by Perdew, Burke and Ernzerhof (GGA-PBE) as implemented in the VASP package [51,52]. We employ the "Hubbard-U" scheme derived by Dudarev et al. to deal with the 3d electrons in Co (Fe) atoms and, as done in previous works, a U value of 6 (4) eV is adopted [9,28,29]. We use the "projected augmented wave" method [53] considering the following electronic states as valence: Co's 4s 1 3d 8 , Fe's 3p 6 4s 1 3d 7 , Mn's 4s 1 3d 6 , Cr's 4s 1 3d 5 , Bi's 6s 2 5d 10 6p 3 , and O's 2s 2 2p 4 . The energy cut-off is truncated at 650 eV and we employ a Γ-centered k-point grid of 4 × 6 × 6 for a 2 × √ 2 × √ 2 supercell containing 20 atoms (that is, four formula units) [50]. Periodic boundary conditions are applied along the three lattice-vector directions. Thin film geometry relaxations are carried out by using a conjugated gradient algorithm that allows to change the simulation-cell volume and atomic positions while constraining the length and orientation of the two in-plane lattice vectors. The geometry relaxations are stopped once the forces acting on the ions are smaller than 0.01 eV/Å. We have checked the vibrational stability of every phase by estimating the lattice phonons at the Γ-point with the small-displacement method [54] and considering central differences for the calculation of atomic forces derivatives.
Non-stoichiometric phases are generated by removing one oxygen atom from an apical or equatorial position in the 20-atoms simulation cell, thus rendering the chemical composition BiCoO 3−x with x = 0.25. Apical and equatorial V O positions are investigated systematically by generating all inequivalent configurations in all competing phases and considering all possible magnetic spin orderings (FM, AFM-G, AFM-G, and AFM-A [29]). The results presented in the main text are obtained by assuming neutral oxygen vacancies since we have found that neutral V O 's are energetically more favourable than charged oxygen vacancies ( Supplementary Fig.8, Supplementary Table 7, and Supplementary Methods). In particular, we have employed the following well-established formula to estimate the ranking of V O formation energies as a function of charge, q, [55]: where E[V q O ] is the energy of the non-stoichiometric system containing the oxygen vacancy, E stoi the energy of the corresponding stoichiometric system, E q corr a finitesize supercell correction, n O the number of created V O 's (typically equal to 1 in our calculations), µ O the chemical potential of oxygen atoms, F the Fermi energy in the non-stoichiometric system, v the top energy in the valence band of the non-stoichiometric system, and ∆V a term used for aligning the electrostatic potentials of the stoichiometric and defective supercells. In order to calculate E q corr and ∆V , we have followed the methods explained in work [55]. According to our estimations, neutral V O 's (q = 0) are energetically more favourable than charged vacancies (q = +2 e) by about ∼ 1 eV per formula unit ( Supplementary Fig.8, Supplementary Table 7, and Supplementary Methods).
We have performed several tests to assess the influence of the adopted DFT exchange-correlation functional and U value on our theoretical predictions (Supplementary Figs.3-5 and Supplementary Methods). Specifically, we repeated most calculations by considering the PBEsol functional [56] and 2 ≤ U ≤ 6 eV values. It is found that the main conclusions presented in the main text are not affected qualitatively by the choice of the U parameter or exchange-correlation functional. At the quantitative level, the a in parameters at which the phase transitions occur and the electronic occupations that are deduced from electronic density plots change a little in some cases (Supplementary Figs.3-5 and Supplementary Methods); however, the structural properties and energy ranking of the competing phases estimated in most a in cases remain invariant. For a detailed discussion on these technical aspects, see Supplementary Methods.
Regarding the estimation of the electric polarization, P , we started by employing the Berry phase formalism [57]. However, the presence of oxygen vacancies induces a notable reduction in the energy band gap of the system ( Supplementary Fig.2) that in some cases frustrates the determination of the corresponding Berry phase (due to the appearance of intermediate metallic phases). In order to overcome such a limitation, we opted for calculating the electric polarization perturbatively. Specifically, we estimate P with the formula [50]: where Ω is the volume of the cell, κ runs over all the atoms, α, β = x, y, z represent Cartesian directions, u κ is the displacement vector of the κ-th atom as referred to a non-polar reference phase, and Z * κ the Born effective charge tensor calculated in the non-polar reference state. In the stoichiometric phases we do not find the technical limitations just explained for non-stoichiometric systems, hence in that case we have been able to compare the P values obtained with the Berry phase approach (exact) and Eq.(2) (approximate). According to our estimations, the electric polarizations calculated perturbatively are accurate to within ∼ 10% of, and systematically larger than, the P values calculated with the Berry phase method. It is reasonable to assume then a similar level of accuracy in the P values estimated in non-stoichiometric thin films that are reported in Fig.2.
Heisenberg model Monte Carlo simulations. To simulate the effects of thermal excitations on magnetic ordering in T , O, and M BCO thin films, we construct several spin Heisenberg models of the formĤ = 1 2 ij J (0) ij S i S j , in which the value of the involved exchange constants are obtained from zero-temperature DFT calculations (see works [9,28,29] for the technical details on the determination of the J (0) ij parameters). We use those models to perform Monte Carlo (MC) simulations in a periodically-repeated simulation box of 20 × 20 × 20 spins; thermal averages are computed from runs of 50, 000 MC sweeps after equilibration. These simulations allow us to monitor the T -dependence of the magnetic ordering through the computation of the AFM-C (in the T phase) and AFM-G (in the O and M phases) order parameters, namely, S C ≡ 1 N i (−1) nix+niy S iz and S G ≡ 1 N i (−1) nix+niy+niz S iz . Here, n ix , n iy , and n iz are the three integers locating the i-th lattice cell, and N is the total number of spins in the simulation box. For the calculation of S C and S G , we considered only the z component of the spins because a small symmetry-breaking magnetic anisotropy was introduced in the Hamiltonian in order to facilitate the numerical analysis [9,28,29].

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author (C.C.) upon reasonable request.