{'fraction_non_alphanumeric': 0.03728523351242841, 'fraction_numerical': 0.01706323433902108, 'mean_word_length': 4.771852104276708, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
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.


INTRODUCTION
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 control of the magnetization with electric fields to design efficient logic and memory devices 3,4 , and the realization of 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][10] .
Through epitaxial strain, it is actually possible to create new multiferroic materials in the laboratory that exhibit giant electric polarization and unexpected magnetic spin ordering [21][22][23] . An illustrative example is given by BiFeO 3 (BFO) thin films, in which large spontaneous polarization and ferromagnetism (FM) have been observed under moderate compressive biaxial strains at room temperature 24,25 . The origins of the net magnetization reported for BFO thin films, however, are not clear yet and from a technological point of view it is crucial to understand them correctly.
Most magnetic ferroelectrics with chemical formula AMO 3 and perovskite-like structure present antiferromagnetic spin ordering along the three pseudo-Cartesian directions (AFM-G), due to the dominant role of oxygen-mediated super-exchange interactions between neighbouring transition metal atoms M [26][27][28][29] . In large axial ratio structures (c/a > 1), also known as "super-tetragonal", the covalency of M-O bonds oriented along the electric polarization direction is significantly reduced and consequently the magnetic exchange interactions, which favour parallel magnetic spins, dominate in that direction. The coexistence of "in-plane" antiferromagnetism and "out-of-plane" ferromagnetism in super-tetragonal phases (AFM-C, Fig. 1a, b) leads anyway to null crystal magnetization 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 [30][31][32][33] .
A plausible explanation for the appearance of FM in AMO 3 perovskite oxide thin films is based on extrinsic causes like point defects 25,34 . The volatile nature of bismuth and the high temperatures involved in the preparation of samples, for instance, make the presence of oxygen vacancies (V O ) almost inevitable in Bibased multiferroics 17,19 . 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 [35][36][37][38][39][40] . However, a number of theoretical works based on first-principles 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 [41][42][43] .
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, (ii) the magnetic effects that we predict could be realized on substrates that are commonly employed for growth of epitaxial perovskite thin films (Fig. 1c), and (iii) BCO naturally exhibits a high concentration of oxygen vacancies (i.e., it decomposes easily at temperatures around 700 K and needs to be synthesised at high pressures 13,14 ). In particular, it is found that oxygen vacancies occupying specific lattice positions can induce the stabilization 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 Bibased multiferroics, hence our conclusions are of broad applicability and significance to the field of multifunctional materials.

