Multiple structural transitions driven by spin-phonon couplings in a perovskite oxide

See allHide authors and affiliations

Science Advances  30 Jun 2017:
Vol. 3, no. 6, e1700288
DOI: 10.1126/sciadv.1700288


Spin-phonon interactions are central to many interesting phenomena, ranging from superconductivity to magnetoelectric effects. However, they are believed to have a negligible influence on the structural behavior of most materials. For example, magnetic perovskite oxides often undergo structural transitions accompanied by magnetic signatures whose minuteness suggests that the underlying spin-phonon couplings are largely irrelevant. We present an exception to this rule, showing that novel effects can occur as a consequence. Our first-principles calculations reveal that spin-phonon interactions are essential to reproduce the experimental observations on the phase diagram of magnetoelectric multiferroic BiCoO3. Moreover, we predict that, under compression, these couplings lead to an unprecedented temperature-driven double-reentrant sequence of ferroelectric transitions. We propose how to modify BiCoO3 via chemical doping to reproduce such marked effects under ambient conditions, thereby yielding useful multifunctionality.

  • Multiferroics
  • spin-phonon coupling
  • density functional theory
  • phase transitions


Most ferroelectric (FE) and ferroelastic perovskite oxides undergo transitions involving structurally similar phases. One might guess that spin-phonon (SP) effects should play a role in these transformations, as it occurs in materials exhibiting more drastic changes (for example, Ni-based superalloys or steel) (1, 2). However, by excepting the especial case of compounds in which a magnetically driven symmetry breaking yields FE order (3), SP couplings tend to have no impact. Even in compounds like room temperature multiferroic BiFeO3 (BFO), in which SP effects significantly affect the free energy of competing phases, their influence on the structural transitions is minor (4).

Figure 1 shows the relevant polymorphs in BFO. The rhombohedral FE phase (Embedded Image) that is stable under ambient conditions displays displacements of the Bi cations and concerted antiphase rotations of the O6 octahedra about the polar [111] axis. (Axes are given in the pseudocubic setting.) There is also a paraelectric (PE) phase characterized by antiphase O6 tilts about [110] and in-phase rotations about [001]. This orthorhombic (Embedded Image) structure is stable above 1100 K (5). Magnetism in these phases is dominated by a strong antiferromagnetic (AFM) superexchange between adjacent irons, and first principles–derived spin models yield a Néel temperature of about 600 K for both of them (4). Furthermore, SP couplings turn out to be very similar in both structures and thereby have a minute impact on their relative stability (4).

Fig. 1 Structural and magnetic properties of energetically competitive polymorphs in bulk BiFeO3 and BiCoO3.

(A) rhombohedral R3c (R), (B) orthorhombic Pbnm (Embedded Image), and (C) tetragonal P4mm (Embedded Image). Patterns of O6 octahedra rotations are expressed in Glazer’s notation (30). The c/a aspect ratio of pseudocubic lattice parameters is approximately 1 for the Embedded Image and Embedded Image phases, whereas it is about 1.3 for the so-called super-tetragonal Embedded Image structure. Sketches of the lowest-energy spin configurations and exchange constants for a Heisenberg spin model of each phase are also shown.

We can conjecture that, for SP couplings to have a strong influence on the structural transitions, the magnetic interactions in the competing polymorphs need to be as different as possible. Perovskite BiCoO3 (BCO)—also a room temperature multiferroic—complies with this requirement (6). Under ambient conditions, BCO presents an FE tetragonal (Embedded Image) phase with polarization along [001] (Fig. 1C) and a very distorted cell with aspect ratio approaching 1.3. Consequently, the magnetic interactions within the [110] plane are stronger than those along the perpendicular direction, rendering a relatively low Néel temperature of about 310 K, according to our first-principles estimate. (See the Supplementary Materials for additional comments on this result and its comparison to the experiments.) At high temperatures, BCO presents a PE Embedded Image phase with c/a ≈ 1 and a three-dimensional spin lattice; our corresponding first principles–based Heisenberg model yields TN ≈ 500 K. Similar to the spin-spin interactions, we expect the SP couplings in BCO’s Embedded Image and Embedded Image phases to differ significantly as well.


