Monte Carlo simulation on rotation of ferroelectric polarization by rotating magnetic field in conical-spin–ordered multiferroics

Xiaoyan Yao¹ and Qichang Li²

¹ Department of Physics, Southeast University - Nanjing 211189, China
² School of Physical Science, Qingdao University - Qingdao 266071, China

E-mail: yaoxiaoyan@gmail.com

Received 16 August 2009, accepted for publication 26 October 2009
Published 24 November 2009

Introduction

Multiferroics, especially magnetoelectric (ME) materials, in which two or more ferroic orders coexist and mutually interact, have attracted considerable attention due to their potential applications and fundamental curiosity. In recent years, the multiferroicity of magnetic origin observed in the frustrated materials renews the interest in this field. Here the ferroelectricity results from the loss of inversion symmetry induced by the spatial variation of magnetization (M) [1, 2]. Recent experiments indicate that both collinear and noncollinear spin orders can lead to ferroelectricity [3, 4]. Especially, the multiferroicity originated from noncollinear spin order has been observed in a lot of compounds, such as RMnO₃ (R = Tb and Dy) [5], MnWO₄ [6] and LiCu₂O₂ [7]. Two possible microscopic mechanisms have been proposed, that is, the spin current model and the theory of inverse Dzyaloshinskii-Moriya interaction [8–10]. It is interesting that both mechanisms give the same relation between ferroelectric polarization (P) and the neighboring canting spins (S_i and S_j), namely,

where e_ij denotes the vector connecting two sites of S_i and S_j, and a is a factor determined by the spin exchange coupling and the spin-orbit interaction.

Due to the magnetic origin of ferroelectricity, the ferroelectric behavior in these materials is highly sensitive to the applied magnetic field (h). Various magnetic controls of ferroelectric behavior have been observed, such as flop, reversal, rotation of P and the anomaly of the dielectric constant [2, 11, 12]. Among them, the rotation of P controlled by h is one of the noticeable phenomena. A lot of interest is attracted due to its advantage for designing novel electric devices, such as sensitive magnetic field sensors. Different from the switch of P between two directions perpendicular to each other observed in TbMnO₃ [13], the direction of P can be continuously controlled by h in some other compounds, such as Ba₂Mg₂Fe₁₂O₂₂ [14], ZnCr₂Se₄ [15] and Eu_0.55Y_0.45MnO₃ [16]. It is argued that the conical spin order is essential to realize this continuous rotation. However, in these compounds the conical spin structure is generally obtained by applying a weak h to screw or cycloidal spiral structures. Different from these h-induced conical spin orders, the conical spin structure observed in the spinel compound CoCr₂O₄ is inherent. Experiments revealed that P can be reversed during the sweep of h (without changing its direction) in this compound [17–19]. It is indicated in our earlier work that the flexible conical spin order is responsible for this ME phenomenon [20]. Therefore, the inherent conical spin order in spinel structure should be a better candidate to realize a h-controlled rotation of P. To our knowledge, the work on the ferroelectric behavior upon application of rotating h in spinel structure has rarely been reported, and this ME behavior together with its microscopic mechanism remains an interesting question.

In this paper, based on the preliminary work [20], a three-dimensional (3D) model with classical Heisenberg spins in spinel structure is employed to investigate the magnetic and ferroelectric responses to rotating h by using a Monte Carlo simulation. It is observed that M always sticks to the direction of h and rotates with h together without decay, while P shows an anisotropic behavior for the different paths in which h rotates. Especially, when h is kept perpendicular to the magnetic modulation vector k, P demonstrates a continuous rotation without any noticeable decay. A detailed discussion on the microscopic structure of spins confirms that these ME behaviors are ascribed to the flexible rotation of conical spin structure. It is revealed that as long as the conical spin structure is formed in the ME poling process, it can be well retained and smoothly rotated during the whole procedure, in whichever path h rotates. Thus, the macroscopic magnetic and ferroelectric properties can be flexibly tuned by h through handling the conical spin structure.

Model and simulation

