Homogenization method for microscopic characterization of the composite magnetoelectric multiferroics

Tuning of magnetization or electrical polarization using external fields other than their corresponding conjugate fields (i.e., magnetic field for the former or electric field for the latter response) attracts renewed interest due to its potential for applications. The magnetoelectric effect in multiferroic 1–3 composite composed of alternating magnetic and ferroelectric layers operating in linear regime consequent to external biasing fields is simulated and analysed theoretically. Two-scale homogenization procedure to arrive at the equilibrium overall physical properties of magnetoelectric multiferroic composite is formulated using variational analysis. This procedure is extended to quantify the underlying local (microscopic) electric, magnetic and elastic fields and thereby compute local distribution of stresses and strains, electrical and magnetic potentials, the electric and magnetic fields as well as the equivalent von Mises stresses. The computational model is implemented by modifying the software POSTMAT (material postprocessing). Computed local stress/strain profiles and the von Mises stresses consequent to biasing electrical and magnetic fields provide insightful information related to the magnetostriction and the ensuing electrical and magnetic polarization. Average polarization and magnetization against magnetic and electric fields respectively are computed and found to be in reasonable limits of the experimental results on similar composite systems. The homogenization model covers multiferroics and its composites regardless of crystallographic symmetry (with the caveat of assuming an ideal and semi-coherent interface connecting the constituent phases) and offer computational efficiency besides unveiling the nature of the underlying microscopic field characteristics.


I. INTRODUCTION
The materials classified as magnetoelectric multiferroics possess both the magnetic and ferroelectric orders 1 . More than having them make it to exploit the functionalities of both the orders, a coupling between the ferromagnetic and FE states enable appearance of novel properties not present in either of the states 2 . Simultaneous occurrence of magnetism and ferroelectricity in multiferroic materials is constrained by the conflicting physical requirements of partially filled d orbitals for magnetism and empty d orbitals for ferroelectricity to manifest 3,4 . These restrictions results in rare occurrence of intrinsic multiferroics.
Nevertheless, substantial ME effect can be derived through fabricating composites of a ferroelectric (FE) and a ferromagnetic (FM) material in the form composites 5 . Crystallites of the two phases are assumed to be in good mechanical contact if the two phases are polycrystalline. The FE phase is poled near the ferroelectric Curie temperature (T C ) in a strong electric field to make the composite piezoelectric while the magnetic poling of the FM phase is accomplished in a similar way, by annealing the composite in a magnetic field near the Néel temperature (T N ). When an electric field is applied to this composite the FE grains elongate parallel to the electric field. The change in shape of the ferroelectric grains causes the ferromagnetic grains to deform, resulting in a change in magnetization 6 . The connectivity too play an important role in property development in multiphase solids where the connection patterns can change some physical properties 7 .
Inclusion of FM nanoparticles into a FE matrix to form bulk particulate composites, where the interfacial strain is transmitted through the grain boundaries constitute one of this routes 8 . Another important way is stacking of FM and FE sheets into composite magnetoelectric laminate where the strain is transmitted between phases relying on the bonding between the layers 9 . Since the voltage control of magnetism in magnetic metals are hampered by the short screening lengths, researchers began to explore ultra thinfilms 10 and magnetic metals/ferroelectric oxide bilayers 11,12 . Electrostatic doping can alter the FE-FM interface and the magnetoelectric effect can penetrate deeper than the screening length into the metallic La 0.67 Sr 0.33 MnO 3 (LSMO) in a LSMO-BTO heterostructure 13 . Single crystal ferroelectric substrates are being used as FE-FM systems for utilizing the anisotropy and symmetry to its full potential in planar geometry in multilayer thinfilms. The magnetic field induced electric polarization (direct magnetoelectric coupling) is so vital in device applications that the magnetic field-produced strain in a magnetic phase transferring into a piezoelectric phase and inducing a charge output, is the main working mechanism of magnetic sensing 14 , and also perceived as a criterion for magnetic detection. Burton and Tsymbal show that magnetic reconstruction induced by switching of electric polarization can be facilitated by the doping of divalent cation near the magnetic phase transition 15 . Uetsuji et al., have studied the coupling mechanical response consequent to the application of a magnetic field in magnetoelectric laminated composite 16 .

