Mesoscopic scale

Early works on modeling and predicting the response of ME materials were carried out using continuum-medium methods at the mesoscale, which not only provided quantitative understanding of the direct ME responses in bulk ME composites, but also predicted new composites with a giant ME effect.^{Reference Nan, Bichurin, Dong, Viehland and Srinivasan2} In recent years, mesoscale modeling of ME materials has focused on nanocomposites, especially heterostructures integrating ferroelectric and magnetic materials. Magnetization switching or domain structure in ME heterostructures can be directly controlled by applying an electric field, offering attractive possibilities for novel potential ME devices such as memories, spintronic devices, and electrically tunable microwave devices.

On the mesoscale, the phase-field method is powerful for understanding electric-field control of magnetism by investigating complex magnetic/ferroelectric domain evolution processes, and thereby designing novel ME devices.^{Reference Hu, Duan, Nan and Chen7,Reference Hu, Chen and Nan8} For example, based on phase-field simulations of electric-field-induced magnetization switching via strain-mediated ME coupling, the combination of a spin-valve or magnetic tunnel junction (a key component in most spintronic devices) with a ferroelectric layer has been designed to construct an electric-field actuated magnetic random-access memory.^{Reference Hu, Li, Chen and Nan9}

Recently, the micromagnetic method combined with the phase-field method (shown in Figure 2) has been further developed to explore the dynamic process of creation/release of strains in ferroelectric materials and strain-enabled magnetization switching in magnets. The temporal evolution of the magnetization vectors in the ME heterostructures is obtained numerically by solving the Landau–Lifshitz–Gilbert (LLG) equation with the total effective field, H _eff. This field is determined by the minimization of free energy, which not only includes magnetic exchange energy and stray field energy, but also includes the elastic energy of the heterostructure system. Using this phase-field method, full 180° reversal of in-plane or perpendicular magnetization in designed ME heterostructures was achieved purely by a perpendicular electric-field pulse^{Reference Hu, Yang, Wang, Huang, Zhang, Chen and Nan10–Reference Peng, Hu, Chen and Nan12} (Figure 2). Such 180° reversal of magnetization driven by an electric field is highly desired and is critical for the efficiency of beyond complementary metal oxide semiconductor (CMOS) devices.^{Reference Manipatruni, Nikonov and Young1}

Figure 2.

A 180° reversal of magnetization is observed in (a) heterostructures integrating ferroelectric and magnetic films. (b) In-plane strain, ε_p, induced by a voltage pulse (shown by the gray step functions) applied to the ferroelectric transfers to a magnet, and thus the magnetization direction reverses (see the direction changes shown by the arrows from up with red color to down with blue color in [b]). (c) The strain is calculated using the finite element method (FEM) and used as an input for phase-field (PF) and Landau–Lifshitz–Gilbert (LLG) calculations in the magnet. Note: M, magnetization; t, time.

Micromagnetic phase-field simulations have also been combined with the finite element method to analyze the speed of the piezostrain-mediated voltage-driven magnetization reversal.^{Reference Peng, Hu, Chen and Nan12} The strain was modeled by solving the elastodynamics equation in a three-dimensional finite element model with a nonzero piezoelectric stiffness damping coefficient, which was then fed into a micromagnetic phase-field model for simulating strain-enabled magnetization switching in the magnet. The homogeneous strain in the micromagnetic phase-field model was set to be equal to the output of the time-dependent finite element model, which is the average local piezostrain on the ferroelectric film. This computational framework can be generally used to analyze the speed of piezostrain-enabled magnetization switching and magnetic domain-wall motion, which critically determines the response time of the corresponding ME devices.

One of the challenges for modeling is how to take into account the effects of realistic topological structures on multilength scales, devices with irregular geometries, and defects (such as grain boundaries) in polycrystalline materials. Miehe et al.^{Reference Miehe and Ethiraj13} presented a geometrically consistent rate-type variation formulation for phase-field models, where the evolution of the displacement, magnetization, and magnetic potential fields were defined, and the magnetization direction was considered as a geometric object with its rotation nature preserved to study structures with irregular geometries or under complex loading conditions.