SP-controlled transitions at ambient pressure

To test this, we compute the free energy of BCO’s Embedded Image and Embedded Image phases as a function of temperature, following the first-principles approach of the study by Cazorla and Íñiguez (4) (see Methods). We obtain a critical temperature at zero pressure of Embedded Image K, in agreement with the experimental value Embedded Image K (Fig. 2A) (7). Note that the Embedded Image-Embedded Image transition occurs at a temperature at which both phases are paramagnetic (PM) and that our method accounts for the contribution of disordered spins to the free energy. If the spins are frozen in their ground-state configuration, the transition is predicted to occur at 2350(50) K, which is unrealistically high (Fig. 2B). Hence, we find that magnetic disorder greatly contributes to the stabilization of the Embedded Image phase and that SP effects are critical to reproduce the experimental Tc.

Fig. 2 P-T phase diagram of bulk BCO calculated from first principles.

(A) SP coupling effects are considered in the calculation of quasi-harmonic free energies. Multiple T-induced multiferroic phase transitions occur in the region colored in yellow. (B) Fixed magnetic spin order (corresponding to the lowest-energy spin arrangement) is imposed in the calculation of quasi-harmonic free energies. Experimental data from Belik et al. and Oka et al. (6, 7) corresponding to FE-to-PE (dots) and AFM-to-PM (triangle) phase transitions are shown for comparison.

To understand this, note how the Γ-phonon frequencies change when considering AFM and FM (ferromagnetic) spin orders in the Embedded Image and Embedded Image phases. These frequency shifts, Δω ≡ ωAFM − ωFM, reflect the magnitude of SP couplings (8), and their sign indicates which phonon eigenmodes are more important to stabilize the corresponding PM phase (4). Figure 3 shows that in the Embedded Image phase, large and positive Δω’s mostly correspond to high-energy phonons (ħω ≥ 60 meV), whereas in the Embedded Image phase, those are associated to relatively low-frequency modes (ħω ~ 30 meV). Consequently, magnetic disorder favors the Embedded Image polymorph. At Embedded Image, for instance, fluctuating spins provide a lattice free energy difference of 0.168 meV/fu (formula unit) between Embedded Image and Embedded Image, which is three times larger than the one obtained when constraining AFM spin order.

Fig. 3 Analysis of SP couplings in the Embedded Image and Embedded Image phases of BCO.

Vibrational frequency shifts between AFM and FM spin configurations calculated at the reciprocal lattice point Γ, Δω = ωAFM − ωFM, are shown as a function of eigenmode energy (where this energy is the one obtained for the AFM ground state). Equivalent AFM and FM vibrational frequencies are unequivocally identified with the largest projection (scalar product) between AFM and FM Γ-phonon eigenmodes. Note that modes with a positive (negative) shift will tend to soften (harden) as T increases. Phonon frequency shifts for the Embedded Image phase span over a smaller energy interval than those for the Embedded Image phase, indicating that the former structure is vibrationally softer than the latter. In both phases, phonon eigenmodes presenting largest spin couplings correspond to medium- and high-energy excitations dominated by Co and O atoms; in contrast, low-energy eigenmodes dominated by Bi and O atoms, including those associated to the polar distortion in BCO, present small |Δω| values.

SP-controlled transitions under compression, reentrant behavior

In most perovskites, hydrostatic pressure (P) favors the Embedded Image phase over competing polymorphs (7, 9). Hence, compression might help to reduce BCO’s Tc and bring it closer to the AFM transition temperatures. To check this, we perform free energy calculations as a function of pressure (see Methods). Our results are shown in Fig. 2A.

