Ke Xu(徐可) Junsheng Feng(馮俊生) and Hongjun Xiang(向紅軍)

1Hubei Key Laboratory of Low Dimensional Optoelectronic Materials and Devices,School of Physics and Electronics Engineering,Hubei University of Arts and Science,Xiangyang 441053,China

2School of Physics and Material Engineering,Hefei Normal University,Hefei 230601,China

3Key Laboratory of Computational Physical Sciences(Ministry of Education),State Key Laboratory of Surface Physics,and Department of Physics,Fudan University,Shanghai 200433,China

Keywords: magnets, ferroelectrics, multiferroics, Monte Carlo simulation, four-state method, DFT calculation,Property Analysis and Simulation Package for materials(PASP)software

1. Introduction

Magnetic materials are widely used in a diverse range of applications, including information storage, electronics,and even in biomedicine.[1]The collinear magnetic materials,e.g., ferromagnet (FM) and antiferromagnet (AFM),have been studied for a long time. However, the exotic magnetic states attract a lot of interest recently in both theoretical and experimental studies, such as spin glasses,[2]spin ice,[3]spin liquid,[4,5]skyrmions,[6,7]and other topological spin textures.[8]In 2017, the ferromagnetism in monolayer Cr2Ge2Te6[9]and CrI3[10]were reported experimentally. After that, the two-dimensional (2D) magnetic materials has received considerable attention, and several related studies have been carried out on the few layers CrI3and Cr2Ge2Te6.[11-14]Other 2D FM materials were also predicted theoretically or synthesized experimentally, such as Fe3GeTe2,[15]VSe2,[16]MnSe2,[17]Fe3P,[18]VI3,[19]and van der Waals(vdW)structures.[20]

Similar to magnetism, ferroelectricity is also an appealing property because it can be utilized to make functional devices, especially electronic devices such as nonvolatile memory, tunable capacitor, solar cell and tunnel junction.[21]Ferroelectric (FE) materials are a large family, including inorganic oxides,[22]organic-inorganic hybrid materials,[23]organic compounds,[24]liquid crystals,[25]and polymers.[26]In the past decade, with the rapid development of thin film technologies, 2D FEs have attracted much attention. Ferroelectricity in 2D materials with stable layered structures and reduced surface energy may be stable even in the presence of a depolarization field, thereby opening a pathway to explore low dimensional ferroelectricity. So far, various 2D van der Waals (vdW) materials with the thickness of few layers or even monolayer have been reported to be intrinsically FE, such asMX2(M= W Mo,X= S,Se, Te) transition metal dichalcogenides (TMDs),[27]group-IV monochalcogenides,[28]metal triphosphates,[29]layer perovskites,[30]indium selenide (In2Se3),[31]and sliding FEs.[32]

Magnetism and ferroelectricity are like the front and back of a coin,which can be combined in a single material dubbed as multiferroic. Multiferroic materials may be ideal candidates for applications in novel multifunctional magnetoelectric (ME) devices and high performance information storage and processing devices.[33]Two breakthrough studies in 2003 bring up the renaissance of multiferroics.In one study,the epitaxial BiFeO3thin film is successfully prepared and demonstrated the room-temperature multiferroic behavior.[34]In another study,Tokura and coworkers discovered the colossal ME effect in TbMnO3.[35]Then multiferroics are classified into two types (i.e., type-I and type-II) according to the different mechanisms of ferroelectric polarization by Khomskii.[36]In fact,the ME effect in type-I multiferroics(e.g.,BiFeO3)[34]is usually weak,since the ferroelectricity and magnetism are independent of each other. While in type-II multiferroics (e.g.,TbMnO3),[35]the ME effect is expected to be much stronger,because the ferroelectricity is induced directly by a special spin order. Very recently,the 2D ferroelastic-ferroelectric,[37]and 2D ferroelectric-ferromagnetic materials[38-40]were also reported. Recently, Songet al.reported[41]the first experimental discovery of type-II multiferroic order in a 2D vdW material NiI2.