In polycrystalline materials, the effects of grain size on the stress-related polarization behavior are discussed by combining phase-field theory and finite element methods. However, atomistic topological defects, such as dislocations, are beyond the traditional phase-field approaches, which typically lack an explicit connection to the crystal structure features in the atomic scale. This could be addressed by introducing an atomic number density-related periodic order parameter to the so-called phase-field crystal method.^{Reference Seymour, Sanches, Elder and Provatas14–Reference Asadi, Asle Zaeem and Baskes16} With this model, it is possible to explore the interaction of defect structures with external magnetic or electric fields, on diffusion time scales, which can be instructive for microstructural engineering via external fields in ME materials.

Atomistic scale

At the atomistic scale, the so-called effective Hamiltonian approach has been successful. It can predict the low-energy states of multiferroic ME materials under a variety of conditions (e.g., BiFeO₃ [BFO] thin films^{Reference Prosandeev, Kornev and Bellaiche17} or nanodots^{Reference Ren and Bellaiche18}) at finite temperatures or under an applied electric field. For a better understanding of phase-transition dynamics and domain evolution, molecular dynamics (MD) can be used along with the effective Hamiltionian approach and LLG method. This combination was first used to study the coupling effect of phonon–magnon and oxygen octahedral distortions in BFO.^{Reference Wang, Weerasinghe and Bellaiche19}

More recently, more sophisticated atomistic potentials were introduced to deal with polarization and spin separately. Two types of atomistic potentials were introduced, bond-valence potential^{Reference Brown20} and shell model potential,^{Reference Graf, Sepliarsky, Stachiotti and Tinte21} with respect to polarization. In the bond-valence potential, the nonsymmetric atomic configuration is mainly determined by the bond-valence energy and the bond-valence vector energy, which reflect the environmental bond order and local distribution of density of states.^{Reference Brown20} With this method, the phase transition of BFO^{Reference Liu, Grinberg and Rappe22} was accurately predicted and a precise description of domain dynamics can be expected in view of the existing success in domain evolution in BaTiO₃,^{Reference Qi, Liu, Grinberg and Rappe23} PbTiO₃,^{Reference Xu, Liu, Grinberg, Karthik, Damodaran, Rappe and Martin24} and other more complicated compounds.^{Reference Takenaka, Grinberg, Liu and Rappe25} In the shell model, the interatomic interactions are treated by an ion core coupled with an “electronic” shell, where the electronic polarization is included and short-range interactions as well as the three-body bond bending terms are described through an anharmonic term.^{Reference Graf, Sepliarsky and Stachiotti26} Dynamical behaviors such as terahertz-range collective vibration of the pinned domain wall^{Reference Hlinka, Pasciak, Korbel and Marton27} as well as supersonic movements of a kink inside the domain wall ^{Reference Salje, Wang, Ding and Scott28} were also discovered.

Consideration of the magnetic moment with respect to the applied stress/strain is absent from the previously discussed atomistic simulation schemes. Recently, so-called spin–lattice dynamics^{Reference Tranchida, Plimpton, Thibaudeau and Thompson29–Reference Hellsvik, Thonig, Modin, Iuşan, Bergman, Eriksson, Bergqvist and Delin31} was developed to address this issue—it mainly couples the dynamics of spin and the lattice. The spin–lattice MD starts with the spin–lattice Hamiltonian, which includes one extra spin Hamiltonian term in addition to the kinetic and potential terms in classical MD. The atom position, velocity, and spin are obtained self-consistently from the equations of motion, as well as the LLG equation. The spin configuration and lattice are coupled through the exchange interaction, or pseudo-dipolar and pseudo-quadrupolar functions.^{Reference Perera, Eisenbach, Nicholson, Stocks and Landau32–Reference Perera, Nicholson, Eisenbach, Stocks and Landau34} With this method, the spin–orbital interaction,^{Reference Perera, Eisenbach, Nicholson, Stocks and Landau32} magnetostriction,^{Reference Ma, Dudarev and Woo30} and magnon–phonon interaction^{Reference Hellsvik, Thonig, Modin, Iuşan, Bergman, Eriksson, Bergqvist and Delin31} can be investigated. Moreover, other dynamic structure factors for spin, displacement, and velocity can be obtained by calculating the space- and time-displaced pair-correlation functions as well as by its Fourier transform.