Our prediction for the Embedded Image-Embedded Image transition in the limit of low temperatures, Embedded Image GPa, is in fair agreement with the experimental value Embedded Image GPa (7). Meanwhile, our first-principles calculations render a sequence of two spin-phase transitions induced by pressure at zero temperature: first, from a high-spin (HS) state to a mixed HS and low-spin (LS) (HS-LS) state, and second, from an HS-LS state to an LS state. However, the estimated transition pressures are well above Pc(0) (that is, ~ 50 GPa); hence, they are not relevant to the present discussion. [See the Supplementary Materials for a detailed discussion on our spin transition results in the light of measurements by Oka et al. (7) and for a comparison with respect to previous first-principles results obtained by other authors on bulk BCO under pressure (1013)]. With regard to the critical pressure and volume drop for the Embedded Image-Embedded Image transition at room temperature, we compute 2.55(0.15) GPa and ~11%, respectively, whereas the experimental values are 2.50(0.25) GPa and ~13%, respectively (7). Then, as shown in Fig. 2A, the agreement is less satisfactory at intermediate pressures. Finally, by comparing Fig. 2 (A and B), we ratify that SP effects are critical to reproduce the experiments.

Our phase diagram is rich in the region where structural and magnetic transitions get close. For P ≈ 2.5 GPa (colored area in Fig. 2A), we predict that BCO presents three temperature-driven transformations: a high-temperature PM Embedded Image phase followed, upon cooling, by a PM Embedded Image phase, a G-type AFM Embedded Image phase, and a C-type AFM Embedded Image phase. We move from a PE phase to an FE phase, back to a PE structure, and finally, to the FE ground state. Note that a PE-FE-PE sequence constitutes a rare reentrant behavior, because it is uncommon to stabilize a PE structure (typically more disordered) by cooling down an FE phase (which is typically more ordered) (1416). Here, we find a double reentrance, because the low-temperature PE phase eventually transforms into the FE ground state.

This unprecedented PE-FE-PE-FE sequence is possible because the Embedded Image and Embedded Image phases display different Néel temperatures and SP couplings. In Fig. 4A, we show the temperature dependence of the quasi-harmonic Gibbs free energy, Embedded Image, of BCO’s polymorphs calculated at 2.5 GPa. We find that, whenever a phase becomes PM, the slope of the corresponding Embedded Image curve changes noticeably; this results in three energy crossings (structural transitions) within an interval of about 325 K. The Gibbs free energy can be split into entropic (Embedded Image; Fig. 4B) and enthalpic (Embedded Image; Fig. 4C) terms, the latter being responsible for the slope changes accompanying the spin transitions. This effect, which is larger in the Embedded Image phase, corresponds to a sizeable decrease in the thermal expansion of the crystal when spins become disordered and is driven by SP couplings (namely, the effect disappears when considering frozen spins; see the Supplementary Materials).

Fig. 4 Calculated quasi-harmonic free energies of BCO’s competing polymorphs, as a function of temperature and at fixed pressure Pf = 2.5 GPa.

(A) Quasi-harmonic Gibbs free energy is calculated as Embedded Image; our estimates are accurate to within 5 meV/fu. (B) Quasi-harmonic Helmholtz free energy, Embedded Image, where Embedded Image represents the vibrational lattice entropy. Magnetic entropy effects stemming exclusively from the spin fluctuations have been safely neglected in our analysis, because their free energy difference among the Embedded Image and Embedded Image phases amounts to less than 5 meV/fu (see the Supplementary Materials). (C) Quasi-harmonic enthalpy, Embedded Image, where Embedded Image. Black arrows indicate the occurrence of structural transitions characterized by the thermodynamic condition Embedded Image. Black and red vertical lines signal magnetic spin order transformations occurring in the Embedded Image and Embedded Image phases, respectively. Black and red dashed lines in (B) and (C) represent results obtained by constraining AFM magnetic spin order in our quasi-harmonic free energy calculations, showing how spin disorder tends to stabilize the Embedded Image phase (B). Note that the temperature dependence of Embedded Image (B) is smooth; in contrast, the slope changes in Embedded Image (C), which are associated to the spin ordering transitions, are the main cause of the successive structural transformations.


Engineering multiferroic effects under ambient conditions

The phase diagram of Fig. 2A suggests interesting possibilities to obtain functional properties. For example, starting from the Embedded Image-AFM phase, one could use an electric field to induce the Embedded Image structure, which would result in either a loss of spin order (if we reach the Embedded Image-PM phase) or a transformation into a different AFM state (if we reach the Embedded Image-AFM phase with C-type order). For applications, one would like to realize such phase-change effects under ambient conditions (1719).