In the past decades, first principle based studies have made significant contributions to the field of magnetism and ferroelectricity. Density functional theory (DFT) calculations have been adopted to investigating various topics including the magnetic domain wall induced electric polarization in Gd3Fe5O12,[42]skyrmions in Ni-halide monolayer,[43]Berezinskii-Kosterlitz-Thouless (BKT) phase existing between the FE and paraelectric (PE) states in BaTiO3ultrathin film[44]and single layer SnTe,[45]the proximate quantum spin liquid states in X[Pd(dmit)2]2,[46]frustrated magnetism in Mott insulating (V1-xCrx)2O3[47]and topological surface states in superlatticelike MnBi2Te4/(Bi2Te3)n.[48]In addition,first principle based effective Hamiltonian models were now widely adopted to investigate the ground state at zero temperature and temperature dependent properties.[44,45,49-52]For example,with the effective Hamiltonian model approach,Bellaicheet al.predicted a magnetic phase diagram in BiFeO3as a function of first- and second-nearest-neighbor interaction strength and find new spin arrangement.[49]The effective Hamiltonian method was also used to investigate the FE-antiferrodistortive coupling in small tolerance factor perovskite.[51]

This remaining sections of this review are organized as follows. In Section 2, we will briefly introduce the microscopic models about magnetism,ferroelectricity,multiferroicity and how to calculate ME coupling tensorαi j. Next,in Section 3 we will introduce the four-state method to compute the various parameters and our Property Analysis and Simulation Package for materials (PASP). In Section 4, we will demonstrate the applications of the models, methods and PASP to various magnetic and ferroelectric systems.[51,53,54]At last,we summarize the review and give a perspective on future computational developments in the field of magnetism and ferroelectricity.

2. Theoretical models

Although DFT played an important role on studying magnetic and ferroelectric materials, it has some limitations: it is computationally demanding and thus is not suitable for largescale systems; it has difficult to consider effects of temperature and external field; it gives the total results, but does not directly provide microscopic insight. Therefore, many theoretical models were developed to address these issues. In particular,effective Hamiltonian models were proposed to understand the microscopic mechanisms of the magnetism, ferroelectricity, and multiferroicity, and to simulate the behaviors of ferroic materials. In this section, we will first discuss the effective spin Hamiltonian models for magnetic materials.Then, we briefly introduce the effective Hamiltonian model for ferroelectrics,which provides an efficient way for computing the thermodynamic properties and simulating the dynamics behaviors of ferroelectrics. Finally, we will discuss two theoretical models related to multiferroics.

2.1. Spin Hamiltonian models

The effective magnetic Hamiltonian can be written as

equation (1) includes the bilinear interactionsHexandHSIA,and the higher order interactionsHhigher. The summation indices〈ijk〉and〈ijkl〉mean that the sum goes over triangles of sites and the sum is take over all distinct four-spin cycles on square plaquettes, respectively. For the bilinear interactions,the Heisenberg exchange parametersJij, the Dzyaloshinskii-Moriya (DM) interaction[55]vectorDi jand the Kitaev interaction parametersKγcan be extracted from the fullJmatrix (as will be defined in Subsection 2.1.2). The single ion anisotropy (SIA) Hamiltonian is characterized by theAmatrix. The calculation of these magnetic parameters will be discussed in Subsection 3.1.

2.1.1. Heisenberg model

By applying the perturbation approach to the Hubbard model,[56]the effective spin Hamiltonians such as celebrated simplest or isotropic Heisenberg model can be obtained by down-folding the fermionic degrees of freedom into the proper low energy subspaces.[56]The Heisenberg model can be written as

For magnetic systems, the Heisenberg spin exchange constants in a given system can be calculated by the four-state energy mapping method within the DFT framework as discussed below.[57,58]WhenJij ＞0 andJij ＜0, antiferromagnetic (AFM) and ferromagnetic (FM) configuration are preferred.

2.1.2. Kitaev model

The exactly solvable Kitaev model[59]is a bonddependent interaction of spin-1/2 model on a 2D honeycomb lattice,in which the spins fractionalize into Majorana fermions and form a topological quantum spin liquid (QSL) in the ground state.α-RuCl3[60]is a 2D magnetic insulators with a honeycomb structure, and it was noticed that they accommodate essential ingredients of the Kitaev model owing to the interplay of electron correlation and spin-orbit coupling(SOC). The Kitaev model consists of spin-1/2 on a honeycomb lattice, with nearest-neighbour Ising interactions with bond-dependent easy axes parallel to thex,y, andzaxes.The orthogonal anisotropies of the three nearest-neighbour bonds of each spin conflict with each other,leading to the spin frustration.[61]In 2018,Xuet al.predicted the existence of Kitaev interactions in 2D CrI3and CrGeTe3-like systems.[53]In 2019, Hae-Young Keeet al. presented a theory of the spin-1 Kitaev interaction in 2D edge-shared octahedral systems.[62]It was recently proposed that spin-3/2 QSL can be realized in 3d transition metal compound CrSiTe3.[63]

