Monte Carlo simulation

Multiferroic response and clamped domain structure in a two-dimensional spiral magnet:

Monte Carlo simulation

Qichang Li, Shuai Dong, and J.-M. Liu

*

Laboratory of Solid State Microstructure, Nanjing University, Nanjing 210093, China

and International Center for Materials Physics, Chinese Academy of Sciences, Shenyang 110015, China 共Received 5 November 2007; revised manuscript received 7 January 2008; published 29 February 2008兲 The Monte Carlo simulation on multiferroic behaviors of the two-dimensional MnO₂lattice in multiferroic manganites is performed based on the model关I. A. Sergienko and E. Dagotto, Phys. Rev. B 73, 094434共2006兲兴 associated with the Dzyaloshinskii-Moriya interaction. The simulated ferroelectric polarization induced by the spiral spin ordering and its response to external magnetic field agree well with reported experimental obser- vations. Furthermore, the coexistence of clamped ferroelectric domains and spiral spin domains is revealed in our simulation.

DOI:10.1103/PhysRevB.77.054442 PACS number共s兲: 75.80.⫹q, 77.80.Bh, 75.30.Kz, 77.80.Dj

Research interest in multiferroics with coexisting ferro- electric 共FE兲 and ferromagnetic 共FM兲 or antiferromagnetic 共AFM兲 orderings has been revived recently due to their promising technological potentials in storage devices and sensors.^1–3In the past few years, the gigantic magnetoelectric 共ME兲effect, that is, FE polarization Pof the order of mag- nitude of⬃1.0␮^C/cm²or less, was observed in perovskite manganites such as RMnO₃ 共R= Tb, Dy, and Gd兲,^4–6 RMn₂O₅ 共R= Y , Tb, Dy, Ho, Er, etc.兲,^7–9 and kagomé-lattice compound Ni₃V₂O₈.^10,11Huge efforts from both experimen- talists and theorists have been made to understand the new mechanism underlying this effect. At least for some multifer- roic manganites such as TbMnO₃, a spiral-spin configuration might play an important role in generating the FE polarization,^6,12and a change in the magnetic structure at the FE Curie pointTc, wherePis initiated, was demonstrated by neutron diffraction.^13,14 The magnetic symmetry consider- ation also suggested a transition from sinusoidal共collinear兲 spin order to spiral 共noncollinear兲 spin order near Tc.¹⁴ On the other hand, the microscopic origin for this FE order, as proposed by Katsuraet al.,¹⁵is related to the Mn-O-Mn bond structure, in which an overlap of the wave functions between the two adjacent 3d Mn ions with canted-spin arrangement favors a small polarization due to the displacement of the sandwiched oxygen ion. Subsequently, it was suggested that transverse-spiral spin modulation along a specific direction could induce a macroscopic FE polarization.^16,17

The possible form of exchange interactions responsible for the spiral-spin-induced ferroelectricity was proposed in several models. Sergienko and Dagotto addressed the signifi- cance of the antisymmetric Dzyaloshinskii-Moriya interac- tion 共DMI兲 in the appearance of the FE ordering for rare- earth perovskiteRMnO₃. In this model, a ME coupling term D·共S_i⫻S_j兲, where Dis a vector correlating to the displace- ment of the oxygen ion between moments S_i and S_j, is included.¹⁸Harriset al.proposed a more complicated theory to explain the ME effect in Ni₃V₂O₈, in which an exchange tensor is introduced to include the contributions from the superexchange, double exchange, and DMI between the nearest neighboring moments.¹⁹Mostovoy developed a phe- nomenological theory for inhomogeneous multiferroics, where the ME effect is associated with the spatial derivatives of magnetic moment.²⁰

Although preliminary calculations based on the above models highlighted the major features associated with the ME response, more detailed investigations on the interplay between the FE and magnetic orders, in the quantitative sense, would be helpful to illustrate this fantastic effect. Al- though nonzero P under zero magnetic field, induced at an incommensurate magnetic共ICM兲transition, was revealed by Monte Carlo共MC兲simulation,¹⁸yet the detailed dependence ofPon magnetic fieldHhas not been reported. In fact, as far as we know, few works on the MC simulation of such effect have been reported.

To perform a MC simulation, one needs to understand the interactions between order parameters. Raman study indi- cates that the in-plane oxygen vibrations are mainly involved in the spin-phonon coupling in orthorhombic RMnO₃ 共R

