Applying a new interatomic potential for the modelling of hexagonal and orthorhombic YMnO3

We develop and apply an interatomic potential for YMnO3, based on the shell model together with the angular overlap model, which can model ligand field effects. The potential parameters accurately reproduce the complex structure of both hexagonal and orthorhombic phases of YMnO3. The rotation of the MnO6 octahedra in o-YMnO3 suggests the E-type AFM order. The potential is further employed to investigate the energies of intrinsic defects in the material. Lower defect energies were found in o-YMnO3. Oxygen Frenkel and Y2O3 partial Schottky are the most favourable defects in h-YMnO3 and o-YMnO3, respectively. The defect models proposed have implications for the properties of the related non-stoichiometric phases.


Introduction
The rare-earth manganite, YMnO 3 , is a well-known example of a multiferroic which possesses magnetism and ferroelectricity simultaneously and thus has attracted considerable attention in recent years due to its multiferroic and magnetoelectric properties and potential interest. 1-3 YMnO 3 prepared under ordinary synthetic conditions crystallizes in a hexagonal structure, which belongs to the space group P6 3 cm, and can be considered as alternating ab-layers of Y 3+ ions and corner-sharing MnO 5 trigonal bipyramids, as shown in Fig. 1(a). Each MnO 5 bipyramid is formed of a central Mn 3+ ion surrounded by three planar oxygen (O3 and O4) atoms and two apical oxygen (O1 and O2) atoms. Whilst there is one unique site for manganese atoms and four for oxygen atoms, there are two inequivalent sites for yttrium atoms. 4 By means of softchemistry synthesis, applying pressure, or epitaxial strain in thin films, the hexagonal structure can be converted into the more dense, albeit metastable orthorhombic structure (Pnma space group), [5][6][7][8][9][10][11][12][13] which contains two distinct oxygen sites (four planar O1 sites and two apical O2 sites), while there is only one inequivalent site for Y atoms and likewise for Mn atoms. 7 The orthorhombic structure is composed of corner-sharing MnO 6 octahedra with one yttrium cation occupying each hole; see Fig. 1(b). Compared to the cubic phase, the unit cell of which is composed of just one regular octahedron, the MnO 6 octahedra are both rotated (Mn-O bonds originally align parallel with the crystallographic axes) and distorted (the six equivalent Mn-O bond distances become three degenerate pairs) as a result of steric and Jahn-Teller effects. 25 Both the hexagonal and orthorhombic phases of YMnO 3 have been subjected to theoretical studies based on Density Functional Theory. 2,6,[14][15][16][17] Many questions, however, remain unsolved that lie beyond the computational limits of applying purely quantum mechanical approaches, especially problems related to complex defect structures, and other configurations that require the inclusion and relaxation of a large number of ions. Atomistic modelling techniques -employing interatomic potentials based on the ionic model which has already been successfully employed to a range of oxide-based materials 18,19 -have been applied to the study of YMnO 3 20,21 but, due to the purely radial feature of the potentials, the calculations showed major limitations in accurately modelling its structure which has a substantial distortion of the MnO 5 polyhedron in the hexagonal phase and the MnO 6 polyhedron in the orthorhombic phase.
In the present paper, we employed a semi-empirical force field approach, which includes the angular overlap model (AOM). 22 The AOM model has been successfully applied to modelling the Jahn-Teller distortions within and the rotations of the MnO 6 octahedra in the orthorhombic perovskite LaMnO 3 , 23 oxygen anion migration in LaMnO 3Àd , 24 and the changes in the structure that occur across the series of lanthanide manganates. 25 Here, we develop a set of potential parameters that accurately reproduce the structure of YMnO 3 in both the hexagonal phase and metastable orthorhombic phase and, moreover, successfully model the thermodynamic parameters relating to the phase transition. For h-YMnO 3 , the lattice distortions are associated with its ferroelectric properties 3,4 and the strong interactions between Mn trimers, which dominate the magnetic and magnetodielectric coupling. 26 For o-YMnO 3 , the magnetic ordering and, consequently, the multiferroic and magnetoelectric properties of this phase depend on the bond angle. 9,27-31 Thus, the capability of the present method to predict accurately the bond length, bond angle, atomic displacement and distortion features of YMnO 3 under hydrostatic pressure and chemical doping is of great interest in the multifunctional properties of RMnO 3 based materials, where R is a rare earth element. Additionally, considering that the large leakage currents with a small quantity of defects induced by the low band gap of YMnO 3 hamper their application, [32][33][34][35] it is of great importance to know the origin of the defects in the material, especially for h-YMnO 3 , which is ferroelectric at room temperature. Consequently, we also perform a detailed investigation of the intrinsic defect properties which also relate to the properties of the related non-stoichiometric phases.