The Kitaev interaction originate spin-orbit coupling(SOC)effect,and its terms can be derived from theJmatrix.TheJmatrix characterizing the magnetic exchange couplings are expressed in the most general 3×3 matrix form as

2.1.3. Higher order interactions

Beyond the bilinear interactions, higher order exchange interactions involving the hopping of two or more electrons sometimes play a crucial role in exotic magnetic materials,[64]especially if the magnetic atoms have large magnetic moments or if the system is itinerant.[1]Applying the fourth order perturbation theory with the L¨owdin’s downfolding technique to the multiorbital Hubbard model, the fourth order interaction can be derived.[56]The fourth order interactions consist twobody fourth-order interactionKi j(Si·Sj)2[i.e., biquadratic exchange (BQ) interaction], three-body fourth-order interactionKijk(Si·Sj)(Sj·Sk)and four-body fourth-order interactionKijkl(Si·Sj)(Sk·Sl)[see Eq.(1)].

Higher order interactions are found to be crucial to the magnetic features in several systems, such as multilayer materials,[65]perovskites,[66,67]and 2D magnets.[68,69]Orthorhombic perovskite o-HoMnO3[66]is reported to be E-type antiferromagnetic (E-AFM) in the experiment. A theoretical study suggested that the E-AFM state can be properly predicted only when the BQ exchange term is included in the model Hamiltonian.[67]In 2D vdW magnets, the higher order magnetic interactions were also found to play a key role on the magnetic properties.[64]In the single layer NiCl2, the nearest-neighbor BQ interaction is necessary for FM ground state prediction.[68]Furthermore,the four site four spin interaction term was found to have a large effect on the energy barrier preventing skyrmion and antiskyrmion collapse into the FM state in Pd/Fe/Rh(111)ultrathin film.[70]

2.2. Effective Hamiltonian models for ferroelectrics

Following the approach of Vanderbilt and his coworkers,[71]the effective Hamiltonian is written in terms of the FE degree of freedom{u}and the homogeneous strains

Parameters in such an effective Hamiltonian can be fitted by the energy mapping method based on DFT calculations with numerous configurations. The effective Hamiltonian method combined with the MC simulations can not only predict FE phase transitions from the temperature dependence of〈u〉and〈η〉,[48,50]but also construct the temperaturepressure phase diagram.[72]

2.3. Theoretical models for multiferroics

As we mentioned above,there exist both magnetism and ferroelectricity in a single multiferroic system. The coupling between magnetism and ferroelectricity is particularly intriguing as it not only is of great scientific importance, but also provides the basis for making novel energy-saving memory devices.Here,we discuss two theoretical models related to the ME coupling. In the first model,we describe a unified model for spin order induced improper ferroelectric polarization,which can explain the polarization induced by any spin structures: collinear,cycloidal spiral,proper screw,etc.[52,73-76]In the second model, we provide a new model for investigating the linear ME coupling effect [i.e., the first order change of magnetization(or electric polarization)response to the external electric(or magnetic)field].

2.3.1. Unified polarization model for spin-order induced ferroelectricity

For type-II multiferroics,the total electric polarizationPtcan be viewed as a function of the spin directionsSiof the magnetic ionsi,the displacementsukof the ionskand homogeneous strainηl:[54,74-76]the pure electronic contributionPearises from the electron density redistribution induced by the spin order,and the lattice deformation contributionPion,latticeoriginates from the spin order induced ion displacements and stress.

To obtain the pure electronic contributionPe,let us consider a spin dimer with spatial inversion symmetry at the center, the distance vector from spin 1 to spin 2 will be taken along thexaxis. When a noncollinear spin configuration is applied on the dimer, the inversion symmetry will be broken,hence ferroelectric polarizationPwill be induced. In general,the polarizationPis a function of the directions of spin 1 and spin 2,i.e.,P=P(S1x,S1y,S1z;S2x,S2y,S2z). ThePcan be expanded as a Taylor series ofSiα(i=1,2;α=x,y,z). The odd terms of the Taylor expansion should vanish due to the time reversal symmetry. To neglect the fourth and higher order terms,thePis written as