The cubic spinel structure has the general formula AB₂O₄, where A and B occupy the tetrahedral and octahedral sites, respectively. Previous investigations indicate that only the nearest-neighbor antiferromagnetic A-B and B-B exchange interactions (J_AB and J_BB) are dominant in CoCr₂O₄, and sufficient to produce a conical spin state [21–23]. Therefore, considering J_AB and J_BB, a classical Heisenberg exchange energy is given by

where S_Bi and S_Ai are i-th normalized classic spins on B and A sites with unit magnitude. [i, j] denotes the summation over all the nearest-neighboring spin pairs. It has been shown in the theory of Lyons, Kaplan, Dwight, and Menyuk (LKDM) that for a cubic spinel structure the magnetic ground state is determined by the parameter [24, 25]

When 8/9 < u < 1.298, it is suggested that the conical spin order has the lowest energy out of a large set of possible spin configurations. In other words,

Therefore J_BB = −2.5 and J_AB = −3 are supposed in the present simulation. Considering the magnetic energy and the electric energy, the total Hamiltonian can be written as

Here E is the external electric field and M is evaluated as

where N is the total number of magnetic ions. P is calculated according to eq. (1). For convenience, only the spin current between spins of the B site is considered and the factor a is assumed to be unity. Therefore P is written as

where e_ij is in the direction of magnetic modulation vector k.

The Monte Carlo simulation is performed on a 3D spinel lattice composed of classical Heisenberg spins. The system size is L × L × L, where L = 36, and the periodic boundary condition is applied. It is assumed that x, y and z axes are, respectively, along the direction of [100], [010] and [001] in the spinel structure. The spin configuration is updated according to the Metropolis algorithm and a trial step is obtained by randomly changing the direction of the spin from its original direction a little. Only the behaviors at low temperature T = 0.01 are focused. The system is initially polarized by E and h. After that, E is removed, and then M and P are calculated with different h directions. For every h direction, the initial 10000 Monte Carlo steps (MCSs) are discarded for equilibration, and then the results are obtained by averaging 1000 data. Each datum is collected every 10 MCSs. The final results are obtained by averaging independent data sets calculated by selecting different seeds for random number generation.

Results and discussion

Similar to the ME poling procedure applied before measurement in the experiments [17, 19], in the simulation the system is initially polarized by a poling magnetic field h = 1.2 along [001] and a poling electric field E = 0.7071 along [110] according to the anisotropy of the spinel structure. This ME poling procedure produces a single ME domain, and fixes the directions of M and P to [001] and [110], respectively [20]. Microscopically, the conical spin structure appears in the spin chains of the B site and the magnetic modulation vector k is aligned along , as shown in fig. 1(a). If all the ionic sites in the chain are moved to one site, then the spin vectors will lie on a surface of cone, namely, the spin cone, as presented in fig. 1(b). The arrow on the edge of the cone shows that the spins in the chain rotate counterclockwise along k. As shown in fig. 1(c), the homogeneous components of the spins along the cone axis produce M along [001] and their spiral parts give rise to P along [110].

After the ME poling procedure, E is turned off. And then the rotation of h, with its magnitude hL = 1.5 kept unchanged, is started from the polarized direction [001]. Here θ is defined as the angle between the z-axis and the orientation of h, as given in fig. 2. It is worthwhile to note that due to the inherent anisotropy in the spinel structure, the ferroelectric behavior should be different when h varies in different paths. Since the direction of k is a key to decide the emergence of P according to eq. (7), three typical processes with h rotating in different planes relative to k are discussed as follows. Here the plane in which h rotates is named R-plane. The R-plane is always normal to the xy-plane because h always rotates from [001].

Case I: R-plane normal to k.