It is expected that MD can be successfully used to consider the ME effects for both pristine structures and heterostructures under external fields. However, the implementation of this method in ME materials requires the precise description of the interatomic interaction, and therefore requires assistance from the electronic scale (Figure 1).

Electronic scale

Understanding and predicting the response of single-phase ME compounds as well as superlattices rely on calculation methods at the electronic scale, such as DFT. By combining the modern polarization theory of the Berry phase (which encodes information about the phase evolution of the ground state along a certain path) and the Landau theory to expand the free energy as a power series of the variant order parameters, it is possible to identify electronic^{Reference Bousquet, Spaldin and Delaney35} and lattice contributions^{Reference Íñiguez36} to the ME coupling of these materials. Moreover, the orbital contribution can be taken into account by the recently developed orbital selective external potential method^{Reference Du, Ding, Sheng, Savrasov, Wan and Duan37} and Chern–Simons theory.^{Reference Liu and Vanderbilt38}

Supported by these methods, connections were set up between ME properties and structural distortion patterns, especially the ferroelectric, antiferroelectric, and antiferrodistortive (AFD) modes of the oxygen octahedron in perovskite ABO₃ materials, where A is a divalent or monovalent metal and occupies the cube corners, and B is a tetra- or pentavalent atom at the face centers for perfect perovskite crystals. These modes interact in different ways depending on the materials, for example, AFD and the ferroelectric mode compete in BFO^{Reference Wojdel and Íñiguez39,Reference Chen, Grisolia, Zhao, González-Vázquez, Bellaiche, Bibes, Liu and Íñiguez40} and can be modified by strain,^{Reference Xiang, Guennou, Íñiguez, Kreisel and Bellaiche41} chemical doping, or an external field.^{Reference Varignon, Bristowe and Ghosez42} To engineer functionalities, accurately treating systems in strongly correlated systems is essential. Recently, the development of dynamical mean field theory combining with DFT calculations presented a detailed description of both local and nonlocal correlations (Figure 1).^{Reference Jiang43} This was successfully implemented in the exploration of Mott multiferroic compounds, which demonstrate magnetoelectric coupling of the polar Mott insulator, such as Sr_0.5Ba_0.5MnO₃,^{Reference Nourafkan, Kotliar and Tremblay44,Reference Puggioni, Giovannetti, Capone and Rondinelli45} and seem to be promising.

The distortion modes can cooperatively interact with other degrees of freedom of the oxygen octahedron, such as the breathing modes, in-phase and anti-phase rotations, and the symmetry breaking distortion of the six-coordinate BO6 octahedral known as the Jahn–Teller effect (Figure 3). More specifically, the elements of the A and B sites can be gradually substituted in an ordered manner where symmetries break both in space and time, and novel ME coupling mechanisms were consequently introduced. For example, enhanced ferroelectric Curie temperature, T _C, and ME coupling strength were predicted by DFT calculations in A-site ordered perovskite LaMn₃Cr₄O₁₂ (AA′B₄O₁₂, where both A′ and B sites accommodate magnetic transition-metal ions),^{Reference Wang, Chai, Zhou, Cao, Cruz, Yang, Dai, Yin, Yuan, Zhang, Yu, Azuma, Shimakawa, Zhang, Dong, Sun, Jin and Long46} a Ruddlesden–Popper compound (AA′B₂O₆),^{Reference Young and Rondinelli47} and superlattices such as (A²⁺/R³⁺TiO₃, R represents rare-earth ions),^{Reference Bristowe, Varignon, Fontaine, Bousquet and Ghosez48} RFeO₃/RCrO₃,^{Reference Zhao, Bellaiche, Chen and Íñiguez49} and BiFeO₃-SrTiO₃.^{Reference Weng, Lin, Dagotto and Dong50} Not only interfaces, but also defects (e.g., dislocations) in ferroelectric PbTiO₃^{Reference Shimada, Xu, Araki, Wang and Kitamura51} can introduce novel ME coupling.