= Tb and Dy兲.²¹Although theabplanes of these materials are antiferromagnetically stacked along thecaxis,¹⁴the direction of inducedPis identical in theseabplanes. It is, therefore, a worthwhile attempt to study the ME effect of a two- dimensional 共2D兲 system in which the oxygen movements are limited within a plane. The spins coupled with the O displacements are treated as in-plane too. Surely, it does not mean that the spins of Mn ion inRMnO₃are limited within one specific plane. However, one can project the spins onto a plane, which will not lose the physical essence.¹⁴ We start from the DMI model developed by Sergienko and Dagotto¹⁸ and perform a detailed calculation of the dependence ofPon H. The calculated results will be compared with recent ex- periments. We further demonstrate how the FE domains and magnetic domains are interactively clamped, illustrating the microscopic pictures of the ME effect in multiferroics.

The MC simulation is performed on a 2D L⫻L square lattice with periodic boundary conditions, representing the MnO₂plane of multiferroic RMnO₃. The plane is schemati- cally illustrated in Fig.1共a兲, where each Mn ion共blue dot兲is surrounded by four O ions共red dots兲. The magnetization M is assumed from the contribution of the Mn spins, andP is attributed to the displacement of O ions. Both of them are restricted within thexy共or ab兲plane. For such a lattice, we directly take the original Hamiltonian¹⁸ based on the orbit- ally degenerate double-exchange model:

PHYSICAL REVIEW B77, 054442共2008兲

1098-0121/2008/77共5兲/054442共4兲 054442-1 ©2008 The American Physical Society

H˜ = −_ia

兺

兺

兺

−␮Bg_L

兺

兺

兺

+␬1

2

兺

兺

兺

whered⁺_i␣␴is the creation operator for the electron carrying spin␴ ^{on sitei, and ta}_␣␤ is the hopping integral. JH is the Hund coupling constant between theegelectrons with spinsi

and three t_2g electrons with a classical spin S_i= 1.4 in moment.^14,18 JAF represents the AFM isotropic superex- change between t_2g spins, H denotes the external magnetic field,␮Bis the Bohr magnon,gLrepresents the Lande factor, Edenotes a small electric field, eis the charge of the elec- tron,rirepresents the O displacement, andD^a共ri兲is a vector concerning the displacement of O ions and scaling the mag- nitude of the DMI.QxiandQyidenote the classical FE pho- non coordinates关shown in Fig. 1共b兲, also referring to Ref.

18兴along the xandyaxes, respectively, associated with the displacement of the four O ions surrounding the Mn site i, while Q_mi represents the remnant phonon coordinates. Fi- nally,␬1and␬2denote the spring constants of the FE modes and the remnant modes, andH_JT is the term from the Jahn- Teller共JT兲distortion.

For theabplane,D^a共ri兲in thexandydirections takes the following forms:¹⁸

D^x共r_i兲=␥共0,0,y_i兲, D^y共r_i兲=␥共0,0,−x_i兲, 共2兲 where␥is a constant scaling the coupling intensity between the spins and the displacements. TermH_JT is redefined to fit the 2D model, based on the three-dimensional form:^22,23

HJT=␭共q1i␳i+q_2i␶xi+q_3i␶zi兲+1

2

兺

共3兲 where␳iand␶iare the orbital pseudospins, andq_1i,q_2i, and q_3irepresent the JT breath modes and the stretch modes. For

simplicity, the contribution of theegelectrons is not consid- ered, referring to highly localized manganites. Hence, the first two terms of Eq. 共1兲 and the first term of Eq. 共3兲 are simply discarded. The exclusion of these terms is physically acceptable since it has been proven that both the double ex- change and superexchange between noncollinear spins are able to induce the FE polarization.¹⁵

For subsequent comparison with experiments, we high- light the choice of the parameters taken for our simulation 共TableI兲. The reported Curie temperatureTcof DyMnO₃and TbMnO₃is of the order of 10 K. The value ofPranges from 0.08 to 0.25␮^C/cm².^4,5Hence the derived magnitude of the FE phonon coordinate Q is of the order of 10⁻³– 10⁻²Å.

Following the routine of Zhou and Goodenough,²⁴ ␬1

=kBTc/Q², which gives␬1/kB⬃10⁵K/ŲwithQ⬃10⁻² Å;

here, kB denotes the Boltzmann constant. ␬2= 10␬1 is taken so that the system favors the ICM-FE phase. Considering the difference of the Néel temperature between LaMnO₃ 共Refs.