As shown in figs. 2(a) and (b), starting from [001], h is slanted towards [110] and rotated by 2π, lying in the plane normal to k. The θ-dependences of M and P are plotted in figs. 2(c) and (d). It is observed that M follows the variation of h closely and sticks to the direction of h all along. As h and M vary together, the components of P also demonstrate well-defined wavelike curves, which imply that P rotates with h simultaneously and smoothly, but its orientation is always perpendicular to h. Moreover, the calculation indicates that both the magnitudes of M and P (ML and PL) remain unchanged in the whole course of rotating h. Therefore, during the rotation of h with its direction always perpendicular to k, both M and P show continuous rotation without decay in the plane normal to k. At the same time, k, M and P are kept perpendicular to each other all along.

Case II: The R-plane is tilted by π/4 from k.

h is slanted from [001] towards [100] (x-axis) and rotated by 2π in the zx-plane, as shown in figs. 2(e) and (f). In this case, M also rotates following h closely with ML unchanged (see fig. 2(g)). The same variations of the x- and y-components of P ensure that P still rotates smoothly in the plane perpendicular to k, but its length varies wavelike, as presented in fig. 2(h). When θ = π/2 or 3π/2, namely both h and k are within the xy-plane, P shows the minimal length. When θ = 0 or π, namely h is perpendicular to k, PL reaches its maximum. Therefore, when h rotates in the plane tilted from k by an acute angle, P rotates in the plane perpendicular to k with its length varying periodically, while M still rotates continuously together with h without decay.

Case III: The R-plane is coplanar with k.

As shown in figs. 2(i) and (j), h is slanted from [001] towards and rotates by 2π within the same plane with k. It is plotted in figs. 2(k) and (l) that M also rotates without decay and keeps to the direction of h, but P dose not rotate anymore. The same variations of the x- and y-components of P together with the zero value of the z-component indicate that P is restricted to the line along [110] or . When θ = π/2 or 3π/2, namely h is parallel to k, P disappears. When θ = 0 or π, namely h is perpendicular to k, PL reaches its maximum. Therefore, during the rotation of h with its direction kept coplanar with k, M still rotates continuously together with h without decay, but P just stretches and shrinks in the line perpendicular to the R-plane.

In order to explore this anisotropic ME behavior upon rotating h, the microscopic spin structure is investigated in detail. To characterize the change of spin structure, four parameters listed below are monitored during the whole procedure of h rotation. As illustrated in fig. 3(a), the radius of the spin cone (r) and the midpoint of bottom C(Cx, Cy, Cz) are used to scrutinize the shape and orientation of the spin cone. In addition, the averaged cross product S_i × S_j (CP) and dot product S_i · S_j (DP) between the nearest-neighboring spins along k are calculated to observe the rotation direction of spins and the average angle between the nearest-neighboring spins. The θ-dependences of r, C, CP and DP are demonstrated in figs. 3(c)–(h). It is seen that r keeps almost invariant for the three cases, which means that the spin cone is well reserved as a whole in the course of h rotating. The three coordinates of C show very similar track curves to those of h. This indicates that the axis of the spin cone rotates smoothly and always coincides with the direction of h. Moreover, the three components of CP also present very similar curves, which means that the orientation of the spin rotation axis also coincides with h. In addition, the almost invariable DP implies that the modulation period of the spin structure hardly changes as h rotates. Therefore, the conical spin structure is always well retained and flexibly rotated by h during the whole procedure, in whichever direction h rotates. On the other hand, the strength of h also influences the ME behavior. As illustrated in fig. 3(b), when hL increases in Case I, M is enhanced and P is suppressed, which is attributed to the microscopic change of spin structure, namely, the spins tend to the direction of h to reduce the Zeeman energy, so the spin cone becomes slimmer.

Based on the aforementioned character of the spin structure, the magnetic and ferroelectric behaviors can be well understood. According to eq. (7), P should be perpendicular to k and CP. Since CP always keeps to the direction of h, P should be normal to k and h simultaneously. In addition, owing to the well-reserved conical spin structure, the magnitude of CP remains unchanged. Therefore, PL should be proportional to |sin |, where is the angle between k and h. If the angle between the R-plane and k is denoted by α, as shown in fig. 4(a), then