Chemical substitution is a practical strategy to reduce BCO’s Tc at ambient pressure. As a simple predictor for Tc, we monitor the enthalpy difference between the Embedded Image and Embedded Image phases at zero temperature, ΔHeq, which is fast to compute from first principles. We thus look for chemical substitutions that yield − 0.08 ≤ ΔHeq ≤ − 0.03 eV/fu to match the results for pure BCO around 2.5 GPa. We find two potential cases—namely, BiCo1/2Fe1/2O3 and Bi3/4La1/4Co03—for which the enthalpy differences (0.075 and −0.033 eV/fu, respectively) lie within the targeted interval. We find both compounds to be vibrationally stable; hence, they are good candidates to reproduce the marked effects predicted for compressed BCO under ambient conditions. (See the Supplementary Materials for more details.)

In conclusion, we have shown that SP effects are critical to understand and reproduce with theory the experimental phase diagram of magnetoelectric multiferroic BCO. SP couplings in the energetically competitive Embedded Image and Embedded Image phases are manifestly different, and consequently, the influence of T-induced spin disorder on the corresponding phase boundaries is tremendous. This is in stark contrast to what is observed in other multiferroic materials, for example, archetypal BFO, where SP effects are very similar in all the relevant polymorphs and therefore have a minor impact on the accompanying structural transformations. We predict that under moderate hydrostatic pressures of 2 ≤ P ≤ 3 GPa, BCO displays an unusual sequence of double-reentrant structural transitions caused by SP couplings. We argue that this series of multiple magnetic-FE transformations can be shifted down to ambient conditions by means of A- and B-cation substitutions involving La and Fe species, respectively.


To calculate the Gibbs free energy of BCO’s Embedded Image and Embedded Image phases as a function of temperature and pressure by considering SP couplings, we used the approach described in the study by Cazorla and Íñiguez (4) and generalized it to the Embedded Image phase along the guidelines described in the study by Escorihuela-Sayalero et al. (20). In particular, we started by expressing the internal energy of the crystal asEmbedded Image(1)where Embedded Image is the effective static energy, Embedded Image is the effective force-constant matrix, and um and un are the atomic displacements; the dependences of the various terms with V and T are explicitly noted. The Helmholtz free energy associated to the lattice vibrations, Embedded Image, is calculated by finding the eigenfrequencies of the dynamical matrix associated to Embedded Image, namely, Embedded Image, and plugging them into the formulaEmbedded Image(2)where Nq is the total number of wave vectors used for integration in the Brillouin zone. The Gibbs free energy of each phase is then estimated asEmbedded Image(3)where the hydrostatic pressure is determined via the formula Embedded Image. Our Gibbs free energy results are accurate to within 5 meV/fu, and the corresponding transition points are determined through the condition Embedded Image.

We note that, as described below, our formalism accounts for the way in which fluctuating spins affect the phonon frequencies; hence, the effect of SP couplings on the free energy is fully considered. In contrast, magnetic entropy contributions stemming exclusively from the spin fluctuations have been neglected in our analysis, because their free energy difference between the Embedded Image and Embedded Image phases can be estimated to be less than 5 meV/fu, that is, the accuracy threshold in our Gibbs free energy calculations (see the Supplementary Materials).

For the Embedded Image phase, Cazorla and Íñiguez showed (4) that the quantities entering Eq. 1 can be calculated asEmbedded Image(4)Embedded Image(5)where γa(V, T) ≡ < SiSj >/|S|2 represents the correlation function between neighboring spins and <…> the thermal average, as obtained from our Monte Carlo simulations of the corresponding spin Hamiltonian (see below). The rest of the parameters in Embedded Image and Embedded Image correspond toEmbedded Image(6)Embedded Image(7)Embedded Image(8)Embedded Image(9)

In the equations above, superscripts “FM” and “G” indicate perfect FM and AFM G-type spin arrangements, respectively. The Embedded Image parameter describes the magnetic interactions when the atoms remain frozen at their equilibrium positions (see Fig. 1B); typically, this captures the bulk of the exchange couplings. Meanwhile, the Embedded Image parameter captures the dependence of the phonon spectrum on the spin configuration (that is, SP coupling effects).