Figure 3.

Schematic of four main possible distortion modes of the oxygen octahedron in perovskite materials, ABO₃ (A = divalent or monovalent metal and B = tetra- or pentavalent atom). The coupling between these distortion modes and the degrees of freedom of charge, orbital, spin, and lattice can induce comprehensive magnetoelectric effects. The distortion modes include contraction or elongation of the equatorial B–O bond lengths due to compressive or tensile strain; the rotation of BO₆ perpendicular to the substrate with two BO₆ octahedra rotating with same angle in the same direction (in phase) and in the opposite direction (antiphase); tilting of BO₆ about an axis that is parallel to the substrate plane with four possible components labeled as 1, 2, 3, and 4; and (breathing mode) a B–O bond length disproportionation with alternating BO₆ octahedra that are all equivalent.

Although spin, orbital, charge, and lattice degree of freedom can be illustrated by DFT calculations to investigate the complex coupling mechanisms crossing interfaces, DFT cannot adequately describe ground states for interfaces, grain boundaries, or dislocations with sophisticated configurations. A so-called second-principles method^{Reference García-Fernández, Wojdeł, Íñiguez and Junquera52,Reference Escorihuela-Sayalero, Wojdeł and Íñiguez53} has been proposed to overcome this limitation as an extension of self-consistent-charge density functional tight-binding approach. With this method, large-scale electronic calculations can be carried out at a modest computational cost, where a minute energy difference on the order of meV per atom can be dealt with quickly and accurately. This method shows the capability of accurately treating complicated electronic structures as well as the accompanying lattice distortion and strong electron–lattice couplings at the interface.^{Reference García-Fernández, Wojdeł, Íñiguez and Junquera52,Reference Escorihuela-Sayalero, Wojdeł and Íñiguez53} The implementation is still limited by the assumption of a fixed bonding topology, and therefore it is restricted to systems in which it is possible to (loosely) define an underlying bonding topology that is to be preserved.

Summary and prospects

Modeling and predicting the responses of ME materials were briefly discussed using different methods such as phase-field methods at the mesoscale, molecular dynamics at the atomic scale, and first-principles calculations at the electronic scale. Despite their success, more efforts are required to bridge the gaps between these different length scales and to avoid complex transferring of parameters between different scales (see Figure 1). The other trend is the consideration of real-time dynamics, such as MD with spin or polarization and the time-dependent phase-field theory, where it is exciting to show the in situ evolution of magnetism/polarization with respect to the application of an external field, revealing the hidden mechanisms that may be absent from previous static or quasi-dynamic simulations.

Ideally, the integration of methods on different scales and the information transfer between them can help to address complicated interactions on all scales. However, this requires significant efforts in developing appropriate potentials (atomistic or phenomenological). Fortunately, the revival of artificial intelligence (AI) or so-called “machine learning” can make things simpler. For example, by using convolutional neural network methods, the fitting process of the explicit potential can be bypassed, and the atomistic behaviors can be accurately predicted by a model trained from a small set of configurations from DFT calculations.^{Reference Botu, Batra, Chapman and Ramprasad54,Reference Zhang, Han, Wang and Car55} The power of AI goes beyond predicting accurate force sets, as it can be used to identify the dynamics of domains,^{Reference Li, Yang, Zhang, Ye, Jesse, Kalinin and Vasudevan56} determine the structure of different topological structures,^{Reference Li, Nelson, Hsu, Damodaran, Li, Yadav, McCarter, Martin, Ramesh and Kalinin57} and search for new ME materials.^{Reference Xue, Balachandran, Yuan, Hu, Qian, Dougherty and Lookman58–Reference Balachandran, Emery, Gubernatis, Lookman, Wolverton and Zunger60} Yet, the requirement of big data sets is the “Achilles heel” of the machine-learning methodology and demands large amounts of accurate calculations on different scales, which will in turn stimulate the further development of existing methods.