Here only the α-range from 0 to π/2 (i.e. from the direction of k to [110]) needs consideration. Other areas are similar to this α-range. |sin |'s evaluated according to eq. (8) for different α and θ are plotted in fig. 4(b). It is seen that the simulated PL curves scaled by their maximum (PL_max) are well consistent with the corresponding curves of |sin |, which confirms that the normalized magnitude of the P vector can be represented by |sin |. When α = π/2, namely the R-plane is normal to k, PL shows a constant behavior with its maximal value. When 0 < α < π/2, namely k is tilted from the R-plane by an angle, all the curves reach maximum at θ = 0 and π where h is normal to k. The minimum emerges at θ = π/2 and 3π/2 where the angle between h and k is minimal. As α is reduced, the minimal angle between h and k decreases, and thus the minimum of PL also decreases. When α reaches 0, the minimum of PL also reaches 0.

Figure 4. (Color online) (a) The sketch of the xy-plane where the angle between the R-plane and k is denoted by α. Here only the α-range denoted by the curved arrow (i.e. from k to [110]) needs consideration. (b) |sin | calculated according to eq. (8) as a function of θ for different α is plotted by symbols, and the simulated PL curves scaled by PLmax are plotted by dotted lines. (c) Polar plot of the P vector with its magnitude and direction for different α and θ. The symbol represents the position that the top of the P vector reaches and the normalized magnitude of P is given by |sin |.

According to the above discussion, both direction and magnitude of P can be decided by the direction of h relative to k. Since P only appears within one plane, i.e. the plane normal to k, the P vector with its magnitude and direction at different θ and α can be plotted in a polar map as fig. 4(c), where the symbol represents the position that the top of the P vector reaches and the length of the P vector is given by |sin |. It is seen that the trace of the P vector is a round circle for α = π/2. When α decreases, the two sides of the trace circle along the z-axis shrink while the other two sides keep unchanged. As α = 0, the whole circle shrinks to a line parallel to [110].

The simulation indicates that the conical spin order exists stably in the spinel structure. In the conical spin state, the homogeneous components of spins along the cone axis produce M while their spiral parts give rise to P, which induces an inherent M-P coupling. When h rotates, M always sticks to the direction of h to lower the Zeeman energy, and then P is rotated by h through this inherent M-P coupling. Although the simulation is performed on the spinel structure, our discussion focuses on the cases with different h rotating paths relative to k. Therefore the obtained results provide a wide-use principle for h-controlled rotation of ME properties, which can be generally applied to the conical spin structure in other materials. For example, it is observed that P can be rotated smoothly by rotating h around k in Eu_0.55Y_0.45MnO₃ [16], which is attributed to the flexible rotation of the h-induced conical spin structure, consistent with our simulated results.

Conclusion

In summary, aiming at the multiferroicity originated from conical spin order, we investigated the variation of the magnetic and ferroelectric behaviors under rotating h by a Monte Carlo simulation based on a classical Heisenberg model in spinel structure. Considering the anisotropy of crystal structure, different positions of the R-plane relative to k are studied. When k is out of the R-plane, a flexible rotation of P controlled by h is observed in the plane normal to k. Especially when k is perpendicular to the R-plane, P rotates by 2π continuously without any noticeable decay. The analysis of the spin structure indicates that the conical spin order, retained very well during the process of h rotating, is responsible for the flexible rotation of P. In order to lower the Zeeman energy, M sticks to the direction of h all along, which makes the axis of the spin cone rotate simultaneously, and hence P is rotated. These results provide more insights into the multiferroicity driven by conical spin order.

Acknowledgments

This work is supported by the research grants from the National Natural Science Foundation of China (10904014), and supported by the computational center from Department of Physics, Southeast University.

References

[1]

Cheong S.-W. and Mostovoy M. 2007 Nat. Mater. 6 13
CrossRefPubMed

[2]

Wang K. F., Liu J.-M. and Ren Z. F. 2009 Adv. Phys. 58 321
CrossRef

[3]

Choi Y. J., Yi H. T., Lee S., Huang Q., Kiryukhin V. and Cheong S.-W. 2008 Phys. Rev. Lett. 100 047601
CrossRefPubMed