wherePintis a matrix in which each element is a vector. For easily understanding, thePintcan be decomposed into the isotropic symmetric diagonal matrixPJ,antisymmetric matrixPD,and anisotropic symmetric matrixPΓ:

Besides the spin order induced atomic displacements,the spin order induced lattice strain may give rise to an additional electric polarization through the piezoelectric effect.[78,79]In general, the total energy of a magnetic system can be written asE(um,ηj,Si)=EPM(um,ηj)+Espin(um,ηj,Si),whereumis the ion displacement of reference structure for all the ions,ηjis the homogeneous strain in Voigt notation, andSiis the spin vector of all the magnetic ions. With neglecting the terms of higher order, the Taylor series expansion of the paramagnetic state energyEPM(um,ηj)is

The first order coefficientsAmandAjrefer to the force and stress, respectively. The second order coefficientsBmn,Bjk, andBmjrepresent the force constant, frozen-ion elastic constant, and internal displacement tensor, respectively.The reference structure is in equilibrium in the PM state:?EPM/?um=?EPM/?ηl=0 withAm=Aj=0.

The ion displacement and lattice deformation caused by the spin order can be obtained by minimizing the total energyE(um,ηj,Si) with respect toumandηj, that is,?E(um,ηj,Si)/?um=0 and?E(um,ηj,Si)/?ηj=0. If the system in the PM state is piezoelectric, the lattice deformation induced by spin order may give rise to an additional electric polarization. The polarization induced by the ion displacementumand strainηjcan also be computed withPα=Zαmum+eα jηj,whereZαmandeα jare the Born effective charge and frozen-ion piezoelectric tensor, respectively.This unified model can not only describe the spin-lattice coupling in 2D CrI3,CrGeTe3[80]and monoxide NiO,[81]but also explain the ferroelectricity of the 2D multiferroic NiI2.[41]

2.3.2. A model for linear magnetoelectric coupling

The linear ME effect refers to a linear response of electric polarization ΔPto an applied magnetic fieldH, or the magnetization ΔMinduced by an electric fieldE,i.e.ΔPi=αijHjand ΔMj=αijEi,whereαis the ME tensor.[82]The underlying mechanism of the linear ME effect comes from the breaking of spatial inversion symmetry (I) and time reversal symmetry (T). Note that a linear ME material (e.g., Cr2O3) may display the joint symmetry TI of spatial inversion symmetry and time reversal symmetry. The linear ME effect are found in Cr2O3,[83]FeS,[84]NaMnF3,[85]BiFeO3,[86]double layered CrI3,[5]and perovskite superlattice.[87]The ME coupling tensor consists of two parts: (i)a“clamped-ion”contribution that accounts for ME effects occurring in an applied field with all ionic degrees of freedom remaining frozen, which is labelαel;(ii)lattice relaxation in response to the external field,namedαlatt. The total response tensor isαtot=αel+αlatt.Although the linear ME tensor has 9 components, the nonzero components depend on the magnetic symmetry of the system.[83-87]The methodology to calculate ME tensor was implemented by modifying the code in VASP package through calculating the electric polarization induced by an applied ZeemanHfield which only couples to the spin component of the magnetization.[83]Here, the response of electron spin is obtained by “clamping” the ions during the calculation; the relaxation of ionic positions in response to theHfield yields the sum of the ionic and electron spin components. The magnetic field is applied self-consistently and the calculations are done with setting a non-collinear spin configuration and including SOC effect. Actually, it is very challenging to calculate the ME coupling coefficient by applied a magnetic field. Because the forces induced by the applied magnetic field on the ions are rather small, the very high convergence criterion of the ionic forces is needed. In order to obtain the accurate ME response,the ionic forces need to be reduced to less than 5 μeV·?A-1(Ref.[87])in practical calculations.

Usually, an external field (εorH) needs to be applied when calculating the ME coupling coefficient. In a pioneering study by Spaldin and coworkers,the polarization was calculated self-consistently by adding a Zeeman magnetic field to the DFT Kohn-Sham potential.[83]The magnetic field will change the spin orientations of magnetic ions, resulting in a change in the electric polarization due to the ME coupling.However, the calculation of the ME response by applying a magnetic field in first-principles calculations requires extremely high precision and this is very time-consuming. Recently,Chaiet al.employed the unified polarization model to the case in the presence of an external magnetic field.[88]By performing magnetic symmetry analysis and calculating the spin susceptibility tensorχm,linear and quadratic ME effects can be obtained. However,a first-principles scheme based on this method is still lacking. In this review,we propose for the first time a new method based on DFT calculations for computing the ME coupling parameters without applying an external field in the DFT calculations. First,we compute the ionic displacement and strain induced by a small external electric field. Then we can estimate the change of the exchange interactions induced by the effect of the external electric field.Finally,we can estimate the change of magnetization by finding the magnetic ground state of the system under the external electric field.