A. Constitutive equations and homogenization
The magnetoelectric effect is described by a term in the thermodynamic potential Φ me that is linear both in the magnetic (H) and electric (E) field, such that where the magnetoelectric coupling coefficient α ik is an unsymmetrical tensor 17 . When H = 0, the y 2 y 1 y 3 Figure S1. The successive angles of anti-clockwise rotations subjected to the crystallites (grains) on the ferroelectric phase of the ME composite as defined by the Euler angles φ • , θ • and ψ • . Here (y1, y2, y3) ≡ y are the local and (x, y, z) ≡ y are the crystallographic coordinate systems.   electric field generates a magnetization and when E = 0, the magnetic field generates an electrical polarization A two-scale asymptotic homogenization analysis combined with a variational formulation is developed for determination of the equivalent material properties of a periodic multiferroic magnetoelectric composite. Local and average (global) electrical, magnetic, and mechanical constitutive behaviour are computed. The model takes into account the ferroelectric and ferromagnetic, and the composite magnetoelectric phases by treating the constitutive behaviour. For a linear magneto-electro-elastic solid, constitutive equations are governed by the electrical, mechanical and magnetic fields. The model quantifies the local electrical and magnetic potential, displacements, electrical and magnetic fields, stress and strain fields, magnetization (through magnetic flux density) and polarization (through electric displacement) and von-Mises stress, besides the effective magneto-electro-mechanical properties. The general homogenization method applied to multiferroics has no limitations regarding volume fraction or shape of the constituents involved and is based upon assumptions of periodicity of the microstructure and the separation of the microstructure scale through a proper asymptotic expansion 18 . The mathematical theory of homogenization accommodates the phase interaction in characterising both the macro-and micro-mechanical behaviours of the composite material. i.e., the method permits the introduction of different field equations in a microscopic scale to each constituent of a composite while following the representative volume element (RVE) notion. This study considers multiferroic materials that respond linearly to changes in the electric field, electric displacement, mechanical stress, strain as well as magnetic field. Let Ω be a fixed domain in xspace. We consider an auxiliary y-space divided into parallelepiped periods Y. For a linear anisotropic magnetoelectric material, generated through the periodic repetition of a unit cell representing the smallest sample of heterogeneity of the material domain Ω, the governng equations are given below; force equilibrium equation, electrical (magnetic) field-electrical (magnetic) potential relations and the quasistatic steady-state Maxwell's equations for electromagnetic phenomena, In a multiferroic solid, the mechanical, electrical and magnetic variables are related by constitutive relations. For small deformations, the linear constitutive laws of multiferroics in the absence of heat flux are given by Here σ, , D, B are stress, strain, displacement, body force, mass density, electric displacement vector, and magnetic flux density respectively. C EH , e, e M , κ H and µ E are stiffness, strain to (electric, magnetic) field coupling constants (or piezo-electric and -magnetic coefficients), permittivity (dielectric) and (magnetic) permeability respectively. Considering the standard homogenization procedure, the material functions C EH , e, e M , κ H and µ E , are considered to be Y-periodic functions in the unit cell Here the two scales x and y are spatial variables where x is a macroscopic quantity and y is a microscopic one 19 . The functions involved in this expansion are assumed to be dependent on these two variables, where one (i.e., x) describing the "global" or average response of the structure and the other (i.e., y) describing the "local" or microstructural behaviour. Here the two variables x and x/ε take into account the two scales of the homogenization; the x variable is the macroscopic variable, whereas the x/ε variable takes into account the microscopic geometry.
The energy functional of the system consists of contributions from magnetic and electric parts and the magnetoelectric coupling component, which describe the effective interaction between order parameters of each phase. Let the surface Γ of a magnetoelectro-elastic body Ω be subjected to prescribed surface traction t k and surface charge per unit area σ and normal magnetic flux density B n . The energy functional 20 for a magnetoelectric multiferroic can be written as, Here ρ ε (≡ ρ(x, y)) is the density and b(≡ b(x, y)) is the body force. The surface traction t k , surface charge density σ and normal magnetic flux density B n are treated to be independent of the spatial scale ε. The spatial derivatives of displacement and potentials are taken and substituted into the Eqs. (S5) and (S6) and will eventually be substituted in the energy functional G. Applying calculus of variations, passing through the limit ε → 0 and advance using asymptotic analysis one obtains the effective magneto-electro-elastic moduli, viz., the homogenized elastic stiffnesses C EH ijkl , piezoelectric coefficients e ijk , piezomagnetic coefficients e M ijk , dielectric permittivities κ H ij , magnetic permeabilities µ E and the magnetoelectric coupling coefficients α ij of the magnetoelectric multiferroic. In other words one could get a good approximation of the macroscopic behaviour of such a heterogeneous material by letting the parameter ε, which describes the fineness of the microscopic structure, tend to zero (ε → 0) in the equations describing phenomena. And for a homogeneous material the physical properties does not depend on x. The detailed theoretical analysis is given elsewhere 18 . (In all the expressions of this paper, it may be noticed to discern the difference in notations of asymptotic scale factor ε and the mechanical strain . All the microscopic quantities will carry the notation of ε to identify their local character). The local strain ε ij (x), electric field E ε j (x) and the magnetic field H ε j (x) too can be obtained in a similar way once the homogenized macroscopic problem is solved.
The microscopic stress σ ε ij (x), electrical displacement D ε i (x) and magnetic flux densities B ε i (x) at each point of the domain can be computed using the constitutive equations (S8)-(S10), and the field equations (S5) and (S6) as The flux represents the average normal component of the field times the surface area and in here the same definition is applicable about magnetic flux density B. The electrical displacement D would be essentially a measure of polarization P consequent to an applied electric field E in a general dielectric material. The equivalence of the average stress and homogenised stress can easily be expressed by applying the average operator denoting ( 1 |Y | Y dY ), where Y is the domain size. Hence the average fields are computed to be These macroscopic equations (i.e. they do not contain y) can be computed once the homogenized solution for u 0 , ϕ 0 and ψ 0 and that of the corresponding fields viz., and ∂ψ 0 (x) ∂xj are prescribed. This postulate is equally applicable for the case with the microscopic fields depicted in equations Eqs. (2)-(4) from the main text and Eqs. (S12)-(S14) for the microscopic stress, electrical displacement and flux densities given in equations (S15)-(S17).