25 and26兲 and RMnO₃ 共R= Dy and Tb兲, JAF/kB= 2 K is a reasonable choice.D共ri兲/kB= 15 K is chosen in order to tune the ICM transition point 共⬃12 K in this work兲 to fit the experimental value. Correspondingly, ␥= 1.5⫻10³K/Å in our simulation. It is noticed that D共r_i兲 is about 1 order of magnitude larger than the estimated value for LaMnO₃.^25,26 However, this discrepancy may be diminished when the model is improved, as illustrated in previous work.¹⁸

The details of the simulation procedure are the same as reported earlier.¹⁸ For every temperature T in the P-T and M-T curves, initial 12 500 Monte Carlo steps 共MCSs兲 are discarded for equilibration and 47 500 MCSs retained for average. For every specificH in the P-Hand M-H curves, initial 500 MCSs are discarded and 1500 MCSs retained for average. The Mn spin and O displacement are updated ac- cording to the Metropolis algorithm, with the displacement restricted between −10⫻10⁻² and 10⫻10⁻²Å in the updat- ing process.

According to the previous report,¹⁸ Palong a specific di- rection共baxis兲is induced due to the ICM phase transition in the 2D lattice. Nevertheless, given the fourfold symmetry of the lattice, four possible directions ofP共parallel and antipar- allel toaaxis andbaxis, respectively兲corresponding to four types of spiral-spin configuration are, in fact, allowed with equal possibility. We present here in Figs.2共a兲–2共d兲all four types of ICM-FE configurations which are achieved by ap- plying a smallEparallel to one of the allowed directions. It is somehow equivalent to the experimental poling procedure.

As a result, the fourfold symmetry is broken and only one direction of P is favored over the lattice. Due to the DMI term in Eq.共1兲, the uniquely oriented Phelps to stabilize a specific spiral spin configuration, characterized by a wave FIG. 1.共Color online兲 共a兲Two-dimensional MnO₂lattice for the

MC simulation, and共b兲the FE classical phonon modes associated with the displacements of oxygen ions.

TABLE I. Parameters chosen for the simulation.

Parameter Value Parameter Value

k_B共J/K兲 1.3807⫻10⁻²³ ␥/k_B共K/Å兲 1.5⫻10³

␮B共J/T兲 9.274⫻10⁻²⁴ ␬1共K/Ų兲 1.0⫻10⁵

g_L 2 ␬2共K/Ų兲 1.0⫻10⁶

J_AF/k_B共K兲 2

LI, DONG, AND LIU PHYSICAL REVIEW B77, 054442共2008兲

054442-2

vector k normal to P, which points to the direction along which the noncollinear spins rotate clockwise. For each con- figuration shown in Fig. 2, E 共Ee/kB= 30 K/Å兲 is applied parallel to P. Unless stated otherwise, our simulations are performed with Ee/kB= 30 K/Å applied along the b axis.

This field is⬃10²kV/cm, on the same order of magnitude with experimentally employed value for poling ferroelectric oxide ceramics at lowT.

Figures3共a兲and3共b兲present the simulatedM andPas a function ofT under various H 共H^储a兲, for both cooling and heating sequences. In this paper, M is the average of spins over the lattice andPis equal to the averaged dipole moment divided by the cell volume of DyMnO₃. It is shown that under H= 0, the lattice exhibits almost zero M over the wholeT range, while a clear FE transition occurs at T=Tc

⬃12 K, at which a transition from the paramagnetic state to the ICM phase共as shown in Fig.2兲is identified. The sudden flop ofP and the weak but clear thermal hysteresis indicate the first order FE transition here. For the cases ofH⬎0,Tc

shifts toward the low-Tside with increasingH. Upon further increasing H up to 12 T, the ICM configuration collapses into the spin-canted weak FM state. Since such a spin order does not favor the FE order, we observe the disappearance of P over the wholeT range.

We present in Figs.3共c兲and3共d兲theHdependence ofM andPat several selectedT共Ts兲 belowTc. The simulation is performed as follows. The lattice is first cooled from T

= 40 K down to Ts 共Ts= 3, 5, and 7 K, respectively兲 under H= −18 T, and then H is cycled between H= 18 T and H=

−18 T. The response ofM andPto varyingHis featured by the double-loop hysteresis. The steep steps of the hysteresis in Fig.3共c兲represent, respectively, the destruction共up arrow兲 and construction共down arrow兲 of the spiral-spin configura- tion. In response to the destruction共construction兲,Pbecomes suppressed共resumed兲. WhenTsis higher, the critical fieldHc