For the Embedded Image phase, we express the corresponding static energy and force constant matrix asEmbedded Image(10)andEmbedded Image(11)where γα(V, T) ≡ < SiSj >/|S|2, with α = a, b, ac, represents the correlation functions between in-plane and out-of-plane neighboring spins, according to the sketch shown in Fig. 1C. The rest of the parameters in Embedded Image and Embedded Image can be obtained asEmbedded Image(12)Embedded Image(13)Embedded Image(14)Embedded Image(15)Embedded Image(16)Embedded Image(17)Embedded Image(18)Embedded Image(19)

In the equations above, superscripts FM, G, “A,” and “C” indicate perfect FM, AFM G-type, AFM A-type, and AFM C-type spin arrangements, respectively. Next, we provided the details of our density functional theory (DFT) calculations (namely, zero-temperature energies and phonon spectra) and Heisenberg spin model Monte Carlo simulations, which are necessary to determine the value of the quantities entering Eqs. 1 to 19.

DFT calculations

We used the generalized gradient approximation to DFT proposed by Perdew et al. (GGA-PBE) (21), as implemented in the VASP package (22, 23). We work with GGA-PBE because this is the DFT variant that provides a more accurate description of the relative stability between the Embedded Image and Embedded Image phases of BCO at zero temperature, as discussed in the Supplementary Materials. A “Hubbard-U” scheme with U = 6 eV was used for a better treatment of Co’s 3d electrons (24). We use the “projector augmented wave” method to represent the ionic cores (25), considering the following electrons as valence states: Co’s 3p, 3d, and 4s; Bi’s 5d, 6s, and 6p; and O’s 2s and 2p. Wave functions are represented in a plane-wave basis truncated at 500 eV. We used a 20-atom simulation cell that can be viewed as a Embedded Image repetition of the elemental five-atom perovskite unit and that is compatible with all the crystal structures of interest here. For integrations within the Brillouin zone, we used Γ-centered q-point grids of 6 × 6 × 6. Using these parameters, we obtained enthalpy energies converged to within 0.5 meV/fu. Geometry relaxations were performed using a conjugate-gradient algorithm that keeps the volume of the unit cell fixed while permitting variations of its shape and atomic positions. The relaxation stops when residual forces fall below 0.01 eV Å−1. Equilibrium volumes were subsequently determined by fitting the sets of calculated energy points to Birch-Murnaghan equations of state (26). To treat the chemical substitutions, we work with a 40-atom cell that can be viewed as a 2 × 2 × 2 repetition of the elemental perovskite cell.

Phonon spectrum calculations

The calculation of phonon frequencies was performed with the direct method (27, 28), in which the force-constant matrix was calculated in real space by considering the proportionality between atomic displacements and forces when the former were sufficiently small. Large supercells need to be constructed to guarantee that the elements of the force-constant matrix have all fallen off to negligible values at their boundaries, a condition that follows from the use of periodic boundary conditions (29). Once the force-constant matrix is calculated, one can Fourier-transform it to obtain the phonon spectrum at any q-point. The impact of long-range interactions on the calculation of long-wavelength phonons is disregarded because we are primarily interested in the computation of quasi-harmonic free energies, and in this context, this factor is known to be secondary (4). The quantities with respect to which our phonon calculations need to be converged are the size of the supercell, the size of the atomic displacements, and the numerical accuracy in the sampling of the Brillouin zone. We find the following settings to provide quasi-harmonic free energies converged to within 5 meV/fu: 160-atom supercells that can be viewed as a 2 × 2 × 2 multiple of the 20-atom unit mentioned above, atomic displacements of 0.02 Å˚, and q-point grids of 12 × 12 × 12. The value of the phonon frequencies and quasi-harmonic free energies was obtained with the PHON code developed by Alfè (29). In using this code, we exploited the translational invariance of the system to impose the three acoustic branches to be exactly zero at the Γ q-point and use the central differences in the atomic forces (that is, positive and negative atomic displacements were considered).

Heisenberg model: Monte Carlo simulations

