Structural and ferroelectric transitions in magnetic nickelate PbNiO3

Density functional calculations have been tremendously useful in understanding the microscopic origin of multiferroicity and in quantifying relevant properties in many multiferroics and magnetoelectrics. Here, we focus on a relatively new and promising compound, PbNiO3. The structural, electronic and magnetic properties of its two polymorphs, i.e. the orthorhombic structure with space group Pnma and the rhombohedral LiNbO3-type structure with space group R3c have been studied by using density functional calculations within DFT + U and hybrid functional schemes. Our data convey an accurate description of the pressure-induced phase transition from the rhombohedral to orthorhombic phase at a predicted critical pressure of 5 GPa in agreement with the measured value of 3 GPa. Both phases show the G-type antiferromagnetic configuration as a magnetic ground state, but differ in the spatial anisotropy associated with nearest-neighbor exchange couplings, which is strongly weakened in the rhombohedral LiNbO3-type phase. The predicted large ferroelectric polarization of the rhombohedral phase (Hao et al 2012 Phys. Rev. B 014116) has been re-explored and analyzed in detail using partial density of states, Born effective charge tensors, charge density difference, electron localization function analysis and distortion mode analysis. The asymmetric bonding between the Pb 6s and O 2p orbitals along the [111]-direction is responsible for the polar cationic displacement, giving rise to a predicted large ferroelectric polarization as high as ∼ 100 μC cm−2.

charge tensors, charge density difference, electron localization function analysis and distortion mode analysis. The asymmetric bonding between the Pb 6s and O 2p orbitals along the [111]-direction is responsible for the polar cationic displacement, giving rise to a predicted large ferroelectric polarization as high as ∼100 µC cm −2 .