Simulation method
Our calculations are based on the Born model, using interatomic potentials that have been modified by the addition of a ligand field term. The contribution from ''spherical'' forces between ions to the lattice energy is described by the standard Coulomb-Buckingham expression: where the summation over i and j includes unique pairs of atoms and parameters q, A, r and C are species dependent parameters.
Three-dimensional periodic boundary conditions were applied to the unit cell, i restricted to atoms within this unit cell and contributions to the energy reduced by a half when atom j is a periodic image so that eqn (1) gives energy per unit cell. The Ewald summation 36,37 was employed to compute the Coulomb term, using formal charges, q, on the ions. The remaining terms constitute the short-range Buckingham potential representing the cation-anion and anion-anion short-range interactions; a cut off of 12 Å was applied. The shell model 38,39 is employed to describe the polarizability of individual ions. In this model, each ion is divided into two coupled parts: a core with mass m and charge X and a massless shell with charge Y (we constrained X + Y to be the formal charge of ion). The interaction between the core and corresponding shell is described by an harmonic function with a spring constant k. In our calculations, A, r, C, k and Y were fitted empirically so that the structural parameters of the appropriate phases for YMnO 3 were reproduced.
To model the irregular coordination geometry associated with Mn 3+ ions, we add a contribution to the lattice energy that corresponds to the Jahn-Teller driving force for the distortions and removal of degeneracy in electronic energy: 40 Here, e t d and O ts d are the energy changes and occupations of the d-orbitals for each transition metal, t, respectively, and s is the spin of an electron. An adaptation of the angular overlap model (AOM) is used to obtain e t d . In this model we compute the eigenvalues, e d , of a 5 Â 5 overlap matrix, H dd 0 for each transition metal ion. H dd 0 is formed by taking the products of the angular contributions to the overlap integrals, G d , 41 between the d-orbital of transition metal ion and the orbitals of any surrounding ligand (l), where the Born-Mayer interaction is used to model the radial dependence of the interaction between transition metal ion and its ligand: The two new parameters, A LF and r LF , depend on ligand type and can be empirically fitted. Upon optimization of the lattice energy with respect to the cell parameters and ionic coordinates, the energy levels may become degenerate and the order of the energy levels may also change. To prevent the energy landscape becoming discontinuous, we employ partial occupancies (or a nonzero probability that an electron can occupy a higher energy state) via the implementation of a Fermi function. The theory of the AOM and its successful application in modelling compounds containing ''non-spherical'' transition metal ions, i.e. manganites with distortions, are reviewed in detail elsewhere. 40 The lattice-energy is minimized by relaxing both cell dimensions and atomic coordinates at constant pressure using a quasi Newton-Raphson procedure together with the BFGS method 42 for updating the Hessian. All our calculations employed the General Utility Lattice Program (GULP). 43,44 The potential parameters derived in our fitting procedure are reported in Table 1.

Results and discussion
Lattice calculations The first challenge, for our newly derived potential model, is to reproduce the crystal structures of both the hexagonal and orthorhombic structures with the same set of potential parameters. The energy minimised structural parameters and experimental data are listed in Table 2. Results calculated from previously published force field parameters are also presented for comparison. We note that the differences between experimental 5,45 and our calculated lattice parameters and bond lengths are all within 2%, and in most cases less than 1%. Hence, our potential parameters reproduced the complex crystal structure of YMnO 3 for both the hexagonal and orthorhombic phases, which is in stark contrast to previous results. 20,21 For example, the largest bond length error was reduced from 4.6% and 6.98% to 1.56% for the hexagonal phase, and from 9% and 7.98% to 2.38% for the orthorhombic phase. It is noteworthy that the bond angle in the orthorhombic phase is now within 2.04% (cf. 3.5% and 4.42% for ref. 20 and 21, respectively). According to the magnetic phase diagram for orthorhombic RMnO 3 as a function of Mn-O-Mn bond angle, [27][28][29][30][31] our results show the E-type AFM order of o-YMnO 3 , while other results 20 determine the magnetic phase to be A-type. Moreover, comparing the predicted structural data obtained from potentials with and without the AOM ligand field term, we note that the ligand term effectively reproduces the asymmetry of the Mn 3+ ions, i.e. the distortion and rotation of both MnO 5 bipyramid in h-YMnO 3 and MnO 6 Jahn-Teller octahedron in o-YMnO 3 .
The robustness of the proposed potential set is further checked through investigation of the pressure dependence of both the orthorhombic and hexagonal structures (Fig. 2). Our calculations showed that the orthorhombic phase became more stable at the pressure of 25.8 GPa, which accords well with the experimental observation that, at room temperature, a pressure-induced hexagonal-orthorhombic phase transition requires a pressure above B22 GPa for YMnO 3 . 46 This result is gratifying in view of the considerable sensitivity to details of the potential model of the thermodynamic parameters associated with this kind of phase transition.

