Local polarization in oxygen-deficient LaMnO$_3$ induced by charge localization in the Jahn-Teller distorted structure

The functional properties of transition metal perovskite oxides are known to result from a complex interplay of magnetism, polarization, strain, and stoichiometry. Here, we show that for materials with a cooperative Jahn-Teller distortion, such as LaMnO$_3$ (LMO), the orbital order can also couple to the defect chemistry and induce novel material properties. At low temperatures, LMO exhibits a strong Jahn-Teller distortion that splits the $e_g$ orbitals of the high-spin Mn$^{3+}$ ions and leads to alternating long, short, and intermediate Mn--O bonds. Our DFT+$U$ calculations show that, as a result of this orbital order, the charge localization in LMO upon oxygen vacancy formation differs from other manganites, like SrMnO$_3$, where the two extra electrons reduce the two Mn sites adjacent to the vacancy. In LMO, relaxations around the defect depend on which type of Mn--O bond is broken, affecting the $d$-orbital energies and leading to asymmetric and hence polar excess-electron localization with respect to the vacancy. Moreover, we show that the Mn--O bond lengths, orbital order and consequently the charge localization and polarity are tunable via strain.

The family of doped (La,Ca/Sr)MnO 3 perovskite manganites has attracted great interest due to its very rich phase diagram and its unusual functional properties (colossal magnetoresistance [1][2][3][4][5], photo-induced infrared absorption [6], and efficient hole conductivity [7]) with strong potential for new applications in the fields of electronics, spintronics, and energy conversion. The orbital order plays a significant role in determining these properties and is strongly coupled with the other electronic, structural, and spin degrees of freedom [8,9]. LaMnO 3 (LMO) as the end-member of this family also exhibits interesting properties such as a pressure-induced insulatorto-metal transition [10] and a dielectric anomaly [11], that also depend on the orbital order [12,13].
Below 750 K, LMO adopts a distorted orthorhombic (P bnm) perovskite structure (see Fig. 1a) with an Atype antiferromagnetic order (A-AFM, Néel temperature T N ≈ 140 K [14]), with Mn atoms coupled ferromagnetic in the ac plane and antiferromagnetic along the b axis [15][16][17]. The stabilization of this orbital-ordered insulating state is a consequence of the high-spin Mn 3+ (t 3 2g e 1 g ) ions inducing a strong cooperative Jahn-Teller (JT) distortion that splits the e g orbitals (d z 2 /d x 2 −y 2 orbitals are alternately occupied within the orthorhombic ac-plane) and lead to alternating long and short in-plane Mn-O bonds, along with intermediate bonds along the b axis (see Fig. 1b). There are two symmetry-distinct oxygen atoms in this structure: the in-plane O (O IP ) in the ac plane with one short and one long Mn-O bond, and the out-of-plane O (O OP ) along the b direction with two intermediate Mn-O bonds (see Fig. 1c).
Unlike most other perovskite oxides, bulk LMO can exhibit oxygen superstoichiometry (LaMnO 3+δ ), which is accommodated through cation vacancies rather than oxygen interstitials [18][19][20][21]. LaMnO 3+δ has an extremely complex magnetic phase diagram exhibiting paramag- netic, ferromagnetic (FM), canted AFM order as well as a spin-glass state as a function of temperature and O superstoichiometry [19]. Topotactic low temperature reac- tions [22,23] were shown to lead to oxygen substoiochiometric LaMnO 3−δ with oxygen vacancies (V O ), where for 0.00 ≤ δ ≤ 0.20, the A-AFM order persists below T N ≈ 140 K [23]. The FM behavior often observed in LMO thin films [24][25][26][27][28][29][30][31] was explained by Mn 3+ -Mn 4+ double exchange due to the cation deficiency generally observed in these films. It was suggested that, below a critical thickness of 6 unit cells, LMO thin films become AFM [32], but several authors also reported thicker films with bulk-like AFM behavior obtained either through growth [33,34] or high-temperature annealing [35] under conditions favoring the formation of V O . V O formation is also expected to be affected by strain via the chemical expansion [36] caused by the reducing defects, tensile strain generally favoring the formation of V O [37][38][39].
In this letter we investigate strained stoichiometric and oxygen-deficient LMO by density functional theory (DFT) calculations within the Quantum ESPRESSO package [43,44] (see supporting information section for details). At the PBEsol+U [45,46] level of theory with a self-consistently computed Hubbard U for Mn-3d states [39,47], we correctly describe the JT distorted structure, which is crucial for accurately predicted magnetic and electronic properties [13,[48][49][50]. As can be seen from Table I, not only are the computed lattice vectors in good agreement with experiments (relative error below 1%), but this approach also provides a satisfactory description of the JT distortion (quantified via the magnitude of the Q 3 mode) and the Mn-O distances reproduced to within 2% of experiment. We find the expected orbitally ordered insulating (see Fig. S1 in the SI) state with d z 2 /d x 2 −y 2 orbitals being alternately occupied within the orthorhom- in previous theoretical studies [13,48], but still larger than experiment [42] (see Table I).
Upon oxygen-vacancy formation, the JT distorted structure leads to excess charge localization different from the one observed for example in CaMnO 3 and SrMnO 3 [37,39], where the two excess electrons reduce Mn sites in nearest neighbor (NN) sites to the defect. We consider neutral oxygen vacancies (V •• O in the Kröger-Vink notation [51]), but for simplicity refer to them as V O in the following. Figure 2c schematically illustrates the changes in the Mn-O framework induced by the relaxations around an out-of-plane oxygen vacancy (V OOP ).
After breaking the two equivalent Mn-O i bonds along the b-axis, the structure relaxes by shortening the Mn-O l bonds along c of the two undercoordinated NN Mn atoms (Mn 1 and Mn 5 in Fig. 2c), significantly expanding the former Mn-O s bonds along that direction of the NNN Mn atoms (Mn 4 and Mn 8 in Fig. 2c). As a result of these relaxations, the energy of the d z 2 orbitals increases for Mn 1 and Mn 5 and decreases for Mn 4 and Mn 8 (see Fig. 2d), localizing the two excess electrons resulting from V O formation on these orbitals and reducing the NNN Mn 4 and Mn 8 atoms (see Fig. 2c).
The case of an in-plane defect (V OIP ), is more complex since both a Mn-O s and a Mn-O l bond are broken. As can be seen from  [23] for LaMnO 2.80 that has a V O concentration very close to the one in our calculations. We note that this behavior was not observed in previous DFT calculations of oxygen vacancies in the high-temperature non-JT distorted G-AFM phase of LMO [52,53], where the reduction of the two NN Mn atoms was reported, thus supporting that our findings are directly related to the JT distortion in the A-AFM phase.
Interestingly, when the reduction does not take place on the Mn sites adjacent to the vacancy, the localization of the excess electrons with respect to the defect can be asymmetric and hence polar. We quantify this polarization using nominal charges, using a charge of +2 on the reduced Mn identified on the base of the structural changes and oxidation states calculated according to Ref. [54]. A small polarization of 4.06 µC/cm 2 in the ac-plane is observed for V OOP , in line with the reduced Mn 4 and Mn 8 (Fig. 2b), being symmetrically arranged with respect to the V O along the b-axis. A larger polarization (with IP and OP components of 16.38 and Based on the preceding discussion, it is reasonable to expect that any external parameter affecting the Mn-O bond lengths and/or the orbital order will have a strong impact also on the charge localization and the resulting polarization. As one such external parameter, we next investigate the effect of isostatic and epitaxial strain on stoichiometric and oxygen-deficient LMO. As shown in Fig. 3a and b isostatic and epitaxial strain in stoichiometric LMO are primarily accommodated by changes in octahedral tilt angles and bond lengths. In particular, while Mn-O s and Mn-i are relatively insensitive to strain, the Mn-l bonds are more strongly affected especially under tensile strain. This results in an increase of the JT mode amplitude Q 3 under tensile strain (Fig. 3c), which in turn induces an increase of the LMO band gap (Fig. 3d). The effect of compressive strain is opposite to that of tensile strain, reducing the Q 3 amplitude and hence the band gap.
As shown in Figure 4a, when isostatic strain is applied in presence of a V OOP , we observe changes in excesscharge localization as a function of the applied strain. For compressive strain larger than -2%, one of the excess electrons localizes on a Mn further from the defect compared to the NNN sites reduced without strain. This is a consequence of the reduced Mn-O l bond lengths (see Fig. 3). For tensile strain larger than 2%, instead, the NN Mn sites are reduced. A similar behavior is observed for V OIP (Figure 4b), but in this case, even at +4% strain, only one of the NN Mn atoms is reduced (the one previously forming a Mn-O s bond). Despite the reduction of NNN Mn sites in some strain ranges, the formation energies of both the V OOP and V OIP are reduced under tensile and increased under compressive strain, with discontinuities when the excess-charge localization changes. This general trend is expected from chemical expansion arguments [37]. Epitaxial strain, induced for example by heteroepitaxial growth on a substrate with a different lattice parameter, also affects the excess-charge localization, but the results shown in Figures 4c and d are qualitatively different compared to isostatic strain. For a V OIP , the charge localization changes only for tensile strain larger than 3% (see Fig. 4d). Even at these large tensile strains, the reduction never takes place on NN but always NNN Mn sites, which is due to the non-isotropic and weaker effect on the bond lengths of epitaxial compared to isostatic strain [55]. Consequently we also do not observe the chemical expansion behavior but instead the opposite trend with a reduced formation energy under compressive strain that was previously reported for ionized oxygen vacancies [55]. Moreover, at variance to isostatic strain, the IP component of the polarization is reduced for +3% and larger strain due to the change of the relative Mn' orientation in the ac-plane, while the OP component stays fairly constant.
Interestingly, for V OOP , simply imposing the constrain of a cubic substrate with an a lattice parameter resulting in the same in-plane area as LMO (0% strain, Fig. 4c) induces a more asymmetric charge localization compared to the bulk, which is maintained also under compressive epitaxial strain. Between +1 and +2% strain, the localization of the two excess electron changes, the reduction still taking places on NNN Mn atoms and resulting in an increase of the formation energy, similarly to V OIP . However, differently from the latter, for larger tensile strains one (+3% strain) or two (+4% strain) excess electrons localize on NN Mn sites and the formation energy remains almost constant.
In summary, we have shown that, in the case of Jahn-Teller distorted materials, orbital order is an additional parameter to take into account during defect-based design of functional material properties. As we have shown for the example of oxygen vacancies in LaMnO 3 , unexpected excess-charge localization results from the orbital order and can lead to the emergence of local polarization. The charge localization and hence the polarization magnitude are tunable by strain, suggesting the possibility to engineer ferroelectricity in LaMnO 3 when local dipoles couple for high enough defect concentrations.

ACKNOWLEDGMENTS
This research was supported by the NCCR MAR-VEL, funded by the Swiss National Science Foundation. Computational resources were provided by the University of Bern (on the HPC cluster UBELIX, http://www.id.unibe.ch/hpc), by the Swiss National Supercomputing Center (CSCS) under projects ID mr26 and SuperMUC at GCS@LRZ, Germany, for which we acknowledge PRACE for awarding us access.
The calculations were performed in the framework of DFT with the Quantum ESPRESSO package [1,2]. PBEsol [3] was used as exchange-correlation functional together with ultrasoft pseudopotentials [4] including La(5s,5p,5d,6s,6p), Mn(3p,4s,3d), and O(2s,2p) states. A kinetic-energy cutoff of 90 Ry for wave functions and of 1080 Ry for spin-charge density and potentials were applied. A gaussian smearing with a broadening of 0.01 Ry was used in all calculations.
LMO (space group P bnm) was studied using a 40-atom supercell of the 5-atom primitive cubic cell. A 4×4×4 Monkhorst-Pack k-point grid was applied to sample the Brillouin zone of this cell. A denser 8×8×8 grid was used for plotting the density of states (DOS). Only the A-type antiferromagnetic (AFM) phase with ferromagnetic coupling in the ac plane was considered. Neutral oxygen vacancies (V •• O ) were created by removing one O atom from the supercell. We studied two inequivalent oxygen-vacancy positions: an out-of-plane (OP) configuration and an in-plane one (IP) in which the oxygen removal breaks Mn-O-Mn bonds along the b-axis and in the ac-plane, respectively. For stoichiometric bulk and isostatic calculations, both lattice parameters and atomic positions were relaxed. Isostatic strain was applied by expanding or shrinking the cell vectors by the same amount in all directions. Thin film geometries with epitaxial biaxial strain in the ac-plane to simulate growth on a cubic substrate were instead computed following the procedure reported in Ref. [5]. Finally, for defective systems, atomic positions were optimized while keeping the lattice vectors fixed at optimized values of the stoichiometric system. In all cases, atomic forces were converged to within 5 × 10 −2 eV/Å, while energies were converged to within 1.4 ×10 −5 eV.
The rotationally invariant formulation by Dudarev et al. [6] was used in all the DFT+U calculations. We imposed a global Hubbard U parameter (U SC ) on all the Mn-3d states. The U SC value of 4.08 eV was computed via DFPT calculations [7] and through the self-consistent procedure introduced in Ref. [8] for stoichiometric bulk geometries of the A-AFM phase. Γ point sampling of the q-space was performed for all DFPT calculations. A convergence threshold of 0.01 eV was applied for the self-consistency of the U values. Atomic orbitals were used to construct occupation matrices and projectors in the DFT+U scheme.
The defect formation energy for a oxygen vacancy in its neutral state (E f,VO ) was computed with the following equation (see Ref. [9]): where E tot,VO and E tot,stoic are the total energies of the defective and stoichiometric systems, is the applied strain and µ O is the oxygen chemical potential. We will show results in the oxygen-rich limit, i.e. with µ O = 1 2 E O2 , E O2 being the energy of an oxygen molecule.
To reduce the computational cost, the polarization P was computed using a point-charge model: where r i is the position of atom i and q i is its nominal charge: +3 for La, -2 for O, and +3 or +2 for stoichiometric-like or reduced Mn sites. The charge applied on each Mn was defined on the base of its oxidation state computed through the method introduced by Sit et al. [10]. The polarization being a multivalued quantity, it has been corrected by an integer number of polarization quanta Q, computed as: with a, b, and c being the lattice parameters, V the volume of the unit cell, and e the elementary charge.
ELECTRONIC PROPERTIES Fig. S1 shows the computed electronic density of states, showing the insulating character of JT-distorted stoichiometric LMO, the Fermi energy lying between the occupied and unoccupied e g states.
FIG. S1. Total and projected density of states (DOS) for bulk LMO in the A-AFM phase. The zero of the energy scale was set at the Fermi energy value