The total energyEis a function of electric fieldε, ionic displacementu, and homogenous strainηwith respect to a reference system(ε=u=η=0),and it can be expanded up to second order with respect to(ε,u,η):

The equilibrium conditions require that the first derivatives of total energyEregarding to the ionic displacementuiαand strainηjare zero:

To calculate the dJ/dεin Eq. (17), we need calculate three first derivative parameters?J/?ε,?J/?uiα, and?J/?ηifirstly (see Subsection 3.1 for the calculation details). After the exchange interaction parameterJobtained under a small external electric fieldε,the total magnetizationM(ε)can be achieved by performing MC simulation or conjugate gradient(CG) minimization. Hence, the ME coupling parameters can be deduced withαij=?Mi(ε)/?εj.

The advantage of our method is that no external field is required, and the ME coupling tensor can be calculated only by calculating the corresponding parameters through the DFT method. Note that the dJ/dεconsists of three parts, the first part can be regarded as the electron contribution:?J/?ε,and the other two parts (?J/?uiαand?J/?ηi) are lattice contributions. Therefore, the ME tensorαcalculated through this method can also be divided into the contributions of electronsαeland ionsαlattfor further analysis.

3. Computational methods and PASP code

3.1. Four-state method

3.1.1. Second-order magnetic parameters

For isotropic Heisenberg exchange interactionJ(without SOC), it can be deduced through setting the following four spin configurations ofi-th andj-th magnetic ions: (i)spin up states for bothSiandSj: (↑,↑); (ii) spin up state forSiand spin down state forSj: (↑,↓);(iii)spin down state forSiand spin up state forSj: (↓,↑); (iv) spin down states for bothSiandSj: (↓,↓). The other spins are set according to experimental spin state or a low-energy state and remain unchanged in the four spin states. By calculating the total energies, the exchange interactionJ(neglect the subscriptijfor short) can be determined byJ=(E1+E4-E2-E3)/4. In addition,the bilinear interactions with the effect of SOC can also be calculated by this four-state method in a similar way,such asJand A matrices.[57]

3.1.2. Parameters in unified polarization model

3.1.3. Parameters of linear ME effect

DFT calculations can give atomic forces and stress for a given structure and a given magnetic configuration. The four spin states setting on thei-th andj-th magnetic ions are same as in the case of the isotropic HeisenbergJ. Differentiating the exchange parameterJwith respect to electric fieldε,ionic displacementukαand strainηk,the corresponding parameters of linear ME effect can be obtained with the following expression:?J/?ε=-(P1+P4-P2-P3)/4,?J/?ukα=-(F1+F4-F2-F3)/4, and?J/?ηk=-(σ1+σ4-σ2-σ3)/4.[75]

3.2. Machine learning method for constructing effective Hamiltonian

In principle, all the magnetic parameters can be derived by the energy mapping method. The four-state method can be used to calculated the bilinear interactions and biquadratic interactions, but it has a limitation for calculating the higher order interactions. However,the higher order interactions can not be neglect in some situations when describing the physical properties,especially in metallic and narrow gap magnets.Thus, a machine learning (ML) approach[89]was developed

where the important interaction parameters can be extracted with multiple linear regression analyses and adopting several ML techniques. In our ML method,we first set manually the truncation distance and the highest order of the interactions.Then,we will obtain all possible interactions according to the truncation order and truncation distance by performing group theory analysis. The important interactions are selected automatically with the ML approach. The decision of the truncation distance and order is mainly based on experience and a small amount of testing. For magnetic systems, it is usually enough to consider spin interactions up to the fourth order within the distance of 10 ?A.For ferroelectrics,it may be necessary to keep the interactions up to the sixth order.

The ML approach is a general approach to construct the effective Hamiltonian models. For example,the ML approach has been well applied to study the magnetic interactions in the multiferroic material TbMnO3,[90]monolayers NiX2(X=C,Br, I)[68]and Fe3GeTe2.[69]Our ML approach can not only reproduce previous results, but also can predict the other important magnetic interactions,such as three-body fourth-order interactions.

