Magnetoelectric properties of the multiferroic CuCrO studied by means of ab initio calculations and Monte Carlo simulations
Motivated by the discovery of multiferroicity in the geometrically frustrated triangular antiferromagnet CuCrO below its Néel temperature , we investigate its magnetic and ferroelectric properties using ab initio calculations and Monte Carlo simulations. Exchange interactions up to the third nearest neighbors in the plane, inter-layer interaction and single ion anisotropy constants in CuCrO are estimated by series of density functional theory calculations. In particular, our results evidence a hard axis along the  direction due to the lattice distortion that takes place along this direction below . Our Monte Carlo simulations indicate that the system possesses a Néel temperature K very close to the ones reported experimentally ( K). Also we show that the ground state is a proper-screw magnetic configuration with an incommensurate propagation vector pointing along the  direction. Moreover, our work reports the emergence of spin helicity below which leads to ferroelectricity in the extended inverse Dzyaloshinskii-Moriya model. We confirm the electric control of spin helicity by simulating - hysteresis loops at various temperatures.
Through the discovery of the mineral CuFeO in 1873, Friedel opened the door to the delafossites ABO (1); (2). Such a family crystallizes in the layered space group, see Fig. 1. The diversity of properties they exhibit raises up an ever increasing interest in this class of compounds. In particular, the discovery of simultaneous transparency and -type conductivity in CuAlO by Kawazoe et al. (3), laid ground for the development of transparent optoelectronic devices. Furthermore, depending on the chemical composition, a plethora of behaviors can be evidenced. For instance, for A in a configuration, e.g., A = Pd or Pt, highly metallic compounds with anomalous temperature dependence of the resistivity have been reported (4); (5); (6); (7). The transport in these compounds has been found to be strongly anisotropic, with a degree of anisotropy that may reach 1000 (4); (5); (8). For A in a configuration, the semi-conducting materials CuBO, with B = Cr, Fe, Rh, may be turned into promising thermoelectric ones through hole doping (9); (10); (11) — in particular, an especially high power factor has been found in the case of CuRhMgO (12), which transport coefficients served as a basis for the Apparent Fermi Liquid scenario (13). Regarding the magnetic compounds CuFeO and CuCrO, many studies point towards a strong coupling of the magnetic and structural degrees of freedom (14); (15); (16); (17); (18); (19); (20); (21); (22), that paves the way to multiferroelectricity.
With its frustrated triangular lattice CuCrO received a lot of attention since it is ferroelectric without applying magnetic fields or doping upon Cr sites, unlike CuFeO (17); (24). The emergence of ferroelectricity in CuCrO is induced by the proper-screw magnetic ordering below the Néel temperature , and the control of this ferroelectricity by an applied magnetic field is very important for new spin-based device applications. CuCrO forms a rhombohedral lattice where the edge-shared CrO layers are alternatively stacked between Cu layers along the c-axis as shown in Fig. 1. Due to the weak inter-layer interaction (Fig. 2), the material behaves as a quasi-2D magnet, which makes it even more interesting.
The magnetic properties of CuCrO have been investigated by neutron diffraction experiments (25); (20); (26); (27); (28). It was shown that the magnetic configuration of CuCrO below is proper screw with an incommensurate propagation vector q = (0.329, 0.329, 0) (28) pointing along the  direction. Such deviation from the commensurate magnetic configuration of q = (1/3, 1/3, 0) is due to the lattice distortion that takes place along the  direction below upon the spiral-spin ordering which leads to anisotropic in-plane exchange interactions and (Fig. 2) (29). Polarized neutron-diffraction measurements on single crystals of CuCrO (20) showed that the spins are oriented in a spiral plane parallel to the (110) plane suggesting that the  direction is a hard axis.
The electric polarization emerges upon the spiral-spin ordering (30); (20); (31), which reflects the strong coupling between non-collinear magnetic ordering and ferroelectricity in CuCrO. Within the spin-current model or the inverse Dzyaloshinskii-Moriya (DM) mechanism (32); (33); (34), the electric polarization produced between the canted spins and , located at sites and , respectively, is given by
where is a unit vector joining the sites and . However, Eq.(1) fails to explain the emergence of ferroelectricity in CuCrO because in the proper-screw configurations, is parallel to ( is along the  direction due to symmetry considerations (30)) unlike the cycloid spin structures.
Based on symmetry considerations, Kaplan and Mahanti (35) introduced an additional contribution to the macroscopic polarization which contributes in both cycloid and proper-screw configurations. Therefore, within this model, now referred to as extended DM model, the total polarization is given by
In this study, we investigate the magnetoelectric properties of CuCrO by means of a combination of Density Functional Theory (DFT) calculations and Monte Carlo (MC) simulations. More precisely, we estimate a set of exchange interactions and anisotropy constants and confront it to the experimental magnetic properties and we verify the appearance of spiral spin ordering at low temperatures which can be related to the ferroelectric polarization.
In Sec. II we detail briefly the DFT method that we used to extract the coupling and anisotropy constants in CuCrO, while the model and MC method are presented in Sec. III. Sec. IV is devoted to the results where we discuss the magnetic and ferroelectric properties of CuCrO. A conclusion is given in Sec. V.
Ii DFT computational method
We performed a series of DFT calculations using full-potential linear muffin-tin orbital (FP-LMTO) method as implemented in RSPt (36) code. An experimental crystal structure (37) was considered, taking into account a small in-plane lattic distortion, suggested in Ref. (29). Our results are in-line with earlier calculations (21). The DFT+ (38) approach was used in order to take into account the effect of strong correlations between Cr electrons. The adopted values of Hubbard and Hundâs exchange were 2.3 and 0.96 eV, which were extracted from first-principles calculations for a similar system LiCrO (23). The same computational scheme was used in a prior study on the magnetic properties of CuCrO (22). The Fully Localized Limit (FLL) (39) formed of the double-counting correction was applied. We calculated the exchange parameters between Cr ions by means of the magnetic force theorem (40); (41). The so-called muffin-tin head projection scheme was applied to construct the set of localized Cr- orbitals (for more details see Ref. (42)). The ’s were extracted from both ferromagnetic and antiferromagnetic configurations. The obtained values turned out to be insensitive to the assumed magnetic order, which implies that they can be used as fixed parameters in a Hamiltonian describing the interacting spins. The spin-orbit coupling was taken into account only for the calculation of the magnetocrystalline anisotropy, which was calculated directly from the total energies.
Iii Model and Monte Carlo simulation
To model the magnetic properties of CuCrO, we note that Cr ions with = 3/2 spins are large enough to be treated classically, so we used the following classical three dimensional (3D) Heisenberg Hamiltonian
where refers to the exchange interactions up to the 4 neighbors (Fig. 2). The -axis corresponds to the  direction and the -axis corresponds to the  direction. and correspond to the hard and easy axes anisotropy constants respectively. The fourth term corresponds to the Zeeman energy where is the applied magnetic field ( is the Bohr magneton and is the Landé factor).
To model the ferroelectric properties of CuCrO and the coupling between the spins and the electric field , we added the following term to the previous Hamiltonian
where the sum runs over the magnetic bonds along the  direction, and is a coupling constant related to the spin-orbit and spin exchange interactions. Adding this contribution leads to the model for multiferroics proposed by Kaplan and Mahanti (35).
Our MC simulations (43) were performed on 3D triangular lattices (Fig. 1 with only Cr ions) with periodic boundary conditions (PBC) using the standard Metropolis algorithm (44) and the time-step-quantified method (45) when needed.
Typically, the first MC steps were discarded for thermal equilibration before averaging over the next MC steps. Note that our results are averaged over 24 simulations with different random number sequences so that statistical fluctuations are negligible.
Iv Results and discussions
It was reported in Ref. (29) that the lattice undergoes a tiny in-plane distortion below with and being the lattice constants along the  and the  directions, respectively. As a first step, we considered (29) to calculate the exchange interactions and anisotropy constants in CuCrO. The extracted values given in Table 1 (line 1) are very close to the ones reported in Ref. (46) concerning and as well as the single ion anisotropy constants. Note that here is very close to 1 (). Knowing that PBC favors the commensurate configuration when is close to 1, large enough sizes are required to obtain an incommensurate magnetic ground state (GS).
However, a MC simulation with 90902 unit cells was not able to reproduce an incommensurate GS with this set of interactions (). Thus larger sizes of the simulation box were required which are not accessible within reasonable computer time (47). Therefore we enhanced the lattice distortion by a factor of 30 (i.e. ). We found that the new set of ’s () and anisotropy constants (Table 1) is a good candidate to reproduce an incommensurate GS for a system of reasonable size 45452 unit cells. It is worth noting that the considered distortion mainly affect the first nearest neighbors interactions while the remaining interactions are not affected. Also it is very interesting to note that the magnitude of the in-plane anisotropy constant () increases when enhancing the lattice distortion reflecting that this anisotropy results from the lattice distortion.
iv.1 Magnetic properties
In order to characterize the GS configuration and to estimate the Néel temperature we performed a first set of simulations without applying an external magnetic field. The following procedure has been retained: we started the simulations from random spin configurations at a high enough temperature () and we then cooled down to = 0.01 K with a constant temperature step = 0.5 K.
In order to estimate the Néel temperature, we calculated the specific heat per spin defined as
where with being the energy of each magnetic configuration, means thermal average, is the number of spins and is the Boltzmann constant. For the parameter set given in Table 1 () the phase transition as signaled by the peak of the specific heat (Fig. 3) takes place at K. This value is in a good agreement with the reported experimental values ( K) (9); (30); (28). This may be taken as a first validation of the extracted exchange interactions of Table 1.
To characterize the nearly GS configuration we considered the spin chirality defined as
where 1, 2 and 3 refer to the spins at the corners of each elementary triangular plaquette in an plane. Then we defined the order parameter per plane to be where is the number of plaquettes per plane, and finally the order parameter of the whole system was defined as where is the average of over the planes. We found that the direction of the vector chirality () of each plane is pointing along the  direction confirming the fact that the spins are oriented in the (110) plane as reported in Ref. (20). Fig. 4 shows the variation of the order parameter as function of temperature. At K, indicates a small deviation from the commensurate (120) configuration of . Moreover, the simulated value of q (0.322, 0.322, 0) confirms that the GS is an incommensurate configuration very close to the reported experimental configuration of q (0.329, 0.329, 0) (28). This good agreement may be taken as another validation of the parameters of Table 1.
On the other hand, the magnetic field dependence of the magnetization calculated along the easy axis (-axis) shows a linear behavior () confirming the antiferromagnetic nature of the GS (not shown here).
Magnetic properties under 0.5 T were simulated between 300 K and 2 K to estimate the CurieWeiss temperature (). Fig. 5 shows the variation of the magnetization and inverse susceptibility measured along the applied magnetic field. It can be seen that obeys well the CurieWeiss law for antiferromagnets (, with is the Curie constant) at high temperatures with K close to the measured experimental values ( K) (9); (48). The curve starts to deviate from the linear behavior at about 100 K. In order to understand the origin of this deviation we calculated the temperature dependence of the spin-spin correlation function defined as along the  direction. As shown in Fig. 6, short-range antiferromagnetic correlations start to develop below 100 K, which leads to the deviation from the CurieWeiss law seen in Fig. 5. Furthermore, these correlation functions exhibit inflection points close to estimated from the specific heat curve (Fig. 3). Besides, an anomaly in the magnetization curve (Fig. 5) appears at K consistent with the estimate of from the specific heat curve. We note that the ratio () reflects the frustrated nature of the GS (49); (50).
iv.2 Ferroelectric properties
In this section, we considered the Hamiltonian . In these simulations, we applied a poling electric field during the cooling process to obtain a single ferroelectric domain. We then turned it off just before statistical averaging to calculate which is associated to the spontaneous ferroelectric polarization (Eq. (2)) according to Ref. (35). Fig. 7 shows the temperature dependence of , the projection of along the  direction, which starts to develop at . It is clearly seen that by switching the poling electric field, can be reversed.
Further insight into the degree of electrical polarization may be gained through the knowledge of the - hysteresis loops, which are shown in Fig. 8 at different temperatures. shows a linear dependence without hysteresis above because the system is in the paraelectric phase, while clear hysteresis loops are seen for temperatures below . This strongly suggests that ferroelectricity is induced by the out-of-plane incommensurate magnetic configuration, in agreement with Ref. (31).
Also, it can be seen that below the saturation field MV/m is independent of the temperature. The hysteresis loop simulated at K shows an electric coercive field for reversal MV/m very close to that measured experimentally ( MV/m (51)). Note that the reversal of results from the reversal of the helicity of each atomic plane. Thus our simulations confirm the electric control of spin helicity in CuCrO as reported in Ref. (20).
In this paper, we proposed estimates of the exchange interactions and single ion anisotropy constants in the multiferroic CuCrO using DFT calculations. They were checked against the experimental Néel and CurieWeiss temperatures as well as the electric coercive field, thereby proving them to be good candidates to model the magnetoelectric properties of CuCrO. We showed that the lattice distortion that takes place below is responsible for the appearance of a weak in-plane hard-axis anisotropy. Regarding the magnetic properties, we obtained a peak in the specific heat curve at K very close to the experimental observations. Furthermore the ground-state has been shown to be an antiferromagnetic incommensurate proper-screw configuration. The estimated K is in a good agreement with experimental data too. Also, our simulated - hysteresis loops confirm the electric control of spin helicity which is related to the ferroelectric polarization below .
We gratefully thank M. Alouani and S. Hébert for stimulating discussions. We are grateful to the Centre Régional Informatique et d’Applications Numeriques de Normandie (CRIANN) where our simulations were performed as project number 2015004. We also acknowledge the computational resources provided by the Swedish National Infrastructure for Computing (SNIC) and Uppsala Multidisciplinary Center for Advanced Computational Science (UPPMAX). The authors acknowledge the financial support of the French Agence Nationale de la Recherche (ANR), through the program Investissements d’Avenir (ANR-10-LABX-09-01) and LabEx EMC3.
- C. Friedel, Sciences Academy 77, 211 (1873).
- R. D. Shannon, D. B. Rogers and C. T. Prewitt, Inorg. Chem. 10, 713 (1971); C. T. Prewitt, R. D. Shannon and D. B. Rogers, Inorg. Chem. 10, 719 (1971); D. B. Rogers, R. D. Shannon, C. T. Prewitt and J. L. Gillson, Inorg. Chem. 10, 723 (1971).
- H. Kawazoe, M. Yasukawa, H. Hyodo, M. Kurita, H. Yanagi, and H. Hosono, Nature 389, 939 (1997).
- H. Takatsu, S. Y. Onezawa, S. M. Ouri, S. Nakatsuji, K. T. Anaka, and Y. Maeno, J. Phys. Soc. Jpn 76, 104701 (2007).
- C. W. Hicks, A. S. Gibbs, A. P. Mackenzie, H. Takatsu, Y. Maeno, and E. A. Yelland, Phys. Rev. Lett. 109, 116401 (2012).
- C. W. Hicks, A. S. Gibbs, L. Zhao, P. Kushwaha, H. Borrmann, A. P. Mackenzie, H. Takatsu, S. Yonezawa, Y. Maeno, and E. A. Yelland, Phys. Rev. B 92, 014425 (2015).
- P. Kushwaha, V. Sunko, P. J. W. Moll, L. Bawden, J. M. Riley, N. Nandi, H. Rosner, M. P. Schmidt, F. Arnold, E. Hassinger, T. K. Kim, M. Hoesch, A. P. Mackenzie, and P. D. C. King, Science Advances 1, 1500692 (2015).
- R. Daou, R. Frésard, S. Hébert, and A. Maignan, Phys. Rev. B 91, 041113(R) (2015).
- T. Okuda, N. Jufuku, S. Hidaka, and N. Terada, Phys. Rev. B 72, 144403 (2005).
- T. Nozaki, K. Hayashi, and T. Kajitani, J. Chem. Eng. Japn 40, 1205 (2007).
- K. Kuriyama, M. Nohara, T. Sasagawa, K. Tabuko, F. Mizokawa, K. Kimura, and H. Takagi, Proc. 25th Int. Conf. Thermoelectrics (IEEE, Piscataway, 2006), p. 97.
- A. Maignan, V. Eyert, C. Martin, S. Kremer, R. Frésard, and D. Pelloquin, Phys. Rev. B 80, 115103 (2009).
- S. Kremer and R. Frésard, Ann. Phys. (Berlin) 524, 21 (2012).
- M. Mekata, N. Yaguchi, T. Takagi, S. Mitsuda, and H. Yoshizawa, J. Magn. Magn. Mater. 823, 104 (1992).
- M. Mekata, N. Yaguchi, T. Takagi, T. Sugino, S. Mitsuda, H. Yoshizawa, N. Hosoito, and T. Shinjo, J. Phys. Soc. Japan 62, 4474 (1993).
- O. A. Petrenko, G. Balakrishnan, M. R. Lees, D. McK. Paul, and A. Hoser, Phys. Rev. B 62, 8983 (2000).
- T. Kimura, J. C. Lashley, and A. P. Ramirez, Phys. Rev. B 73, 220401(R) (2006).
- F. Ye, Y. Ren, Q. Huang, J. A. Fernandez-Baca, P. Dai, J. W. Lynn, and T. Kimura, Phys. Rev. B 73, 220404(R) (2006).
- V. Eyert, R. Frésard, and A. Maignan, Phys. Rev. B 78, 052402 (2008).
- M. Soda, K. Kimura, T. Kimura, M. Matsuura, and K. Hirotam, J. Phys. Soc. Jpn. 78, 124703 (2009).
- A. Maignan, C. Martin, R. Frésard, V. Eyert, E. Guilmeau, S. Hébert, M. Poienar, and D. Pelloquin, Solid Stat. Comm. 149, 962 (2009).
- J. Xue-Fan, L. Xian-Feng, W. Yin-Zhong, and H. Jiu-Rong, Chin. Phys. B 21, 077502 (2012).
- I. I. Mazin, Phys. Rev. B 75, 094407 (2007).
- J. T. Haraldsen, F. Ye, R. S. Fishman, J. A. Fernandez-Baca, Y. Yamaguchi, K. Kimura, and T. Kimura, Phys. Rev. B 82, 020404(R) (2010).
- H. Kadowaki, H. Kikuchi, and Y. Ajiro, J. Phys.: Condens. Matter 2, 4485 (1990).
- M. Soda, K. Kimura, T. Kimura, and K. Hirota, Phys. Rev. B 81, 100406(R) (2010).
- M. Frontzek, G. Ehlers, A. Podlesnyak, H. Cao, M. Matsuda, O. Zaharko, N. Aliouane, S. Barilo, and S. V. Shiryaev, J. Phys.: Condens. Matter 24, 016004 (2012).
- M. Poienar, F. Damay, C. Martin, V. Hardy, A. Maignan, and G. André, Phys. Rev. B 79, 014412 (2009).
- K. Kimura, T. Otani, H. Nakamura, Y. Wakabayashi, and T. Kimura, J. Phys. Soc. Jpn. 78, 113710 (2009).
- S. Seki, Y. Onose, and Y. Tokura, Phys. Rev. Lett. 101, 067204 (2008).
- K. Kimura, H. Nakamura, K. Ohgushi, and T. Kimura, Phys. Rev. B 78, 140401(R) (2008).
- Y. Tokura and S. Seki, Adv. Mater. 22, 1554 (2010).
- Y. Tokura, S. Seki, and N. Nagaosa, Rep. Prog. Phys. 77, 076501 (2014).
- N. Terada, J. Phys.: Condens. Matter 26, 453202 (2014).
- T. A. Kaplan and S. D. Mahanti, Phys. Rev. B 83, 174432 (2011).
- J. M. Wills, O. Eriksson, M. Alouani, and D. L. Price, in Electronic Structure and Physical Properties of Solids, Lecture Notes in Physics, Vol. 535, edited by H. Dreysse (Springer Berlin Heidelberg, 2000) p. 148-167.
- Y. Ono, K. I. Satoh, T. Nozaki, and T. Kajitani, Jpn. J. Appl. Phys. 46, 1071 (2007).
- A. I. Liechtenstein, V. I. Anisimov, and J. Zaanen, Phys. Rev. B 52, R5467(R) (1995).
- M. T. Czyzyk and G. A. Sawatzky, Phys. Rev. B 49, 14211 (1994).
- A. I. Liechtenstein, M. I. Katsnelson, V.P. Antropov, and V. A. Gubanov, J. Magn. Magn. Mater. 67, 65 (1987).
- M. I. Katsnelson and A.I. Lichtenstein, Phys. Rev. B 61, 8906 (2000).
- Y. O. Kvashnin, O. Grånäs, I. Di Marco, M. I. Katsnelson, A. I. Lichtenstein, and O. Eriksson, Phys. Rev. B 91, 125133 (2015).
- D.P. Landau and K. Binder, A Guide to Monte Carlo Simulations in Statistical Physics (Cambridge University Press, Cambridge, England, 2008).
- N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, J. Chem. Phys. 21, 1087 (1953).
- U. Nowak, R. W. Chantrell, and E. C. Kennedy, Phys. Rev. Lett. 84, 163 (2000).
- H. Yamaguchi, S. Ohtomo, S. Kimura, M. Hagiwara, K. Kimura, T. Kimura, T. Okuda, and K. Kindo, Phys. Rev. B 81, 033104 (2010).
- Note that systems larger than 90902 unit cells require more than 12.5 days of simulation which is not accessible at the super-computer of CRIANN.
- T. Okuda, R. Kajimoto, M. Okawa, and T. Saitoh, Int. J. Mod. Phys. B 27, 1330002 (2013).
- A. P. Ramirez, Annu. Rev. Mater. Sci. 24, 453 (1994).
- J. E. Greedan, J. Mater. Chem. 11, 37 (2001).
- K. Kimura, H. Nakamura, S. Kimura, M. Hagiwara, and T. Kimura, Phys. Rev. Lett. 103, 107201 (2009).