## Abstract

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 BiCoO_{3}. 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 BiCoO_{3} via chemical doping to reproduce such marked effects under ambient conditions, thereby yielding useful multifunctionality.

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

## INTRODUCTION

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 BiFeO_{3} (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 () that is stable under ambient conditions displays displacements of the Bi cations and concerted antiphase rotations of the O_{6} octahedra about the polar [111] axis. (Axes are given in the pseudocubic setting.) There is also a paraelectric (PE) phase characterized by antiphase O_{6} tilts about [110] and in-phase rotations about [001]. This orthorhombic () 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*).

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 BiCoO_{3} (BCO)—also a room temperature multiferroic—complies with this requirement (*6*). Under ambient conditions, BCO presents an FE tetragonal () 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 phase with *c*/*a* ≈ 1 and a three-dimensional spin lattice; our corresponding first principles–based Heisenberg model yields *T*_{N} ≈ 500 K. Similar to the spin-spin interactions, we expect the SP couplings in BCO’s and phases to differ significantly as well.

## RESULTS

### SP-controlled transitions at ambient pressure

To test this, we compute the free energy of BCO’s and 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 K, in agreement with the experimental value K (Fig. 2A) (*7*). Note that the - 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 phase and that SP effects are critical to reproduce the experimental *T*_{c}.

To understand this, note how the Γ-phonon frequencies change when considering AFM and FM (ferromagnetic) spin orders in the and 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 phase, large and positive Δω’s mostly correspond to high-energy phonons (*ħ*ω ≥ 60 meV), whereas in the phase, those are associated to relatively low-frequency modes (*ħ*ω ~ 30 meV). Consequently, magnetic disorder favors the polymorph. At , for instance, fluctuating spins provide a lattice free energy difference of 0.168 meV/fu (formula unit) between and , which is three times larger than the one obtained when constraining AFM spin order.

### SP-controlled transitions under compression, reentrant behavior

In most perovskites, hydrostatic pressure (*P*) favors the phase over competing polymorphs (*7*, *9*). Hence, compression might help to reduce BCO’s *T*_{c} 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 - transition in the limit of low temperatures, GPa, is in fair agreement with the experimental value 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 *P*_{c}(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 (*10*–*13*)]. With regard to the critical pressure and volume drop for the - 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 phase followed, upon cooling, by a PM phase, a G-type AFM phase, and a C-type AFM 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) (*14*–*16*). 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 and phases display different Néel temperatures and SP couplings. In Fig. 4A, we show the temperature dependence of the quasi-harmonic Gibbs free energy, , of BCO’s polymorphs calculated at 2.5 GPa. We find that, whenever a phase becomes PM, the slope of the corresponding 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 (; Fig. 4B) and enthalpic (; Fig. 4C) terms, the latter being responsible for the slope changes accompanying the spin transitions. This effect, which is larger in the 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).

## DISCUSSION

### Engineering multiferroic effects under ambient conditions

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

Chemical substitution is a practical strategy to reduce BCO’s *T*_{c} at ambient pressure. As a simple predictor for *T*_{c}, we monitor the enthalpy difference between the and phases at zero temperature, Δ*H*_{eq}, which is fast to compute from first principles. We thus look for chemical substitutions that yield − 0.08 ≤ Δ*H*_{eq} ≤ − 0.03 eV/fu to match the results for pure BCO around 2.5 GPa. We find two potential cases—namely, BiCo_{1/2}Fe_{1/2}O_{3} and Bi_{3/4}La_{1/4}Co0_{3}—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 and 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.

## METHODS

To calculate the Gibbs free energy of BCO’s and 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 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 as(1)where is the effective static energy, is the effective force-constant matrix, and *u*_{m} and *u*_{n} 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, , is calculated by finding the eigenfrequencies of the dynamical matrix associated to , namely, , and plugging them into the formula(2)where *N*_{q} is the total number of wave vectors used for integration in the Brillouin zone. The Gibbs free energy of each phase is then estimated as(3)where the hydrostatic pressure is determined via the formula . Our Gibbs free energy results are accurate to within 5 meV/fu, and the corresponding transition points are determined through the condition .

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 and 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 phase, Cazorla and Íñiguez showed (*4*) that the quantities entering Eq. 1 can be calculated as(4)(5)where γ_{a}(*V*, *T*) ≡ < *S*_{i}*S*_{j} >/|*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 and correspond to(6)(7)(8)(9)

In the equations above, superscripts “FM” and “G” indicate perfect FM and AFM G-type spin arrangements, respectively. The 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 parameter captures the dependence of the phonon spectrum on the spin configuration (that is, SP coupling effects).

For the phase, we express the corresponding static energy and force constant matrix as(10)and(11)where γ_{α}(*V*, *T*) ≡ < *S*_{i}*S*_{j} >/|*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 and can be obtained as(12)(13)(14)(15)(16)(17)(18)(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 and 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 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 and phases, we constructed several spin Heisenberg models of the form , 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 phase) and AFM(G) (that is, in the phase) order parameters, namely, and . Here, *n*_{ix}, *n*_{iy}, and *n*_{iz} are the three integers locating the *i*-th lattice cell, and *N* is the total number of spins in the simulation box. For the calculation of *S*^{C} and *S*^{G}, we considered only the *z* component of the spins because a small symmetry-breaking magnetic anisotropy was introduced in the Hamiltonian to facilitate the analysis [see the supplementary materials in the study by Escorihuela-Sayalero *et al*. (*20*)].

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/3/6/e1700288/DC1

table S1. Calculated zero-temperature total energy differences between the tetragonal *P*4*mm*, 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 BiCoO_{3}.

fig. S2. Enthalpy energy differences between several crystal structures and the tetragonal *P*4*mm* 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 *P*4*mm*, 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 *P*4*mm* 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 *P*4*mm* 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 *P*4*mm* 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 *P*_{f} = 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 *P*_{f} = 2.5 GPa.

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

fig. S13. Spin free energy difference between the *P*4*mm* 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 “O_{3}”).

fig. S15. Phonon spectrum calculated in the Bi_{3}/_{4}La_{1}/_{4}CoO_{3} compound at zero pressure in the tetragonal *P*4*mm* phase.

fig. S16. Phonon spectrum calculated in the Bi_{3}/_{4}La_{1}/_{4}CoO_{3} compound at zero pressure in the orthorhombic *Pbnm* phase.

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.

## REFERENCES AND NOTES

**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.

- Copyright © 2017 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).