III. NUMERICAL IMPLEMENTATION
In order to solve the microscopic system of equations resulting form homogenization we developed a finite element (FEM) formulation and the details are given elsewhere 18 . A three-dimensional (3D) multiferroic finite element is conceived with five degrees of freedom (DOF)-three DOFs for spatial displacements and one each for electric and magnetic potentials. Eight-noded isoparametric elements with 2 × 2 × 2 Gauss-point integration are used obtain solutions. Altogether there were nine microscopic equations that should be solved for as much number of unknowns namely the characteristic functions χ mn i , η mn , λ mn , R m , Φ m i , Θ m , Q m , Ψ m and Γ m k , where the indices m, n = 1, 2, 3 18 . The problem is reduced to standard variational FEM, after the usual approximations of finite element formulation and can be expressed concisely as where K is the global stiffness matrix, u is the vector of unknown functions and f is the load vector. Magnetoelectric multiferroic material, in general, can be considered as an aggregate of single crystalline crystallites/grains and hence the system altogether would be polycrystalline nature. The unit cell or the representative volume element (RVE) conceived in this work is a volume containing a sufficiently large number of crystallites or grains that its properties can be considered as equivalent to that of the macroscopic sample. The electric polarization P as well as the magnetization M interspersed inside the crystallites could be mapped using using some coordinate system. In a sense the underlying crystal orientation can encompass the orientations of P or M. Thus we introduce the Euler angles (ϕ, θ, ψ) to quantify the crystal orientations of a multiferroic polycrystal (Fig. S1), as the crystallites in an as-grown sample are randomly oriented in the lattice space and hence require three angles to describe its orientation with reference to a fixed coordinate system. Here we use the so called x-convention, where the first and third rotation is through the y-axis (here it is y 2 -axis)and the second rotation is through the intermediate x-axis (here it is y 1 -axis). Thus all the physical quantities λ ijklmn... (y ) expressed in a crystallographic coordinate system y would be coordinate-transformed to the local coordinate system y according to the following scheme erty data entered into the homogenization program are obtained with respect to the crystallographic coordinates.) Here e µν are the Euler transformation matrices 21 .
In the FEM simulation of the homogenization, a multiferroic crystallite is represented by a finite element in the unit cell. Thus we have a polycrystalline unite cell having as much number of crystallites as the number of finite elements by which it is discretized. As-grown FE (or for that matter FM) polycrystal, often ends up in a near complete compensation of polarization (or magnetization) and the material consequently exhibit very small, if any, electric (or magnetic) effect until they are poled by the application of an electric (magnetic) field. The orientation distribution of the crystallites (grains) in such a polycrystalline material would be uniform with a standard deviation σ → ∞ before poling (application of electric/magnetic field) and that after poling would best be represented by a distribution function with σ → 0. Thus, any pragmatic configuration of orientation distribution of grains in multiferroic material would fit in a Gaussian distribution defined by the probability distribution function where µ and σ are the the mean and the standard deviation of the angles α (which stands for the Euler angles (φ, θ, ϑ)). The convergence of magnetoelectric properties with unit cell size allows us to determine the simulation-space independent, equivalent magnetoelectric properties of the composite. Convergence analyses, on magnetoelectric composites reveal that accuracy one derives from descretizing the unit cell (in other words sampling of the unit cells with more number of grains or less number) is minimal above 1000 elements (grains) 22 . Consequently, we kept unit cells' sizes greater than 1000 finite elements in this study.