for the destruction 共construction兲 becomes lower, and the loop area for theM-HandP-H hysteresis loops shrinks.

If the ICM transition point atH= 0 is adjusted to⬃12 K, comparable to the so-called lock-in temperature of DyMnO₃ 共⬃20 K兲 共Ref. 5兲 and TbMnO₃ 共⬃28 K兲,⁴ the simulated maximal P is ⬃0.3␮C/cm² 共Fig. 3兲, in good agreement with the measured value in DyMnO₃.⁵ The thermal hyster- esis illustrated in Fig.3 was also reported experimentally.^4,5 The simulated minimalHrequired for the suppression of P is about 12 T, a little larger than the experimental value 共⬃5 – 9 T兲for DyMnO₃ and TbMnO₃, but still in a compa- rable range.

The simulated results presented above refer to the re- sponse of uniquely orientedPtoHat different temperatures.

However, at E= 0, four equivalent directions of P are al- lowed; it is, therefore, reasonable to expect that the FE do- mains and spin-ordered domains coexist and interactively clamped owing to the strong mutual influence between the FE and magnetic orders. To illustrate this feature, we per- form the simulation under zeroE. First, the lattice is cooled fromT= 40 K toT= 0.01 K under zeroHandE, and thenH increases gradually to H= 18 T in the isothermal condition.

As shown in Fig. 4共a兲, P is gradually suppressed 共path 1兲 down to zero around H⬃13 T. Nevertheless, upon subse- quent decrease ofHfromH= 18 T down to zero,Pdoes not return back to the initial state but evolves along path 2 with P= 0. The reason behind this anomaly is the existence of the 90° FE domains in the lattice compatible with the spiral- spin-ordered regions, as shown in Fig.4共c兲. The directions of P are normal to each other between the adjacent domains, and opposite between the cross domains. Hence, the macro- scopic FE polarization remains zero. The wave vectors for the local spin ordered regions associated with the FE do- mains exhibit similar characters. If one defines the magnetic domains by the wave vector, one may find that the magnetic domain walls are coincident with the FE ones. In other words, they are mutually clamped owing to the strong ME coupling. In addition, the 180° FE domains are observed in this work, although they are not shown here.

As shown above, the formation of multiferroic domains FIG. 2.共Color online兲Snapshots for four types of the spiral spin

and FE dipole structure of a 4⫻4 cluster atT= 0.01 K. The Mn and O ions are located at the sites of the MnO₂layer of the ideal cubic perovskite lattice. Spin of Mn ions and displacement of O ions are denoted by blue and red arrows, respectively.

FIG. 3.共Color online兲 共a兲MagnetizationMand共b兲polarization Pas a function ofTat variousH, and共c兲Mand共d兲Pas a function ofHat different T. The up arrow and down arrow in共a兲and 共b兲 denote the heating and cooling sequences, respectively, and those in 共c兲 and 共d兲 denote the varying paths ofH. Numbers in the figure denote the sequence of the MC simulations.

MULTIFERROIC RESPONSE AND CLAMPED DOMAIN… PHYSICAL REVIEW B77, 054442共2008兲

054442-3

can be realized in a properly designed system. Actually, ob- servation of the clamped FE domains and magnetic domains has been reported in recent work on multiferroic manganites.²⁷The MC simulation provides an extra proof for the experimental observation. For the present system, due to the fourfold symmetry, both the 90° and 180° domains are

allowed, while different domain configurations may be pos- sible depending on the symmetry of the system. The meta- stability of the multidomain configuration in this work, evi- denced by Fig. 4共b兲, is attributed to the exclusion of the depolarization energy in our simulation, which can be low- ered by the formation of a domain structure in compensation of the increase of the domain wall energy.

In conclusion, we have simulated the presence of FE po- larization in a 2D noncollinear multiferroic lattice as a result of the DM interaction. The response of the FE polarization and magnetization to external magnetic field at different tem- peratures has been investigated. The simulation results reveal that the FE polarization relies essentially on the spiral-spin structure. Given the sufficiently high magnetic field, the FE polarization may be destroyed due to the broken spiral-spin configuration. The interactively clamped FE domains and spiral-spin domains are demonstrated in this simple model system. The gigantic ME effect as demonstrated is in good agreement with previous experiments.

We thank the Natural Science Foundation of China共No.

50601013 and No. 10674061兲and National Key Projects for Basic Research of China 共No. 2002CB613303 and No.

2006CB921802兲.

*Author to whom any correspondence should be addressed;

[email protected]

1M. Fiebig, J. Phys. D 38, R123共2005兲.