3.3. Introduction to PASP

In order to study the properties of complex condensed matter systems and overcome the insufficiency of the DFT calculations (e.g., size limit and unable to calculate finite temperature properties), we have developed a software package which is named PASP(Property Analysis and Simulation Package for materials).[90]Almost all models and methods discussed above were implemented in PASP.Based on the first principle calculations,PASP is able to simulate the thermodynamic properties of complex systems and provide insight into the microscopic mechanisms of the coupling between multiple degrees of freedoms.

The PASP package can be used in conjunction with first principle calculation packages,e.g., VASP and Quantum Espresso (QE),etc. Before the DFT calculations, our package can automatically generate a general form of the Hamiltonian according to the symmetry analyses and can generate a set number of configurations with a given supercell. Then these configurations will be prepared for running the first principle calculation. When the DFT calculations done, the parameters in Hamiltonian can be fitted with PASP code and the effective Hamiltonians will be constructed.Based on this effective Hamiltonian,the PASP package can be adopted to perform parallel tempering Monte Carlo (PTMC)simulation[91,92]which can predict the ground state and thermodynamic properties. Our current PASP package includes several functionalities as follows:

(I)Symmetry and group theory analysis

The PASP software can identify the symmetry(including point group,space group,and magnetic group)of the system.It can also identify the irreducible representation (IR) of the Bloch wave function, which is convenient for further analysis, such as constructing thek·pHamiltonian. In addition,we added some modules to facilitate the usage of four-state method(Subsection 3.1). All nonequivalent magnetic pair interactions within a given cutoff radius can be identified automatically. For each magnetic pair, PASP will output the four magnetic configurations for the use in subsequent DFT calculations.

(II)Global structure searching

In our package, we implemented the basin-hopping method[93]and genetic algorithms (GA) to predict the structures of clusters, 2D and 3D crystals, and interfaces. Recently, we implement the GA based global optimization approach with explicitly considering the magnetic degree of freedom.[94]

(III)Effective Hamiltonian methods

The effective Hamiltonian method can be applied to study the different kinds of physical properties. Correspondingly,several different types of effective Hamiltonian methods are implemented in this module. (i)Tight binding(TB)Hamiltonian: The TB model is often used to calculate the electronic states of a material and understand the mechanism of magnetism and ferroelectricity. (ii) Effective atomistic Hamiltonian: In PASP,we implemented the effective Hamiltonian approach [see Section 2, Eq. (1)] including the spin Hamiltonian model (Subsection 2.1), and effective Hamiltonian for ferroelectrics (Subsection 2.2). Besides, the unified polarization model is implemented in order to compute the spin order induced polarization (Subsection 2.3). (iii) Machine learning method for constructing realistic effective Hamiltonian: Besides the simple four-state method,the ML method discussed in Subsection 3.2 is integrated to PASP to calculate the higherorder magnetic interactions (Subsubsection 2.1.4). (iv) Machine learning potential:Instead of the usual polynomial form,here the artificial neural network(ANN)is used to describe the potential-energy surface of magnets and ferroelectrics.

(IV)Monte Carlo simulation module

In our package, we implement not only the usual Metropolis MC algorithm but also the parallel-tempering(PT) algorithm, which can be regarded as a parallelized version of the simulated tempering but with different extended ensemble.[95,96]Our effective PTMC method was adopted to predict the ground state and phase transition in different kinds of systems.[50-52,69]

4. Applications

In this section, the application of the effective Hamiltonian method in the PASP package to three typical systems(2D FE SnTe,[50]the 2D FM CrI3and CrGeTe3[53]and the layered multiferroic MnI2[54]) will be described to demonstrate the power of the method and software.

4.1. 2D ferroelectric SnTe

In the perfectly defect-free SnTe thin films, in order to investigate howTcintrinsically changes with the layer thickness,an effective Hamiltonian is constructed to estimate itsTc.A simple bulk structure is considered where the dipoles form a tetragonal lattice and the FE polarization is along the in-plane[110]direction.

Fig. 1. Tc as a function of film thickness estimated with the simple effective Hamiltonian. The bulk structure and 3-layer thin film structure are shown in panels (a) and (b), respectively. The results for ksur 2 =k2 =0.05 and ksur2 =-0.1 ＜k2 =0.05 are shown in panels (c)and(d),respectively. The other parameters are k4=0.4 and J=-0.1.The horizontal dashed line represents the Tc for the bulk structure. Reproduced with permission from Ref.[50].