Frenkel and Schottky disorder
We next present results of our defect calculations for YMnO 3 . We use the supercell method rather than the Mott-Littleton approach in our study of defects within both the hexagonal and orthorhombic structure as the use of potential functions Table 1 Buckingham parameters for the interaction between the shells of the ions; the force field parameters for radial part of the AOM that acts between each manganese core and the surrounding oxygen cores; and the parameters for the shell model. Subscripts s and c indicate shell and core, respectively. The cut-off radius for all short range potentials is 12 Å  needed to describe Mn 3+ -O 2À interactions in the non-cubic phase is not currently implemented in the latter. Supercells were generated from the relaxed cell parameters and ionic coordinates of bulk YMnO 3 . After removing or adding one ion from or into the supercell, we added a neutralizing uniform charge background. These structures were then optimized with the cell parameters being fixed, since one isolated defect should not change the bulk lattice constants. To obtain the energy required to form an isolated defect, E f , we first calculate the energy difference, DE, between the lowest local energy minimum obtained from the optimization of the defect containing cell and the lattice energy of the perfect crystal. We require this energy to be converged with respect to the size of the supercell, or shortest distance between the defect and the unwanted images of this defect. In order to speed up this convergence, we reduce the unwanted long-range interactions between the defect (vacancy and local distortions) and its images, by adding a correction term E c to DE. 24,47 where a is the Madelung constant and e r the dielectric constant of the perfect crystal. This term refers to the Coulomb energy of a point charge Q (charged defect) immersed in a structureless dielectric, within a cubic unit cell of length L with a neutralizing uniform charge background. 47 With increasing size of the supercell, the contribution given by E c will dominate the energy difference between the lattice with a periodic array of defects and that with an isolated defect, i.e. the energy difference caused by the unwanted interaction mentioned above.
In order to estimate E c , we first compute For a single point charge Q, E D is very straight forward to obtain from GULP as it is the lattice energy of the cation or anion, which was removed in order to create the defect, in an empty unit cell that has the same dimensions as the supercell it was removed from. Note that a charged uniform neutralizing background is also required. This procedure can be generalised to include more point charges, and will be discussed elsewhere. The interactions of the charge associated with the defect and that of its images will be reduced by the screening effect of the remaining ions, and so the correction term, E c , is estimated as: where e r is the average of the diagonal components of the diagonalized static dielectric constant tensor. The defect formation energy, E f , is therefore Fig. 3 shows the results of isolated oxygen vacancy formation energy of hexagonal and orthorhombic YMnO 3 . We investigated all unique oxygen sites. By adding E c , the term in eqn (7), the convergence of E f was improved considerably in both phases. In our final results, we extrapolate the value for E f at 1/R = 0. This supercell method has also been compared with the Mott-Littleton approach before and proved to be reliable and effective. 24 The energies of isolated vacancy and interstitial defects were calculated first for both the hexagonal and orthorhombic structures. Regarding the vacancies, we calculated both the yttrium sites and the four oxygen sites in the hexagonal structure and the two oxygen sites in the orthorhombic structure. For h-YMnO 3 , the planar oxygen O3 and O4 sites have a lower energy -about 1 eV lower than the apical sites, which is consistent with experimental observation. 48 It is reported that the host h-YMnO 3 shows a preference for removing the equatorial oxide anions (at O3 and O4) to yield phases of composition YMnO 3Àd , which may be rationalized by considering the coordination requirement of cations. 49 The formation energy of the Y2 vacancy is slightly lower than Y1 by 0.25 eV. For the orthorhombic phase, the difference in E f between different sites is much smaller. The O2 vacancy is a little more energetically favourable by 0.65 eV. Only the lowest energies were used in the following calculations.
As to interstitial defects, for the hexagonal structure, we tested several possible positions to confirm the optimal position of the interstitial site; for the orthorhombic phase, we considered two possible sites: the octahedral interstice at (0, 1/2, 0) and the tetrahedral interstice at (1/4, 1/4, 1/4), and lower energies were found for interstitials placed at the octahedral site. Only the lowest formation energies were employed in the calculations below.
Then, Frenkel disorder, full Schottky disorder and partial Schottky disorder energies were calculated by combining individual defect energies and lattice energies. These defect reactions are described by the following equations, where Kröger-Vink notation is used.
Yttrium Frenkel disorder: Manganese Frenkel disorder: Oxygen Frenkel disorder: YMnO 3 full Schottky disorder: Y 2 O 3 partial Schottky disorder: Mn 2 O 3 partial Schottky disorder: The corresponding energies for all these types of intrinsic defects are listed in Table 3. In order to enable comparisons between different types of disorder reactions, the energies are given as defect formation energies per defect, i.e. E F /2 for Frenkel disorder and E S /5 for Schottky disorder. For both systems, the cation Frenkel disorder energies are much higher than other disorder energies. However, in the hexagonal structure, the oxygen Frenkel energy is the lowest defect type; while in orthorhombic structure, it is slightly higher than but comparable to Schottky disorder energies, which indicates that vacancies, not interstitials, will be the dominant structural defects. For both structures, however, defect energies are so high that intrinsic disorder would not be expected to dominate the defect chemistry, although we find that very significantly lower energies were obtained for the orthorhombic structure.
Essentially, from our results obtained using the supercell method, it is predicted that the oxygen Frenkel disorder and Schottky disorder involving Y are found to be the most energetically favourable intrinsic defect in h-MnO 3 and o-YMnO 3 , respectively, although the defect energies are still high. The greatest significance of our calculations is perhaps not in the low levels expected for intrinsic disorder, but their predictions for related non-stoichiometric phases. The variable valence of Mn should result in both oxygen deficient and oxygen excess phases depending on the oxygen partial pressure. Our calculations predict that for oxygen deficient phases the reduction of Mn should, as expected, lead to oxygen vacancy formation. The results for oxidation are more interesting as there are more possibilities. And we will discuss this in more details below.