RESULTS
We have performed first-principles simulations of stoichiometric and non-stoichiometric BCO thin films based on density functional theory (DFT) techniques (see "Methods" section). The V O -induced effects explained in this section have been observed for compositions BiCoO 2.75 and BiCoO 2.875 . Detailed discussions on the influence of the selected DFT exchange-correlation energy functional and other technical parameters on our general conclusions can be found in the "Methods" and Supplementary Discussion sections. Overall, the results presented next are robust and qualitatively accurate.
V O -induced effects on phase competition and functionality It is known both from experiments and first-principles calculations that bulk BCO presents a polar tetragonal T phase (space group P4mm) with a large axial ratio of c/a ≈ 1.3 and relatively small lattice paramater a = 3.76 Å 9,14,31,44 (Fig. 1a). The competing structures are a non-polar orthorhombic O phase (space group Pnma) and a polar monoclinic M phase (space group Pc) (Supplementary Fig. 1, Supplementary Tables 1-3) 9,14,31,44 ; 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, recent firstprinciples calculations have shown that phases O and M exhibit G-type antiferromagnetism (AFM-G) with a quite high Néel temperature, T N ≈ 500 K, whereas the T phase presents Ctype antiferromagnetism (AFM-C) with a relatively low T N of ≈ 310 K 31 . In stoichiometric BCO thin films, and by completely neglecting temperature effects, a multiferroic T ! M phase transition involving large structural, polar, and magnetic changes has been predicted to occur at in-plane parameter a in = 3.91 Å 9 (Fig. 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 calculated for T -C(Eq) are significantly lower than for the analogous stoichiometric phase, in particular, we estimate differences of ΔP ≈ −75 μC cm −2 and ΔT N ≈ −50 K when considering identical in-plane parameters.
At a in = 3.77 Å, an unusual magnetic phase transition from AFM-C to AFM-G spin ordering occurs that is accompanied by the appearance of a small in-plane electric polarization (P xy~1 0 μC cm −2 ) and change in V O symmetry. The Néel temperature of the T -G(Ap) phase is lower than that of T -C(Eq) by approximately 50 K. Furthermore, at a in ≥ 3.91 Å the system adopts a monoclinic ferrimagnetic (FiM) phase with oxygen vacancies preferably in equatorial positions, M-FiM(Eq), and a considerable net magnetization (i.e., 0.5 μ B per formula unit for BiCoO 2.75 - Fig. 2d-and 0.25 μ B for BiCoO 2.875 ). The Néel temperature of the M-FiM(Eq) phase is larger than that of T -G(Ap) and remains close to room temperature almost independently of a in .
Interestingly, the formation energy of neutral oxygen vacancies in BCO thin films (see "Methods" section) is found to decrease steadily for a in ≥ 3.77 Å as compared to the bulk equilibrium system ( Supplementary Fig. 3). For instance, at a in = 3.95 Å the formation energy of the most stable type of V O is found to be 0.25 eV smaller than obtained at 3.76 Å (per formula unit). The difference in formation energy between Ap and Eq oxygen vacancies can be deduced directly from the energy curves shown in Fig. 2a (see "Methods" section). For instance, in the M phase the creation of an equatorial V O is~0.1 eV more favourable than an apical vacancy, almost independently of the in-plane lattice parameter.
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. First, a large change in the electric polarization orientation involving a rotation of ≈ 60 ∘ is observed during the transition (Fig. 2e); as a consequence, and in analogy to what has been observed in Pb(Zr 1−x Ti x )O 3 alloys 45,46 and Bi(Fe 1−x Co x )O 3 thin films 47 , it should be possible to realize large piezoelectric responses under small electric bias for 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 thin film dimension along the direction perpendicular to the epitaxial plane shrinks by an impressive ≈11% (Fig. 2e) hence there is the Fig. 1 Sketch of multiferroic BiCoO 3 (BCO). a Representation of the super-tetragonal phase (T , space group P4mm) 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 (bottom); other well-known materials exhibiting super-tetragonal phases are shown for comparison (top) 21,24 . "SE" stands for magnetic super-exchange interactions. Bi, Co, and O atoms are represented with magenta, blue, and red spheres, respectively.
possibility of realizing giant piezomagnetic responses 5,6 and multicaloric effects 8-10 through the application of external bias. As it will be shown later, similar magnetic phenomena than found for BCO are likely to occur also in other Bi-based multiferroic thin films, including BiFeO 3 .
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 squarepyramidal 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 magnetic moment of 3.1 μ B (in absolute value), in good agreement with the available experimental data 14 .
In the non-stoichiometric T -C(Eq) phase (Fig. 3b), the relative splitting of electronic d levels remains invariant with respect to the stoichiometric case and the electronic occupation of the two cobalt ions nearest to the neutral V O , which are reduced and are electronically equivalent (indicated with an orange plane in Fig.  3b), changes slightly to b 2 2g e 3 g a 1 1g b 1 1g . Interestingly, when V O is created in an apical position the magnetic spin ordering of the T phase changes to AFM-G at specific a in 's ( Fig. 2a). Figure 3c shows the electronic orbital occupation deduced from our DFT calculations for the two reduced cobalt ions in the T -G(Ap) phase (indicated with an orange plane in Fig. 3c), 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 Fig. 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 Fig. 3d, e) exhibits the more usual distribution b 2 2g a 2 1g e 2 g b 1 1g .

Fig. 2
Effects of oxygen vacancies on the phase competition and functionality of BiCoO 3 (BCO) thin films. Results were obtained for composition BiCoO 2.75 . 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 the lowest energies, either apical "Ap" or equatorial "Eq", are indicated within parentheses. "G" stands for G-type antiferromagnetism, "C" for Ctype 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 a in ; the green arrows represent the magnetic spin ordering in each phase. e Change in volume, ΔΩ, change in electric polarization orientation, Δα, and change in magnetic moment per formula unit, M i , associated with the multiferroic T -G → M-FiM phase transition.
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 Fig. 3c, e), which tends to restore (partially) the missing magnetic superexchange interactions along the out-of-plane direction. This superexchange restoration mechanism leading to the stabilization of AFM-G spin ordering, which is accompanied by a covalency increase in the out-of-plane Co1-O-Co3 bonds, is clearly observed in the electronic density plots shown in Fig. 3d-e (see yellow arrows therein). We have explicitly checked that the stabilization of the T -G(Ap) phase is robust with respect to the concentration of oxygen vacancies and presence of finite-size bias in our simulations (see "Methods" and Supplementary Discussion).
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 systems. For instance, first-principles works agree in predicting AFM-C spin ordering for T BFO thin films 48-50 whereas in most experimental studies AFM-G is reported for the same system 51,52 . As it will be shown later, by explicitly considering the presence of oxygen vacancies in BFO thin films such disagreements between theoretical and experimental works can be reconciled.
V O -induced stabilization of a ferrimagnetic M phase Figure 4 summarizes the electronic, structural, and magnetic properties of stoichiometric and non-stoichiometric 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 Å 53 lifts further the degeneracy in the e g and t 2g manifolds 54 , promoting the electronic occupation state b 2 2g e 2 g a 1 1g b 1 1g (Fig. 4a). Remarkably, when specific V O 's are created in equatorial positions the lowest-energy magnetic spin ordering becomes FiM. In particular, the net magnetization per formula unit associated with the M-FiM(Eq) phase is considerable: 0.5 and 0.25 μ B for the off-stoichiometries BiCoO 2.75 and BiCoO 2.875 , respectively. The two Co ions that are reduced by a neutral vacancy present the same magnetic moment orientation, equivalent electronic orbital occupancy, and sit within the [111] pc plane (indicated with an orange plane in Fig. 4b).
Interestingly, the 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).
For the pair of distant ions Co1 and Co4 to get reduced (Co 3+ → Co 2+ , Fig. 4c), the crystal needs to undergo sizable structural distortions involving the Bi and O atoms that are close to the equatorial V O (Supplementary Tables 6-7). For example, the nonmagnetic oxygen in the apical position above Co2 is displaced significantly out of Co2-Co4 vertical (Fig. 4c). Such a concerted electro-structural process has a clear imprint on the electronic density calculated for the plane oriented perpendicular to the substrate that contains Co1, Co2, Co3, and Co4 (Fig. 4d, e, see yellow arrows therein). As can be appreciated in Fig. 4d, e, the covalency of the Co4-O bond in the M-FiM(Eq) phase is reduced significantly as compared to that in the stoichiometric system. The vertical tilt of the Co2-O bond, which brings the apical oxygen closer to the equatorial V O , is also evidenced in Fig. 4e. We have carried out several numerical tests to assess the influence of the dimensions of the simulation cell and concentration of oxygen vacancies on the findings just described. Our results indicate that the revealed V O -induced stabilization of the M-FiM(Eq) phase is robust and general (see "Methods" and Supplementary Discussion).
The discovery of FiM spin ordering in the non-stoichiometric M phase motivated us to search for similar magnetic states, even if metastable, in the T and O BCO thin film geometries. The presence of highly magnetized O atoms was acknowledged in both stoichiometric phases ( Supplementary Fig. 7), however, upon removal of those oxygens no net magnetization emerged (neglecting small spin canting effects). These results suggest that the structural distortions just described are key to the stabilization  of FiM spin ordering. For instance, in the T phase the nonmagnetic oxygens appear in equatorial positions that are fairly clamped to the substrate ( Supplementary Fig. 7), hence they hardly can move off-symmetry. Meanwhile, in the O phase the lack of polar order and high dielectric permittivity appears to impose high energy penalties to any symmetry-lowering distortion (Supplementary Table 7), thus frustrating FiM spin ordering.
Based on the previous outcomes, we hypothesize that the following three conditions may be necessary for the stabilization of V O -induced FiM spin ordering in magnetic AMO 3 oxide perovskites: (1) lack of inversion symmetry leading to polar order and structural deformation ease, (2) moderate axial ratio structures allowing for the reduction of second-nearest transition metal ions, and (3) the existence of highly magnetized oxygen ions.