Introduction: a revised panorama on multiferroic compounds
Magnetic and ferroelectric materials are ubiquitous in modern science and present-day technology for their diverse applications, ranging from data storage to magnetic and ferroelectric random-access memories and spintronic devices [1][2][3]. Both ferromagnets and ferroelectrics are characterized by a spontaneous symmetry breaking that causes the appearance of a switchable long-range magnetic or dipolar ordering below a critical temperature, and they are commonly referred to as ferroics [4]. In recent years, the quest for (magnetoelectric) multifunctional integration within a single material has motivated a renewed interest in the class of so-called multiferroic materials, displaying the simultaneous presence of two or more spontaneous ferroic phases. The intrinsic combination of magnetism and ferroelectricity in this class of compounds could in principle allow for novel device paradigms exploiting their cross-coupled effect, e.g. allowing to control magnetization (polarization) by an external electric (magnetic) field. The intense research activity on this class of materials, both from the experimental and theoretical sides, is testified by the number of excellent review papers on the topic, opening promising prospects for applications also beyond the field of magnetic and ferroelectric materials, e.g. for information and energy harvesting technologies [5][6][7][8][9][10].
From the theoretical point of view, the intrinsically coupled and multifunctional nature of multiferroics poses numerous fundamental challenges. In fact, if on the one hand the theoretical understanding of magnetic insulators is rather well established, on the other hand a rigorous microscopic theory for ferroelectrics has been formulated only in the last few decades [11][12][13]. Naïvely, the difficulty arises from the fact that magnetic properties can be qualitatively traced back to well-defined localized quantities, i.e. the spin moments; conversely, electric dipoles in extended-state systems are well defined only in the very special case of a fully ionic material, while in general the electric polarization is a global property of the matter that cannot be decoupled in localized contributions, and is subject to change as a function of temperature. As a matter of fact, a rigorous definition of ferroelectric bulk polarization has been provided only in the early 90s via a quantum-mechanical approach based on Berry phases [11][12][13]. Furthermore, the microscopic origin of ferroelectricity is not unique, and the coexistence and interplay between structural distorsions and electronic degrees of freedom (spin, charge, and orbital) has been proven to play a central role in the diverse mechanisms devised so far. In this respect, first-principles calculations based on the density-functional theory (DFT) have played an important role in the description, understanding as well as prediction-via identification of suitable material design rules-of magnetic, ferroelectric and magnetoelectric properties of mutiferroics, due to its ability to describe the many active degrees of freedom within a comparable level of accuracy [3,[14][15][16].
Despite the advances reached in recent years, two major accomplishments still need to be fully achieved, namely (i) how to obtain an efficient coupling between ferroelectric and magnetic orders in the multiferroic phase and (ii) how to realize a room-temperature multiferroic. The first issue is also directly related to a more fundamental question, which indeed lied at the origin of the renewed interest in multiferroics and still lacks a completely satisfactory answer: which conditions and constraints can determine the possibility to combine magnetism and ferroelectricity in a single homogenous material? For a long time it was thought that the two phenomena were mutually exclusive, and the formulation of an empirical exclusion rule claiming a chemical incompatibility between magnetism and ferroelectricity essentially slowed down the study of multiferroics by the 1970s. This exclusion rule was motivated by an early observation: most of the long-known (conventional) ferroelectrics are perovskite transition-metal oxides (BaTiO 3 , PbTiO 3 , KNbO 3 , and LiNbO 3 ), in which the polar distortion mainly involves an off-centering of the perovskite B-site transition-metal cation showing an empty d shell [17]. As rationalized later, the driving force behind the polar distortion in these conventional ferroelectrics comes from the tendency of empty d states to establish some degree of covalency with the surrounding oxygens [18][19][20]. On the other hand, magnetism requires partially filled d states to appear in perovskite oxides, thus apparently explaining the small overlap between the large families of ferroelectric and magnetic materials with the perovskite structure [21,22]. A major breakthrough in the early 2000s was realizing that (i) the covalent bonding-mediated mechanism may not be the only possible origin of ferroelectricity and (ii) in some cases an electric polarization (called improper) may appear even if the polar displacement mode is not the only or the primary order parameter of the paraelectric-ferroelectric transition. A case in point is when a specific spin ordering breaks the inversion symmetry, thus causing the appearance of an improper ferroelectric phase with an intrinsically large magnetoelectric coupling. The latter phenomenon is observed in rareearth manganites (TbMnO 3 and HoMnO 3 ) [8,9]; unfortunately, these systems are usually frustrated (anti)ferromagnets, typically with rather low critical temperatures (T c 40 K). Furthermore, their frustrated magnetic configurations result in complex magnetoelectric effects (such as an electrical-induced swap of the spin-rotation planes in TbMnO 3 ) which may not be straightforwardly exploitable in applications.
It is well worth noting, however, that very recent theoretical advances pointed out that the aforementioned exclusion rule is less strict than expected even in the class of proper ferroelectric oxides. In fact, the conditions of multiferroicity in d n perovskites have been recently derived from the pseudo Jahn-Teller effect, which considers ferroelectric displacements as triggered by vibronic coupling between ground and excited electronic states of opposite parity but same spin multiplicity [23]. According to this somewhat simplified approach, only some specific d n and spin configurations allow for a coexistence of magnetism and ferroelectricity. A particularly interesting case is that of high-spin d 3 configuration, a situation realized in alkaline-earth manganites such as CaMnO 3 or SrMnO 3 ; these systems were shown to possess a weak ferroelectric instability mediated by a covalency-driven mechanism (such as in BaTiO 3 ), which at ambient conditions is hindered by other energetically favorable distortions, namely non-polar rotation/tiltings of BO 6 octahedra [24][25][26]. Applying strain or chemical pressure could in principle tune the balance between different structural instabilities, allowing for ferroelectricity to emerge, as indeed recently shown for strained CaMnO 3 films [27] and bulk Sr 1−x Ba x MnO 3 [28]. Since both magnetic and ferroelectric instabilities are related to the same manganese B cations, the magnetoelectric coupling is expected to be strong in these systems [28][29][30], even though antiferromagnetic phases are apparently still favored; on the other hand, critical temperatures may be far larger than in frustrated magnets (e.g., T c ∼ 185 K in Sr 0.5 Ba 0.5 MnO 3 [28]). Another prime example of the coupling between magnetism and ferroelectricity is represented by B-site magnetic ferroelectric perovskites such as MnMO 3 (M = Ti, Sn) in which the ferroelectric instabilities are driven by the chemical activity of the B-site atom [31].
Interesting theoretical advances were recently achieved noting that usually competing distortions in transition-metal perovskites, such as tilting/rotations of BO 6 octahedra, can in principle be tuned in order to trigger hybrid improper ferroelectricity [32,33]. Such rotations in simple AB X 3 crystal structures generally preserve inversion symmetry, causing an 'antiferroelectric' ordering of A-site ions; these instabilities usually compete with covalencydriven B-site off-centerings, as it happens in e.g. CaTiO 3 or SrTiO 3 . However, several strategies have been proposed from a material design approach, which would allow the realization of ferroelectric systems arising from centric rotations, e.g. by layering perovskite blocks or by partial substitution of A-site ions [32,[34][35][36][37][38][39][40][41][42]. Interestingly, the magnetic properties may be also tuned and magnetoelectric effects are likely to appear, since rotations of corner-sharing octahedra buckle the inter-octahedral B-O-B angles which mediate the interplay of electronic and magnetic degrees of freedom. As an example, a theoretical calculation has shown that a linear magnetoelectric effect, allowing for electric-field control of magnetization, is predicted to appear in a layered Ca 3 Mn 2 O 7 [37]. Due to the generally high temperatures at which these non-polar distortions are expected to appear, as compared to frustrated magnetic phases, in the near future the proposed approach could lead to the identification of new multiferroic materials at room temperature. Another possibility is to explicitly dope with magnetic dopants (e.g. Ti-V substitution) the layered La 3 Ti 2 O 7 structure [39,40] which in itself displays robust ferroelectricity due to the breaking of the octahedral rotation pattern in the direction parallel to the stacking [38]. Furthermore, other typical distortions in perovskites, such as non-polar Jahn-Teller ones, can in principle play a role in the appearance of both hybrid improper ferroelectricity and magnetoelectric coupling, as recently shown in a new class of materials, namely metal-organic frameworks with perovskite geometry [43,44]. As opposed to their inorganic counterparts, metal-organic frameworks offer the possibility of varying their additional degrees of freedom, arising from organic/inorganic duality, which could be exploited for a rational design of new materials with enhanced functionalities.
With respect to the issue of having a room-temperature multiferroic, finally, at the moment the most efficient mechanism in oxide perovskites is that realized in A-site driven ferroelectrics such as BiFeO 3 (and possibly BiMnO 3 and PbVO 3 ), where the presence of lone pairs in Bi 3+ or Pb 2+ plays a central role in the origin of ferroelectric polarization.
Such a mechanism, first theoretically proposed for BiMnO 3 [45] and then experimentally verified in the widely investigated BiFeO 3 [46][47][48][49], originates from the on-site sp rehybridization of the two 6s electrons of bismuth or lead that do not participate in chemical bonds, thus showing a high polarizability which results in a very large ferroelectric polarization (∼100 µC cm −2 [50] in BiFeO 3 ). On the other hand, magnetism is guaranteed by the B-site magnetic Fe 3+ (d 5 ) ions. Due to the independent origin of ferroelectricity and magnetism in this type of multiferroics, the magnetoelectric coupling has been long thought to be rather small; however, a number of recent observations indicate that magnetoelectric effects in BiFeO 3 can be significant, suggesting also an important role played by the large ferroelectric polarization in modulating the magnetic ordering [51][52][53][54]. These recent findings suggest potential relevance also in this class of proper multiferroics, where the inherently giant ferroelectric polarization together with unusual magnetoelectric effects can find interesting technological applications.
It should be noted, however, that the ferroelectric nature of BiMnO 3 has been the subject of controversial studies. After the initial claim of lone-pair ferroelectricity [45], computational structural optimizations revealed that the ferroelectric phase is not stable and indeed converges to the centrosymmetric C2/c structure [55]. This result was later confirmed by the experimental structural data of Belik et al [56]. Subsequent studies tried to reconcile these two seemingly conflicting pictures by proposing that, even considering a centrosymmetric crystal structure, BiMnO 3 can behave as an improper ferroelectric (i.e. where ferroelectricity is induced by the magnetic order, which breaks the inversion symmetry) [57]. Moreover, Solovyev and Pchelkina [57] suggests that BiMnO 3 should be a rare example of materials, where ferroelectricity indeed coexists with (canted) ferromagnetism, and demonstrates how the polarization can be controlled by the magnetic field.
In this latter context of A-site-driven ferroelectrics with large polarization, we can place a relatively new multiferroic, PbNiO 3 , which has recently attracted much interest and which will be the focus of the remaining part of the paper. By means of density-functional simulations, complemented by state-of-the-art hybrid functionals and DFT+U approaches, we have explored the structural, electronic and magnetic properties of two different phases of PbNiO 3 , i.e. the orthorhombic and the LiNbO 3 -like rhombohedral polymorphs. The outcomes of our simulations allow for a thorough understanding of the PbNiO 3 electronic structure, in turn leading to a microscopic understanding of the origin of its ferroelectric polarization, pressure-induced phase transition and magnetism, thereby constituting an additional example of how ab initio simulations can be useful in the framework of new multiferroics.

Introduction
The perovskite structure is ubiquitous in materials science, especially among oxide materials. They show very different properties, like superconductivity, colossal magnetoresistance, ionic conductivity and dielectric properties, which are of great importance in numerous technological applications [58,59]. The perovskite structure, with general formula AB X 3 , gives the possibility to accommodate a large variety of A and B cations with different A−X and B−X bond lengths. Due to this great flexibility, it is possible to realize different types of distortions from the ideal cubic structure. These include tilting of the octahedra, off-centering displacements of the cations out of the coordination polyhedra, and deformations of the octahedra driven by electronic factors (i.e. Jahn-Teller distortions) [23,60,61]. The electronic, magnetic and dielectric properties of perovskites depend crucially on the details of these distortions [58,59]. A useful structural related parameter, the so-called Goldschmidt tolerance factor t, can be used to predict the possible low-symmetry structural distortion, starting from the ideal AB X 3 cubic structure. It is defined as where r A , r B and r O denote the ionic radii of A, B and oxygen, respectively [62]. A range of t values approximately between 0.9 and 1.0 corresponds to stable cubic perovskite structures. For t > 1, the equilibrium structure shows face sharing of BO 6 octahedra, resulting in hexagonal perovskites. On the other hand, for t < 0.9 cooperative rotations of the octahedra yield lower symmetry perovskites such as tetragonal, orthorhombic, rhombohedral, etc [62]. External perturbations, like temperature and/or pressure, can be used to induce phase transitions giving access to other metastable phases. The change of octahedral distortions have been observed in many perovskite-related systems upon heating or cooling, and application of pressure and alloying (i.e. chemical pressure). For example, recently Inaguma et al have successfully synthesized two polymorphs of PbNiO 3 and CdPbO 3 , the orthorhombic (space group Pnma) and the LiNbO 3 -type rhombohedral structure (space group R3c), both shown in figure 1. In situ energy dispersive x-ray diffraction experiments revealed that a phase transition between the two phases could be realized irreversibly in PbNiO 3 if the appropriate thermodynamic conditions, 400 • C and 3 GPa are established [63]. The structural phase transition in CdPbO 3 strongly depends on heating temperature and there are temperature ranges in which both phases are seen to coexist [64]. Additionally, Goldschmit diagram analysis [63,64] suggests a boundary between the orthorhombic and rhombohedral phase in the vicinity of t = 0.85. On the basis of this argument, it has been predicted that several compounds could be stabilized in both phases under different pressure conditions [64]. Both phases show interesting electronic and magnetic properties [63][64][65][66]. The magnetic susceptibility measurements showed that the orthorhombic and rhombohedral phases in PbNiO 3 undergo an antiferromagnetic transition at 225 and 205 K, respectively. In addition, both phases have semiconducting behavior. Bond valence sum based on the refined structures and XPS data revealed that the valence state in both phases should be Pb 4+ Ni 2+ O 3 [63], but DFT-based calculations have clarified that this formal valence view is perturbed by the strong Pb 6s-O 2p hybridization which in turn results in enhanced ferroelectric instabilities [65]. Of particular interest is that the non-centrosymmetric rhombohedral structure allows a relative displacement of Pb, Ni and O atoms along the threefold [111]-axis, i.e. out of the centers of their coordination polyhedra, thus inducing a spontaneous polarization, which put forward the LiNbO 3 -type PbNiO 3 as a potential candidate for proper multiferroic behavior near room temperature.

Methods
We performed first-principles calculations using the Vienna ab initio Simulation Package (VASP) [67]. All the results were obtained using the projector-augmented plane-wave method [68] by explicitly treating 14 valence electrons for Pb (5d 10 6s 2 6p 2 ), 10 for Ni (3d 8 4s 2 ) and 6 for oxygen (2s 2 2p 4 ). Spin-orbit coupling was not included. Brillouin zone integrations were performed using a Gaussian broadening [69] of 0.1 eV during all relaxations, and the ions were relaxed until the Hellmann-Feynman forces were less than 10 −2 eV Å −1 . We used a 8 × 6 × 8 and 8 × 8 × 8 Monkhorst-Pack [70] k-point mesh centered at for orthorhombic and rhombohedral PbNiO 3 , respectively. We used a large basis set defined by plane-wave cutoff of 600 eV. These values ensure good convergence for all the investigated properties.
It is well known that local and semi-local density functionals (generalized gradient approximation (GGA)) [71] usually underestimate the size of the band gap in 'strongly' correlated systems, e.g. those containing localized d orbitals, and they even incorrectly predict metallic behavior for materials that are known to be insulators. Therefore, in order to deal with these shortcomings, we used the GGA+U method [72,73] where the strong Coulomb repulsion between localized d states is treated by adding a Hubbardlike term to the effective potential. This usually gives an improved description of correlation effects in magnetic transition-metal oxides [74][75][76][77]. The GGA+U method requires two parameters, the Hubbard parameter U and the exchange interaction J.
In this work, we adopted the rotationally invariant approach in which the parameters U and J do not enter separately, but rather as an effective Hubbard parameter given by their difference [78,79]. The on-site correction was set to be 4.6 eV for Ni 2+ 3d electrons adapted from previous self-consistent calculations [74]. However, we tested that the results are not strongly dependent on the choice of the U eff parameter chosen within a reasonable range of U values. As a further check we have also applied the Heyd-Scuseria-Ernzerhof (HSE) screened hybrid functional [80]. In the HSE method the short-range (sr) part of the exchange interaction (X) is constructed by a mixing of exact non-local Hartree-Fock exchange and approximated semi-local Perdew-Burke-Ernzerhof (PBE) exchange. The remaining contributions to the exchange-correlation energy, namely the long-range (lr) exchange interaction and the electronic correlation (C), are treated at the PBE level only, as defined in the following: The partitioning between sr and lr interactions is achieved by properly decomposing the Coulomb kernel (1/r , with r = |r − r |) through the characteristic parameter µ, which controls the range separation between the short (S) and long (L) range part We have used here µ = 0.20 Å −1 and α = 0.25, which is known to provide a very accurate description of the structural properties of PbNiO 3 [65]. A more detailed discussion of the role of α on the ground state properties of magnetic oxide perovskites can be found in [81]. In the last few years, the application of hybrid functionals to a wide class of solid state systems, including ferroelectrics [82][83][84], multiferroics [16,65,66] and magnetic perovskites [81,85] has shown that an improved description of ground state structural, electronic and magnetic properties with respect to the standard GGA and GGA+U can be achieved, in general. Due to the high computational cost of HSE calculations, the Fock exchange was sampled using the twofold reduced k-point grid.
Phonon dispersions were calculated within the density functional perturbation theory using the code PHONOPY [86].
In order to determine the magnetic ground state, we considered ferromagnetic, A-type, C-type as well as G-type antiferromagnetic orderings. The calculations show that the G-type configuration is the magnetic ground state for both the systems considered hereafter. Therefore, we will focus on the G-type ordering from now on, unless otherwise specified.

Results and discussion
2.3.1. Equilibrium structures and phase stability. As reported in our previous work [65], the GGA+U and HSE schemes are essential to reproduce the structural properties of rhombohedral PbNiO 3 : the HSE and, to a lesser extent, GGA+U reproduce the structural properties in good agreement with experiments, while the GGA fails to do so. The collection of the structural data for the LiNbO 3 -type rhombohedral phase is proposed again in table 1.
Results for the high pressure orthorhombic phase are shown in table 2. The GGA functional could not reproduce the atomic positions accurately, giving rise to incorrect bond lengths and angles, although there seems to be good agreement for the lattice parameters. The inclusion of either the on-site U or the Fock exchange improves the GGA description.
To explore the phase stability, the structures were fully relaxed for all volumes using force as well as stress minimization within GGA and GGA+U methods. The total energy versus volume is shown in figure 2. Note that the GGA method (figure 2(a)) predicts the orthorhombic Pnma phase as the ground-state structure for PbNiO 3 , which is in contrast with experiments [63]. On the other hand, GGA+U (figures 2(b), (c)) predicts the rhombohedral R3c phase as the ground state for PbNiO 3 although with a slightly larger volume, and correctly describes the relative stability of these two phases: the predicted GGA+U transition pressure is P≈5 GPa, slightly larger than the experimental value, ≈3 GPa [63] (the deviation can be caused by temperature effects, which are not taken into account in the present work), whereas the calculated volume collapse, ≈4.5%, is in agreement with the reported measured value of ≈3.5%.
To investigate the possible role of phonon instabilities in the observed structural transition from the low pressure R3c phase to the high pressure Pnma phase, we have computed the phonon spectra of both structures at their equilibrium volume (V R3c = 244.32Å 3 , V Pnma = 233.39Å 3 ) as well as at compressed (R3c, V = 244 Å 3 ) and expanded (Pnma, V = 250 Å 3 ) volumes. The results shown in figure 3 indicate that none of the considered phases exhibit negative modes in the off-equilibrium regime, suggesting that the pressure-induced structural transition is not phonon-mediated, but rather connected with electronic modifications of the chemical bonds. This implies that the R3c/Pnma phases are metastable structures at high/low pressures. Due to the high computational workload and the good performance of the GGA+U method, we did not perform a similar study with the HSE approach. opposite spin direction, is the most stable configuration, with an energy gain with respect to ferromagnetic, C-type antiferromagnetic and A-type antiferromagnetic configurations of 55.7, 25.6 and 32.7 meV per formula unit, respectively. By mapping these energies onto a Heisenberg model Hamiltonian H = 1 2 i j J i j S i S j with normalized spins, we can estimate the nearest-neighbor exchange couplings J and J ⊥ , parallel and perpendicular to the ac-plane, respectively, as well as the next-nearest neighbor isotropic (additional-and expected to be smaller-anisotropy of possible next-nearest neighbor interactions has been neglected) exchange coupling J 2 . We find a significant anisotropy of nearest-neighbor exchange couplings, namely J = 7.85 meV and J ⊥ = 12.15 meV, that can be most probably ascribed to the different buckling of Ni-O-Ni angles (by ≈ 3 • ) mediating the superexchange interaction between neighboring Ni ions. Furthermore, the next-nearest-neighbor exchange interaction is predicted  to be weakly ferromagnetic, J 2 = −0.16 meV. The critical temperature for the magnetic transition can be then predicted in a mean-field approach as T c ∼ 2 3 (2J + J ⊥ − 6J 2 ), which gives T c ∼ 223 K, in perfect but perhaps fortuitous agreement with the experimental value T c = 225 K. Indeed, the mean-field value is known to overestimate the exact solution. A more accurate classical Monte Carlo calculation on a 20 × 20 × 20 cell, using a Metropolis algorithm, predicts T c = 186 K, which is 20% lower than the experimental value, well within the typical error associated with ab initio predictions of critical temperatures [87]. Such underestimation is most probably due to the artifacts of the GGA+U approach [88].
With respect to the electronic properties of the orthorhombic phase, an interesting issue is related to the possible existence of lone-pairs on Pb ions. In order to evaluate the nominal oxidation state of the atoms, a valence bond sum based on the Rietveld structural refinement as well as XPS measurements were used [63], suggesting that the ionic configuration of PbNiO 3 is Pb 4+ , Ni 2+ , and O 2− . The orthorhombic PbNiO 3 was then reported as the first example of the perovskite containing a tetravalent A-site cation without lone pair electrons [63].
To gain more insights into the valence state in orthorhombic PbNiO 3 , the calculated total and partial density of states for G-type configuration within the corresponding optimized structure via GGA+U and HSE schemes are given in figure 4. The top of the valence band is set to zero, as indicated by a dotted vertical line. Note that the electronic structure obtained within GGA+U and HSE is very similar. The band centered at −7 eV is dominated by the Pb 6s orbital.  An energy overlap between the Pb 6s and O 2p orbitals suggests a covalent bond between them. The Ni 3d states are located starting at −6 eV to the Fermi level. The energy overlap of the Ni 3d and O 2p states in the whole valence band indicates the presence of a strong covalent bonding between these two ions. The electronic configuration of Ni ions is t 6 2g e 2 g , as can be deduced from the total occupation of the t 2g states in both spin channels and of the e g states in one spin-channel, and empty e g states in the other spin channel. This results in 2+ state of Ni atom and 1.66 µ B magnetic moment, very close to the nominal d 8 states. However, the valence state for the Pb ion cannot be determined easily. As shown before for the rhombohedral R3c phase [65], the significant hybridization between Pb 6s and O 2p orbitals suggests a reduction of the O valence state, i.e., the valence state of the Pb ion is smaller than the nominal 4+. Our analysis indicates that the simple concept of nominal valence charge is not sufficient to capture the subtle electronic properties of this system.
According to the density of states shown in the upper panel of figure 4, the GGA+U calculation results in a semiconducting behavior with a narrow band gap of 0.05 eV, which does not depend much on the chosen U . This is not unexpected since the band gap is between the bottom of Pb 6s conduction band and the top of O 2p valence states, which are states not affected by the Hubbard term.
Although the electronic structure obtained via GGA+U and HSE method are very similar, some differences can be found. In HSE, the spectral distribution of O 2s and Pb 5d states are shifted considerably downwards. Also, the bonding Pb 6s states are shifted downwards, while the antibonding states are shifted upwards, increasing the band gap to 1.05 eV, which is much larger than that of GGA+U . It has been argued that the downward shift of d states in transitionmetal oxides is an artifact of the LDA+U scheme [88].

The rhombohedral
LiNbO 3 -type phase. The low-pressure phase of PbNiO 3 has a rhombohedrally non-centrosymmetric structure with the R3c space group, which is isostructural with the most common multiferroic material BiFeO 3 [89]. More interestingly, the rhombohedral phase undergoes antiferromagnetic transition at 205 K, and exhibits semiconducting transport properties [63]. This therefore suggests a possible multiferroic behavior, i.e. co-existence of ferroelectric and magnetic properties. From the energy differences between the considered ferromagnetic, C-type and A-type configurations with respect to the G-type antiferromagnetic ground state, we could estimate the relevant exchange coupling constants J 1 = 8.69 meV, J 2 = 0.27 meV and J 3 = 0.19 meV (where J 3 represents now the isotropic exchange interaction between the third nearest neighboring Ni ions). The nearest-neighbor exchange interactions are now isotropic, consistently with the R3c symmetry of the rhombohedral structure. A small anisotropy could in principle be present in the next nearest neighbor exchange interactions, which is however unattainable with the considered magnetic configurations. More interestingly, however, the next nearest neighbor interaction J 2 is found to be weakly antiferromagnetic, implying a weak magnetic frustration that can well explain the lower critical temperature as compared to the orthorhombic phase. From our estimated J i j , the mean-field formula predicts T c ∼ 195 K (the MonteCarlo approach gives 143 K). Again, the predicted T c significantly underestimates the experimental value; however, the trend we find, i.e. the lowering of the magnetic Néel temperature upon structural transition, is in qualitative agreement with experimental observations, being related to both the weak magnetic frustration introduced by the antiferromagnetic J 2 as well as to a volume effect (responsible for a slight decrease of the nearest-neighbor exchange interaction). In the next section, we will mainly focus on the ferroelectric properties of rhombohedral PbNiO 3 and the apparently mysterious origin of a strong polarization induced by a nominally Pb 4+ ion.

Ferroelectric properties.
We evaluated the spontaneous polarization P of rhombohedral PbNiO 3 from first principles considering as centrosymmetric structural reference the R-3c symmetry, and we linearly interpolate the atomic positions between the centric and the polar phase, i.e. the so-called adiabatic path [11][12][13].
For the paraelectric phase, we used the same lattice constant and rhombohedral angle of the ferroelectric one, but we relaxed the internal atomic positions. It is worth mentioning that the GGA+U calculations were carried out with U eff = 7.6 eV using the relaxed structure obtained by U eff = 4.6 eV, since this value gives an insulating ground state for both polar and non-polar systems.
In general, it is useful to represent the spontaneous polarization P tot as a sum of the ionic contribution, P ion , and the electronic, P ele : P tot = P ion + P ele [90]. Although the ionic and electronic components are sometimes useful tools for qualitative discussions, one has to keep in mind that only the total polarization difference between two structural configurations has a physical meaning [15]. In table 3, we report the ionic and electronic contributions to the difference of ferroelectric polarization between the polar and non-polar configuration, both at GGA+U and HSE level. The computed total polarization P tot reaches 95.7 and 99.5 µC cm −2 at GGA+U and HSE level, respectively, comparable to that of BiFeO 3 [16,46]. The spontaneous polarization of LiNbO 3 and ZnSnO 3 was reported to be 80 µC cm −2 [91] and 56.9 µC cm −2 [92], respectively, which is much smaller than the present results for rhombohedral PbNiO 3 . As far as the relative stability is concerned, the energy barrier between the centric and the polar phase is ∼0.55 eV per formula unit for DFT+U and HSE, thus supporting the thermodynamic stability of the polar R3c phase with respect to the centrosymmetric R-3c phase, and indicating that the LiNbO 3 -type PbNiO 3 may be ferroelectric with a relatively high Curie temperature.
In order to separate purely electronic and structural effects on the calculation of the polarization, we performed the HSE calculations using the GGA+U geometry, which we denote as P HSE(GGA+U ) tot . Note that the ionic contribution is of course the same in both cases, while P GGA+U tot = 95.7 µC cm −2 and P HSE(GGA+U ) tot = 93.3 µC cm −2 . Therefore, keeping the volume unchanged and including the Fock exchange, the polarization decreases, while the band gap increases from 0.37 eV (GGA+U ) to 0.89 eV (HSE). This scenario agrees with the intuitive expectation that the larger the energy gap, the lower the polarization should be [16]. Additionally, it is sometimes useful to evaluate the polarization within the PCM. It is evaluated using the following equation where , Z v i and r i denote the unit-cell volume, the formal charge and the displacement of the ion i, respectively. A large difference between the PCM value and the P tot obtained from a fully quantum mechanical calculation signals that covalency effects are at play. Our PCM polarization value is 87.6 µC cm −2 and 95.5 µC cm −2 for GGA+U and HSE atomic structures, respectively. The closeness of HSE-structure PCM and full HSE values of polarization indicates that (i) the HSE description of PbNiO 3 is more ionic than that of GGA+U (equivalently, we can say that nominal and effective atomic charges are similar in HSE); (ii) LiNbO 3 -type PbNiO 3 should be a proper displacive ferroelectric. Note that even in the absence of highly anomalous charges, covalent effects can play an appreciable role in the determination of large polar r i displacements.
Density of states analysis. The comparison of the density of states between ferroelectric R3c and paraelectric R-3c phase [65] indicates that the major variation is the ∼2 eV downshift and the broadening of the Pb 6s-O 2p spectral weight from the paraelectric to the ferroelectric phase. This highlights the role of the Pb 6s-O 2p hybridization in stabilizing the non-centric ferroelectric phase. A similar broadening of the bonding Pb 6s-O 2p states was found in BiFeO 3 [46,93] where the stereochemical activity of the 6s lone pair of Bi ions is the driving force for the off-centering distortions.
Born effective charge analysis. We have calculated Born effective charge tensors using the linear-response formalism within GGA+U (U eff = 7.6 eV). Table 4 lists the Born effective charge tensors of Pb and Ni in PbNiO 3 in the ferroelectric and paraelectric phases. For Pb it was shown in [65] that its static charge (2.15) is much smaller than the nominal 4+; on the other hand, its Born effective charge is much larger than the static charge, indicating the presence of a large anomalous contribution due to covalent bonding between Pb and O ions. This situation is analogous to that of the highly polarized BiFeO 3 [46,93].
Note also the larger zz component of Pb effective charge that can be traced back to the anisotropy of bond lengths in the PbO 6 octahedra in the ferroelectric phase, and indicates that Charge density difference analysis. In figure 5(a) we plot the charge density difference ρ(r) = ρ ferro (r) -ρ para (r) between the ferroelectric R3c and paraelectric R-3c state.
In the paraelectric phase the electron density contour is symmetric along the [111]direction, while in the ferroelectric phase the lack of the inversion center allows shifts of Pb and O atoms along the [111]-direction, giving rise to an asymmetric distribution of the electron density. As shown in figure 5(a), the difference of charge density contours clearly shows a strong asymmetric electron re-distribution around the Pb atoms along the [111]-direction in the ferroelectric phase. Again, this supports an asymmetric covalent bonding interaction between the Pb 6s and O 2p orbitals.
ELF analysis. We use the ELF analysis, which is a useful tool for describing the chemical bonding character in extended systems [94]. ELF is a measure of the probability of finding an electron near another same-spin electron. It is a ground state property and can be used to discriminate between different kinds of bonding [94]. The ELF is represented as a contour plot in real space, spanning values between 0 and 1. A region where ELF is 1 means there is no possibility to find two electrons with the same spin. This usually occurs in regions occupied by bond pairs or lone pairs; ELF = 0 corresponds to vacuum, while in metals ELF∼0.5 [95]. In figure 5(b) we present the difference in ELF between the ferroelectric and paraelectric phases of PbNiO 3 , i.e. E L F(r) = ELF ferro (r) − ELF para (r). Thus, positive values of ELF(r) represent the area where the electron localization is higher, i.e. the covalent bonding is strengthened; negative values of ELF(r) represent the area where the electron localization is lower and the covalent bond weaker [90,96].
The ELF(r) isosurface clearly indicates that the ferroelectric distortion leads to a strongly asymmetric electron localization along the [111]-direction. The positive region mainly resides in the pyramid formed by the Pb atom and the three oxygen atoms above the Ni (0.5 0.5 0.5) atom; in contrast, the negative counterpart is predominately located in the pyramid formed by the Pb atom and the other three oxygen atoms beneath. This asymmetric distribution is associated with the asymmetric Pb and O atoms bonding interaction along the [111]-direction. The existence or absence of ferroelectricity is determined by a delicate balance between long-range Coulomb forces (which favor the ferroelectric state) and short-range repulsions which are minimized in the symmetric structure (favoring the non-polar paraelectric state [20]). The tendency of additional bonding between the Pb 6s and O 2p orbitals stabilizes the ferroelectric distortion.
All these results suggest that the asymmetric bonding interaction between the Pb 6s and O 2p orbitals due to the presence of the lone-pair on Pb is the primary driving force for ferroelectricity in rhombohedral PbNiO 3 . Clearly, a naïve picture just based on the assumption of nominal charges is unable to capture this complexity.

Summary and conclusions
We have reported structural, electronic and magnetic properties for orthorhombic and rhombohedral PbNiO 3 obtained from ab initio density functional calculations using GGA, DFT+U and HSE. Besides providing a detailed analysis of the pressure-driven R3c-Pnma structural transition, we have re-explored and discussed the ferroelectric properties of rhombohedral PbNiO 3 , which is characterized by a large spontaneous polarization, comparable to that of BiFeO 3 . Density of states, Born effective charge, charge density and ELF analysis have been used to study the ferroelectric transition.
Our main conclusions are: (1) The calculated structural parameters for rhombohedral PbNiO 3 via the GGA+U and HSE scheme are in agreement with experimental values, especially for the HSE method, while GGA (PBE) fails to reproduce the relevant properties. (2) A pressure induced structural transition from the rhombohedral to the orthorhombic structure occurs around 5 GPa at the GGA+U level. This scenario is in agreement with experimental observations. On the other hand, the GGA (PBE) scheme predicts the orthorhombic phase to be the ground state, in sharp contrast with the experiments. (3) The study of the magnetic interactions reveals a significant anisotropy of the nearestneighbor exchange couplings in the orthorhombic phase. This magnetic anisotropy is considerably attenuated (virtually vanishes) in the ferroelectric rhombohedral LiNbO 3 -type phase. The computed magnetic transition temperature is in line with the measured one.
(4) Detailed analysis of the density of states indicates that the valence state of the Pb ions is smaller than the nominal charge 4+, in contrast to naïve expectations based on the fact that Ni is found to be in a high spin 2+ state. (5) Our study confirms that rhombohedral PbNiO 3 might be a potential candidate for multiferroic applications. The spontaneous polarization is 100 µC cm −2 along the [111]direction, comparable to that of BiFeO 3 . (6) The asymmetric hybridization interaction between Pb 6s and O 2p orbitals appears as the driving force of atomic polar displacements along the [111]-direction, thus giving rise to a large spontaneous polarization.