IV. RESULTS
The elemental averages of local fields computed from the nodal values for each element of the FEM of the unit cell are shown in Fig. S2. It would be interesting to see the strain curves of this analysis as it is shown in Fig. S3. The strains are a consequence to the response of the material to stress and here it is caused by the external electric field E. The anomaly at the boundary between the phases is obvious here as was seen in the stress analysis. The tensile component along the y 1 -axis ( ε 11 ) acquires positive value in the FM fibre phase of the composite while it is seen oscillating from tensile to compressive strain in the FE matrix.
The local and global fields using a honeycomb lattice unit cell of the ME composite is studied to verify the robustness of the POSTMAT developed in this work. The constitution of the composite material is not changed and the homogenization is run first to characterise the equivalent magneto-electroelastic properties and the macroscopic strain, electric and magnetic potentials and their derivatives. We have kept the stoichiometry of the composite at 0.65BaTiO 3 -0.35CoFe 2 O 4 first as in the above study. (Here the local distribution diagrams are shown slanted to enhance the visibility of the picture.) Here only the local distributions are drawn omitting symmetric and marginal ones. For instance, σ ε 11 ≡ σ ε 22 and hence σ ε 22 is not shown in the picture.
We have conducted a numerical experiment to simulate the magnetoelectric composite system described in Kuo and Bhattacharya 23 . To mimic their model we first take an ME composite with BaTiO 3 fibres embedded in CoFe 2 O 4 matrix with fibre volume v f = 0.2 where an external average magnetic field H x = 1 A/m is applied to the composite along the in-plane x-direction. The resulting local displacement (u ε ), electric potential (ϕ ε ) and magnetic potential (ψ ε ) contours are plotted and is shown in Fig. S7(a-c). The next column of subplots i.e., Fig. S7(d-f) pertains to the case where the composite is composed of CoFe 2 O 4 fibres (with CFO v f = 0.2) embedded in BaTiO 3 matrix. Here in both cases a transverse 2D section of the profiles at the mid-point of y 3 -axis of the composite microstructure is given. The profiles are in good accord with the ones shown in Kuo and Bhattacharya 23 while admitting that we have used single crystal BaTiO 3 as matrix instead of the ceramic matrix used by them. It is noticed that their model treats only antiplane shear with in-plane electromagnetic fields and contrasts with ours where we have admitted all the shear components and electromagnetic fields.