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

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.

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) (, ). However, by excepting the especial case of compounds in which a magnetically driven symmetry breaking yields FE order (), 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 ().

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 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 () structure is stable above 1100 K (). 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 (). Furthermore, SP couplings turn out to be very similar in both structures and thereby have a minute impact on their relative stability ().

Fig. 1

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

(A) rhombohedral R3c (R), (B) orthorhombic Pbnm (), and (C) tetragonal P4mm (). Patterns of O6 octahedra rotations are expressed in Glazer’s notation (). The c/a aspect ratio of pseudocubic lattice parameters is approximately 1 for the and phases, whereas it is about 1.3 for the so-called super-tetragonal 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 (). 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 TN ≈ 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 () (see Methods). We obtain a critical temperature at zero pressure of K, in agreement with the experimental value K (Fig. 2A) (). 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 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. (, ) 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 and phases. These frequency shifts, Δω ≡ ωAFM − ωFM, reflect the magnitude of SP couplings (), and their sign indicates which phonon eigenmodes are more important to stabilize the corresponding PM phase (). 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.

Fig. 3

Analysis of SP couplings in the and 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 phase span over a smaller energy interval than those for the 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 phase over competing polymorphs (, ). 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 - transition in the limit of low temperatures, GPa, is in fair agreement with the experimental value GPa (). 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. () and for a comparison with respect to previous first-principles results obtained by other authors on bulk BCO under pressure (–)]. 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 (). 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) (–). 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).

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 ; our estimates are accurate to within 5 meV/fu. (B) Quasi-harmonic Helmholtz free energy, , where 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 and phases amounts to less than 5 meV/fu (see the Supplementary Materials). (C) Quasi-harmonic enthalpy, , where . Black arrows indicate the occurrence of structural transitions characterized by the thermodynamic condition . Black and red vertical lines signal magnetic spin order transformations occurring in the and 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 phase (B). Note that the temperature dependence of (B) is smooth; in contrast, the slope changes in (C), which are associated to the spin ordering transitions, are the main cause of the successive structural transformations.

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 (–).

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 and 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 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 () and generalized it to the phase along the guidelines described in the study by Escorihuela-Sayalero et al. (). In particular, we started by expressing the internal energy of the crystal aswhere is the effective static energy, is the effective force-constant matrix, and u and u 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 formulawhere N is the total number of wave vectors used for integration in the Brillouin zone. The Gibbs free energy of each phase is then estimated aswhere 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 () that the quantities entering Eq. 1 can be calculated aswhere γ(V, T) ≡ < SS >/|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

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 asandwhere γα(V, T) ≡ < SS >/|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

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) (), as implemented in the VASP package (, ). 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 (). We use the “projector augmented wave” method to represent the ionic cores (), 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 (). 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 (, ), 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 (). 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 (). 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è (). 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, n, and n 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. ()].