DISCUSSION
The results obtained in BCO thin films raise the natural question: are there any other multiferroic materials in which similar V Oinduced magnetic phenomena may occur? To answer this question we investigated the special 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 non-stoichiometric BFO thin films. In order to be consistent with the notation employed heretofore, we label the usual rhombohedral-like monoclinic phase of BFO 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) 50,55 . 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 for the T phase. (We recall that AFM-G spin ordering is metastable in the corresponding stoichiometric system 48,50 ). The causes of the stabilization of the T -G(Eq) phase in BFO are very similar to those explained previously for T -G(Ap) in BCO thin films (Fig. 3). As mentioned earlier, these results may shed new light on the origins of some unresolved discrepancies between theory and experiments regarding the identification of antiferromagnetic ordering in T BFO thin films [48][49][50][51][52] .
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 for BiFeO 2.75 . 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 M-G(Ap) (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 8-9). 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 could be entropically stabilized over the M-G (Ap) under increasing temperature due to the fact that the zerotemperature energies of the two states are similar and the vibrational free energy of the former phase appears to be more favourable 9,31 . Our results, therefore, can be interpreted as theoretical evidence showing that the observation of "ferromagnetic" behaviour in BFO thin films 24 may be caused by the presence of oxygen vacancies (as it was suggested previously by other authors 25 ). Figure 5d shows spin-up spin-down surface charge densities calculated for non-stoichiometric M and T BFO thin films considering in-plane and out-of-plane orientations. The reason why we could only find one FiM solution when generating neutral V O 's in stoichiometric M phases is now clear: all oxygen atoms in the M-G phase are non-magnetic whereas all apical O in the M-C phase are highly magnetized (Fig. 5d). The stoichiometric T -C phase also displays highly magnetized oxygen atoms in apical positions however, as we have explained before, it seems difficult to observe FiM spin ordering in c/a >1 geometries. Furthermore, we repeated the same type of calculations and analysis for monoclinic-like BiMnO 3 and BiCrO 3 thin films, which fulfill conditions (1) and (2) above. We found that the appearance of FiM spin ordering is also correlated with the presence of highly magnetized oxygen atoms (Supplementary Fig. 8). In particular, a M-FiM(Ap) phase displaying a net magnetization of ≈0.2 μ B is found in BiCrO 2.75 . These last results confirm that the simple rules provided in this work for understanding the appearance of FiM phases in non-stoichiometric multiferroic thin films may be general and robust.
In summary, by using first-principles calculations we have disclosed a number of previously overlooked electro-structural 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 specific 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 multifunctional 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 multifunctional materials. The present work shows that oxygen vacancies should be considered as a design opportunity to create new funcionalities in multiferroic thin films, especially to what concerns magnetism.