[4]

Yao X. and Lo V. C. 2008 J. Appl. Phys. 104 083919
CrossRef
Yao X., Lo V. C. and Liu J.-M. 2009 J. Appl. Phys. 105 033907
CrossRef

[5]

Kimura T., Goto T., Shintani H., Ishizaka K., Arima T. and Tokura Y. 2003 Nature 426 55
CrossRefPubMed

[6]

Taniguchi K., Abe N., Takenobu T., Iwasa Y. and Arima T. 2006 Phys. Rev. Lett. 97 097203
CrossRefPubMed

[7]

Park S., Choi Y. J., Zhang C. L. and Cheong S.-W. 2007 Phys. Rev. Lett. 98 057601
CrossRefPubMed

[8]

Katsura H., Nagaosa N. and Balatsky A. V. 2005 Phys. Rev. Lett. 95 057205
CrossRefPubMed

[9]

Sergienko I. A. and Dagotto E. 2006 Phys. Rev. B 73 094434
CrossRef

[10]

Li Q. C., Dong S. and Liu J.-M. 2008 Phys. Rev. B 77 054442
CrossRef

[11]

Hur N., Park S., Sharma P. A., Ahn J. S., Guha S. and Cheong S.-W. 2004 Nature (London) 429 392
CrossRefPubMed

[12]

Strempfer J., Bohnenbuck B., Mostovoy M., Aliouane N., Argyriou D. N., Schrettle F., Hemberger J., Krimmel A. and Zimmermann M. V. 2007 Phys. Rev. B 75 212402
CrossRef
Aliouane N., Argyriou D. N., Strempfer J., Zegkinoglou I., Landsgesell S. and Zimmermann M. V. 2006 Phys. Rev. B 73 020102
CrossRef

[13]

Abe N., Taniguchi K., Ohtani S., Takenobu T., Iwasa Y. and Arima T. 2007 Phys. Rev. Lett. 99 227206
CrossRefPubMed

[14]

Ishiwata S., Taguchi Y., Murakawa H., Onose Y. and Tokura Y. 2008 Science 319 1643
CrossRefPubMed

[15]

Murakawa H., Onose Y., Ohgushi K., Ishiwata S. and Tokura Y. 2008 J. Phys. Soc. Jpn. 77 043709
CrossRef

[16]

Murakawa H., Onose Y., Kagawa F., Ishiwata S., Kaneko Y. and Tokura Y. 2008 Phys. Rev. Lett. 101 197207
CrossRefPubMed

[17]

Choi Y. J., Okamoto J., Huang D. J., Chao K. S., Lin H. J., Chen C. T., van Veenendaal M., Kaplan T. A. and Cheong S.-W. 2009 Phys. Rev. Lett. 102 067601
CrossRefPubMed

[18]

Kim I., Oh Y. S., Liu Y., Chun S. H., Lee J.-S., Ko K.-T., Park J.-H., Chung J.-H. and Kim K. H. 2009 Appl. Phys. Lett. 94 042505
CrossRef

[19]

Yamasaki Y., Miyasaka S., Kaneko Y., He J.-P., Arima T. and Tokura Y. 2006 Phys. Rev. Lett. 96 207204
CrossRefPubMed

[20]

Yao X., Lo V. C. and Liu J. M. 2009 J. Appl. Phys. 106 073901
CrossRef

[21]

Lawes G., Melot B., Page K., Ederer C., Hayward M. A., Proffen T. and Seshadri R. 2006 Phys. Rev. B 74 024413
CrossRef

[22]

Ederer C. and Komelj M. 2007 Phys. Rev. B 76 064409
CrossRef

[23]

Kaplan T. A., Dwight K., Lyons D. and Menyuk N. 1961 J. Appl. Phys. 32 S13
CrossRef

[24]

Kaplan T. A. 1960 Phys. Rev. 119 1460
CrossRef

[25]

Lyons D. H., Kaplan T. A., Dwight K. and Menyuk N. 1962 Phys. Rev. 126 540
CrossRef