To simulate the effects of thermal excitations on the magnetic order of the Embedded Image and Embedded Image phases, we constructed several spin Heisenberg models of the form Embedded Image, in which the value of the involved exchange constants is obtained from zero-temperature DFT calculations (see discussion in the Supplementary Materials). We used these models to perform Monte Carlo simulations in a periodically repeated simulation box of 20 × 20 × 20 spins; thermal averages were computed from runs of 50,000 Monte Carlo sweeps after equilibration. These simulations allowed us to monitor the T dependence of the magnetic order through the computation of the AFM(C) (that is, in the Embedded Image phase) and AFM(G) (that is, in the Embedded Image phase) order parameters, namely, Embedded Image and Embedded Image. Here, nix, niy, and niz 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 SC and SG, we considered only the z component of the spins because a small symmetry-breaking magnetic anisotropy was introduced in the Hamiltonian to facilitate the analysis [see the supplementary materials in the study by Escorihuela-Sayalero et al. (20)].


Supplementary material for this article is available at

table S1. Calculated zero-temperature total energy differences between the tetragonal P4mm, orthorhombic Pbnm, and monoclinic Pc phases at zero pressure using several DFT exchange-correlation functionals.

I. Supplementary discussion: The role of the exchange correlation energy functional

II. Supplementary discussion: Spin-phase transitions under pressure

III. Supplementary data

fig. S1. Spin properties and spin-phase transitions under pressure calculated in BiCoO3.

fig. S2. Enthalpy energy differences between several crystal structures and the tetragonal P4mm phase calculated at T = 0 K with the PBE+U functional and expressed as a function of pressure.

fig. S3. Quasi-harmonic Helmholtz free energy calculated in the tetragonal P4mm, orthorhombic Pbnm, and monoclinic Pc phases at zero pressure; the corresponding magnetic spin orders are AFM(C), AFM(G), and AFM(G), respectively.

fig. S4. Analysis of SP coupling in the tetragonal P4mm phase at the special point Γ.

fig. S5. Analysis of SP coupling in the orthorhombic Pbnm phase at the special point Γ.

fig. S6. Analysis of SP coupling in the tetragonal P4mm phase over the corresponding Brillouin zone.

fig. S7. Analysis of SP coupling in the orthogonal Pbnm phase over the corresponding Brillouin zone.

fig. S8. AFM(C) magnetic spin order parameter calculated in the tetragonal P4mm phase as a function of pressure and temperature.

fig. S9. AFM(G) magnetic spin order parameter calculated in the orthorhombic Pbnm phase as a function of pressure and temperature.

fig. S10. Volume per formula unit calculated in BCO’s polymorphs at a fixed pressure of Pf = 2.5 GPa and expressed as a function of temperature.

fig. S11. Calculated quasi-harmonic free energies of BCO’s polymorphs at a fixed pressure of Pf = 2.5 GPa.

fig. S12. Contributions to the Gibbs free energy difference, Embedded Image, calculated in bulk BCO as a function of temperature and pressure.

fig. S13. Spin free energy difference between the P4mm and Pbnm phases calculated at P = 2.5 GPa as a function of temperature.

fig. S14. Chemical doping strategy proposed to bring the multiferroic phase competition disclosed in bulk BCO down to zero pressure (oxygen atoms, not indicated in the figure, always enter the formulae as “O3”).

fig. S15. Phonon spectrum calculated in the Bi3/4La1/4CoO3 compound at zero pressure in the tetragonal P4mm phase.

fig. S16. Phonon spectrum calculated in the Bi3/4La1/4CoO3 compound at zero pressure in the orthorhombic Pbnm phase.

References (3136)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: Funding: This research was supported under the Australian Research Council’s Future Fellowship funding scheme (project number FT140100135), the Israel Science Foundation through grants 1814/14 and 2143/14, and the Luxembourg National Research Fund through the programs PEARL (grant P12/4853155 COFERMAT) and CORE (grant C15/MS/10458889 NEWALLS). Computational resources and technical assistance were provided by RES (Red Española de Supercomputación) and by the Australian Government and the Government of Western Australia through Magnus under the National Computational Merit Allocation Scheme and the Pawsey Supercomputing Centre. Author contributions: All authors contributed equally to the present work. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article