Oxidative nonstoichiometry
As mentioned earlier, oxidation in the hexagonal phase and orthorhombic phase might happen via different mechanisms, consequently we have examined the energetics of four different oxidation reactions (eqn (15)-(18)): The hole state in both phases was modelled as a Mn 4+ ion (i.e. as a small polaron), as observed in oxygen excess systems 5 and hole doped systems. 50 In addition to the individual defect energies and the lattice energies, we also included the dissociation energy of an oxygen molecule (5.16 eV), the first and second electron affinities of oxygen (À1.46 and +8.75 eV respectively) and the fourth ionisation energy of manganese (52 eV). 51 This is clearly an approximate approach but should still provide useful comparisons between the different processes listed above. The calculated reaction energies are listed in Table 4. For the hexagonal phase, the results show that the energy for reaction (15) is the lowest, which suggests that oxidation will result in the formation of oxygen interstitials. This prediction is in agreement with experiment studies in which hexagonal YMnO 3+d was observed to accommodate interstitial oxygen. 52,53 Moreover, h-Dy 1Àx Y x MnO 3+d has been investigated as a possible material for storing/releasing oxygen; as high as 2000 mmol-O per g in air is observed when cycling between stoichiometric and interstitial oxygen rich structures by changing temperature between 250 and 350 1C. 54 For the orthorhombic system, we predict yttrium vacancy formation during oxidative reactions. The difficulty of preparation of YMnO 3 in its perovskite form results in the relative scarcity of reports devoted to the description of its basic properties. Insight perhaps can be gained from studies of its larger A-site cation isomorphous LaMnO 3 counterpart, which have reported that LaMnO 3 has the ability to accommodate oxidative nonstoichiometry via cation vacancies. 55,56 There is some disagreement, however, in the literature concerning whether A-site or B-site cation vacancies dominate in oxygen excess compositions. Mitchell et al. 56 used NPD methods and found that there was a tendency towards more A-site lanthanum vacancies present in orthorhombic LaMnO 3 structures, which accords with our calculations. Importantly, the lower reaction energies for the orthorhombic phase would lead us to expect greater non-stoichiometry than for the hexagonal phases. It would be interesting to see whether this prediction can be verified experimentally.
To summarise, the oxidation reaction is found to occur via oxygen interstitial and yttrium vacancy formation for h-MnO 3 and o-YMnO 3 , respectively, and their ability to accommodate oxidative nonstoichiometry is to be expected.

Conclusions
In conclusion, we have developed a new set of interatomic potential parameters based on the AOM model, 23 which successfully reproduces the structure of both hexagonal and orthorhombic structure of YMnO 3 . Furthermore, we have used this potential to examine basic defect energies. In the hexagonal structure, the oxygen Frenkel defect is the most energetically favourable, while Schottky defects involving Y ion have the lowest energies in the orthorhombic phase. Significant intrinsic disorder of either Frenkel or Schottky type is, however, unlikely, but the results have considerable significance for the nature of related non-stoichiometric phases. The capability of the proposed method to estimate accurately the structural properties of both h-YMnO 3 and o-YMnO 3 , and the detailed investigation of the defects in the present work have implications for modified RMnO 3 based multifunctional materials.