4.2. 2D ferromagnetic CrI3 and CrGeTe3

CrI3and CrGeTe3monolayers were discovered to be ferromagnetic.[9-14]The valence state of chromium ion in these two compounds are all +3, with the 3d3configuration andS= 3/2. The FM state is induced by the super exchange interaction between nearest-neighbor Cr ions which are linked by I or Te ligands with the nearly 90°angles. CrI3has been demonstrated to be well described by the Ising behavior with the spins pointing parallel or antiparallel to the out-of-planezdirection. In contrast, the magnetic anisotropy of CrGeTe3was determined to be consistent with the Heisenberg behavior,for which the spins can freely rotate and adopt any direction in the three-dimensional space. The Hamiltonians to describe magnetic properties in CrI3or CrGeTe3is:H=Hex+HSIA=∑〈i j〉Si·Ji j·Sj+∑i Si·Aii·Si, whereJijandAiiare 3×3 matrices gathering exchange interaction and SIA parameters, respectively. Three different coordinate systems{xyz},{XYZ},and{αβγ}are shown in Fig.2. In CrI3,the calculatedJmatrix is symmetric,and theJxx,Jyy,andJzzare-2.29,-1.93, and-2.23 meV, respectively. The diagonal formJα,Jβ,andJγare-2.46,-2.41,and-1.59 respectively. TheJαandJβin CrGeTe3are about-6.65 meV,whileJγ=-6.28 meV is about 0.4 meV smaller in magnitude than the other two exchange coefficients.

Fig. 2. Schematization of the CrI3 and CrGeTe3 structures, as well as the different coordinate systems indicated in the text. The planes in blue,green and red indicate the easy plane form Kitaev interaction for Cr0-Cr1,Cr0-Cr2,and Cr0-Cr3 pairs, respectively. Note that Ge of CrGeTe3 is not shown for simplicity. Reproduced with permission from Ref.[53].

4.3. Layered multiferroics MnI2

Fig. 3. (a) The 5×5×1 supercell of MnI2. The left inset illustrates the layered structure of MnI2. (b) The electric polarizations predicted by the KNB and gKNB models for three different spin configurations of the Mn-Mn dimer. (c)The polarization of the Mn-Mn pairs extracted by the gKNB model(line)and direct density functional calculation(dot).Reproduced with permission from Ref.[54].

5. Summary and perspective

In this review, we summarized the recent theoretical and computational advances in magnetics,ferroelectrics and multiferroics. We first discuss the effective Hamiltonian of magnetics and ferroelectrics, and then describe the unified polarization model for multiferroics and the method for calculating the ME coupling tensor in details. Next,we discussed the“four-state”method within the DFT framework,and the PASP package, and present three examples to illustrate the typical applications of the PASP package.

Although first principle based approaches have made important progress in understanding the mechanisms of magnetism and ferroelectricity,there are still big challenges to be overcome when simulating realistic materials and predicting new materials. For example, the magnetic properties depend on the adopted HubbardUparameter which is hard to estimate accurately; local density approximation (LDA) will usually underestimate the ferroelectricity, while the widely adopted generalized gradient approximation Perdew-Burke-Ernzerhof(PBE) functional will overestimate the ferroelectricity. DFT calculations cannot directly simulate the flipping dynamics of ferroelectric, magnetic and multiferroic domains under external fields. These issues can be solved by developing new DFT functionals or developing more efficient and more accurate effective Hamiltonian approaches.

Besides the limits of current DFT in studying ferroic materials to be overcome, in our opinion some other research directions in this field remains to be explored: (I) Predicting and designing high-performance(especially high critical temperature, low dimensional)ferroic materials is still an important topic; (II) Rare earth compounds may host exotic states(e.g.,quantum spin liquid,multiferroicity,or heavy-fermions),while reliable model Hamiltonians applicable to the f-electron(lanthanides and actinides)systems are lacking;(III)The interaction between strong optical fields and ferroic materials may give rise to unusual phenomena, which deserve much more research attention.

Acknowledgements

Project supported by the National Natural Science Foundation of China (Grant Nos. 11825403, 12188101,and 11804138), the Natural Science Foundation of Anhui Province,China(Grant No.1908085MA10),and the Opening Foundation of the State Key Laboratory of Surface Physics of Fudan University (Grant No. KF201907). We thank Dr.Zhang H M for useful discussions.