Density functional theory calculations
First-principles 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 56,57 . 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,30,31 . We use the "projected augmented wave" method 58 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 ffiffi ffi 2 p ffiffi ffi 2 p supercell containing 20 atoms (that is, four formula units) 55 . The reported net magnetizations correspond to sums over the local magnetic moments, which are accurate to within 0.01 μ B per formula unit.
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 (that is, |a| = |b| and γ = 90°). Periodic boundary conditions are applied along the three lattice-vector directions, thus the influence of surface and interface effects are systematically disregarded in our simulations. This type of calculations are known as "strained-bulk" geometry relaxations and typically are considered to be a good approximation for thin films with thicknesses of at least few nanometers 48,55 (i.e., for which surface and interface effects can be safely neglected). In our simulations, we assume systems that are elastically coupled to a substrate hence we disregard the existence of possible stress relaxation mechanisms in thin films. The geometry relaxations are stopped when the forces 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 59 and considering central differences for the calculation of atomic forces derivatives. Full phonon spectra calculations have not been pursued in this work due to the huge computational expense associated with large and low-symmetry offstoichiometric unit cells.
In this study, most non-stoichiometric configurations are generated by removing one oxygen atom from an apical or equatorial position in a 20atoms simulation cell, thus rendering the chemical composition The results presented in the main text apply to neutral oxygen vacancies since we have found that the creation of neutral V O 's in energetically more favourable than that of charged oxygen vacancies ( Supplementary Fig. 9, Supplementary Table 10, and Supplementary Discussion). In particular, to estimate the ranking of V O formation energies as a function of charge, q, we have employed the following well-established formula 60 : 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 finite-size 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 60 . 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. 9, Supplementary Table 10, and Supplementary Discussion).
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. 4-6, Supplementary Tables 11-13, and Supplementary Discussion). Specifically, we repeated most calculations by considering the PBEsol functional 61 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 energy ranking of the competing phases change in some cases (Supplementary Figs. 4-6, Supplementary Tables 11-13, and Supplementary Discussion). For a detailed discussion on these technical aspects, see Supplementary Discussion. It is also worth noting that lattice thermal excitations have been neglected systematically throughout this study 35 . Nevertheless, we have performed a test to assess the influence of this approximation in our main findings ( Supplementary Fig. 11 and Supplementary Discussion) and found that the conclusions reported in this study appear to be robust and physically meaningful also for T ≠ 0 conditions.
Regarding the estimation of the electric polarization, P, we started by employing the Berry phase formalism 62 . 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 55 : 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 a non-polar reference state (Z Ã κ has been calculated for each defective system). 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 for the P values estimated for non-stoichiometric thin films that are reported in Fig. 2.

Heisenberg model Monte Carlo simulations
To simulate the effects of thermal excitations on the magnetic ordering of T , O, and M BCO thin films, both stoichiometric and non-stoichiometric, we construct several spin Heisenberg models of the formĤ ¼ 1 2 P ij J ij S i S j , in which the value of the involved exchange constants are obtained from zero-temperature DFT calculations (following analogous methods to those explained in detail in previous works of ours 9,30,31 ).
For the T phase (c > a), we consider spin couplings between first nearest in-plane and out-of-plane neighbours (J a and J c , respectively) and between second nearest neighbours (J ac ). The value of such magnetic exchange constants are obtained with the formulas 9,30,31 : where E represents energy and superscripts "FM", "G", "A", and "C" ferromagnetic, antiferromagnetic G-type, antiferromagnetic A-type, and antiferromagnetic C-type spin arrangements, respectively. For the O and M phases (c ≈ a), we consider only spin couplings between first nearest neighbours and obtain their value with the expression 9,30,31 : The value of some characteristic magnetic exchange constants estimated in the manner just described are reported in Supplementary Table 14.
We use the Heisenberg models described above 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 Tdependence 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 P i ðÀ1Þ nix þniy S iz and S G 1 N P 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,30,31 .