2J. F. Scott, Nat. Mater. 6, 256共2007兲.

3W. Eerenstein, N. D. Mathur, and J. F. Scott, Nature 共London兲 442, 759共2006兲.

4T. Kimura, T. Goto, H. Shintani, K. Ishizaka, T. Arima, and Y.

Tokura, Nature共London兲 426, 55共2003兲.

5T. Goto, T. Kimura, G. Lawes, A. P. Ramirez, and Y. Tokura, Phys. Rev. Lett. 92, 257201共2004兲.

6T. Arima, T. Goto, Y. Yamasaki, S. Miyasaka, K. Ishii, M.

Tsubota, T. Inami, Y. Murakami, and Y. Tokura, Phys. Rev. B 72, 100102共R兲 共2005兲.

7N. Hur, S. Park, P. A. Sharma, J. S. Ahn, S. Guha, and S. W.

Cheong, Nature共London兲 429, 392共2004兲.

8L. C. Chapon, P. G. Radaelli, G. R. Blake, S. Park, and S. W.

Cheong, Phys. Rev. Lett. 96, 097601共2006兲.

9D. Higashiyama, S. Miyasaka, and Y. Tokura, Phys. Rev. B 72, 064421共2005兲.

10G. Lawes, A. B. Harris, T. Kimura, N. Rogado, R. J. Cava, A.

Aharony, O. Entin-Wohlman, T. Yildirim, M. Kenzelmann, C.

Broholm, and A. P. Ramirez, Phys. Rev. Lett. 95, 087205 共2005兲.

11K. Taniguchi, N. Abe, T. Takenobu, Y. Iwasa, and T. Arima, Phys.

Rev. Lett. 97, 097203共2006兲.

12T. Arima, A. Tokunaga, T. Goto, H. Kimura, Y. Noda, and Y.

Tokura, Phys. Rev. Lett. 96, 097202共2006兲.

13R. Kajimoto, H. Yoshizawa, H. Shintani, T. Kimura, and Y.

Tokura, Phys. Rev. B 70, 012401共2004兲.

14M. Kenzelmann, A. B. Harris, S. Jonas, C. Broholm, J. Schefer, S. B. Kim, C. L. Zhang, S. W. Cheong, O. P. Vajk, and J. W.

Lynn, Phys. Rev. Lett. 95, 087206共2005兲.

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

16Y. Yamasaki, S. Miyasaka, Y. Kaneko, J. P. He, T. Arima, and Y.

Tokura, Phys. Rev. Lett. 96, 207204共2006兲.

17Y. Tokura, Science 312, 1481共2006兲.

18I. A. Sergienko and E. Dagotto, Phys. Rev. B 73, 094434共2006兲.

19A. B. Harris, T. Yildirim, A. Aharony, and O. Entin-Wohlman, Phys. Rev. B 73, 184433共2006兲.

20M. Mostovoy, Phys. Rev. Lett. 96, 067601共2006兲.

21J. Laverdiere, S. Jandl, A. A. Mukhin, V. Y. Ivanov, V. G. Ivanov, and M. N. Iliev, Phys. Rev. B 73, 214301共2006兲.

22E. Dagotto, T. Hotta, and A. Moreo, Phys. Rep. 344, 1共2001兲.

23T. Hotta, M. Moraghebi, A. Feiguin, A. Moreo, S. Yunoki, and E.

Dagotto, Phys. Rev. Lett. 90, 247203共2003兲.

24J. S. Zhou and J. B. Goodenough, Phys. Rev. Lett. 96, 247202 共2006兲.

25M. Tovar, G. Alejandro, A. Butera, A. Caneiro, M. T. Causa, F.

Prado, and R. D. Sánchez, Phys. Rev. B 60, 10199共1999兲.

26F. Moussa, M. Hennion, J. Rodriguez-Carvajal, H. Moudden, L.

Pinsard, and A. Revcolevschi, Phys. Rev. B 54, 15149共1996兲.

27M. Fiebig, T. Lottermoser, D. Frohlich, A. V. Goitsev, and R. V.

Pisarev, Nature共London兲 419, 818共2002兲. FIG. 4. 共Color online兲 H-dependence of共a兲 P and共b兲 average

Hamiltonian per site atT= 0.01 K.共c兲Typical snapshot of the 90°

domains separated by dashed lines. Numbers in共a兲 and共b兲denote the sequence of the MC simulations.

LI, DONG, AND LIU PHYSICAL REVIEW B77, 054442共2008兲

054442-4