Solutions for maximum coupling in multiferroic magnetoelectric composites by material design

Electrical control of magnetization offers an extra degree of freedom in materials possessing both electric and magnetic dipole moments. A stochastic optimization combined with homogenization is applied for the solution for maximum magnetoelectric (ME) coupling coefficient α of a laminar ME composite with the thickness and orientation of ferroelectric phase as design variables. Simulated annealing with a generalized Monte Carlo scheme is used for optimization problem. Optimal microstructure with single and poly-crystalline configurations that enhances the overall α is identified. It is found that juxtaposing a preferentially oriented ferroelectric material with a ferromagnetic ferrite into a composite would result in manifold increase in magnetoelectric coupling. The interface shear strains are found to be richly contributing to the ME coupling. The preferential orientation of the ferroelectric phase in the optimal ME composite laminate is demonstrated using the optimal pole figure analyses.


Electrical control of magnetization offers an extra degree of freedom in materials possessing both electric and magnetic dipole moments. A stochastic optimization combined with homogenization is
applied for the solution for maximum magnetoelectric (ME) coupling coefficient α of a laminar ME composite with the thickness and orientation of ferroelectric phase as design variables. Simulated annealing with a generalized Monte Carlo scheme is used for optimization problem. Optimal microstructure with single and poly-crystalline configurations that enhances the overall α is identified. It is found that juxtaposing a preferentially oriented ferroelectric material with a ferromagnetic ferrite into a composite would result in manifold increase in magnetoelectric coupling. The interface shear strains are found to be richly contributing to the ME coupling. The preferential orientation of the ferroelectric phase in the optimal ME composite laminate is demonstrated using the optimal pole figure analyses.
Magnetoelectric (ME) multiferroics have recently drawn increasing interest due to their potential applications in multifunctional devices such as nonvolatile memory elements, nano-electronics, etc. 1,2 . Natural multiferroic single-phase compounds are rare and their ME coupling responses are either weak or occurs at very low temperatures 1,3 . Nevertheless, substantial ME effect can be derived through fabricating composites of a ferroelectric (FE) and a ferromagnetic (FM) material in the form of thin film multilayered heterostructures, laminates, particulate or fiber-matrix composites etc 1,4 . Selecting materials possessing highly contrasting lattice parameters (thereby ensuring spontaneous immiscibility of phases) such as the perovskite-spinel system of electrostrictive and magnetostrictive phases in a laminar composite architecture could yield better ME coupling than single-phase materials or particulate composites 1 . Various efforts to improve the value of ME coupling coefficient α have been made by modifying preparation techniques of the samples, by the proper choice of materials or different structures, or by choosing different thickness of the samples 5,6 .
Theoretical characterization of equivalent material properties of complex material systems such as multiferroics is pursued as an alternative to experiments which are hampered by a host of factors affecting the sample such as demagnetization, debonding, microcracks, interdiffusion, depoling etc. [7][8][9][10] . The idea behind the development of ME multiferroic composites is to generate the desired magneto-electric effect as a strain induced product property 11,12 . The product property in a composite is defined as an effective property which is not present in either of the constituent phases 13 . Several analytical and computational models have been developed to quantify the product property of magnetoelectric coupling in ME multiferroic composites 8,[12][13][14][15][16] consequent to the strain-mediated two-phase model proposed by Harshe et al. 4,17 . Green's function technique was developed by Nan 18 for the solution of the constitutive equations and thereby compute the effective properties of ME composites. Micromechanics models too were developed as an alternative method to accomplish the same [19][20][21] . Theory of low-frequency magnetoelectric coupling in magnetostrictive-piezoelectric bilayers was developed by Bichurin et al. 15 . The ME equivalent circuit and corresponding ME coefficients have been derived for laminates by Dong et al. 22 . Phase-field models have been developed for predicting and interpreting both the equilibrium domain structures and domain switching kinetics by simulating the mesoscale microstructure of the magnetoelectric materials 23 . Nonetheless, most of the theoretical developments confine to two-phase materials which constitutes a subset of the multiferroics. The homogenization framework used in this paper is not specifically limited to two-phase composites alone, but provides a general platform to treat the equivalent (linear) static coupling interactions of all forms of multiferroics viz. single crystals, polycrystals, as well as all kinds of multiferroic composites irrespective of the crystallographic symmetry of the materials. The two-scale asymptotic analysis combined with a variational formulation of the underlying electrical, mechanical and magnetic fields would unveil their interaction in the microscopic scale.
Large ME coupling is indispensable for practical applications 1,8 . Recently multiferroic laminate composites of the transversely poled Poly(vinylidene fluoride-co-hexafluoropropylene-P(VDF-HFP) and the iron-based Metglas 24 was found to show ME voltage coupling α E33 as high as 320 V/cmOe. Large ME effects in a host of composites, including thin film heterostructures on nickel zinc ferrite, thick films of nickel ferrite and a piezoelectric deposited by aerosol deposition, in multiferroic composite bimorph structures are reported recently 9,[25][26][27] . Recent advances in manufacturing single crystal oriented multilayer films through heteroepitaxial growth and multilayering are found to enhance the magnetoelectric coupling 6 . If large ME coupling is achieved by upscaling the existing materials and geometries, it would cater to the economy of resources. We model the optimal material design ideal to accomplish this goal by combining the homogenization procedure with a stochastic optimization method.
In this paper, we are optimizing the microstructure (MS) of a laminar ME composite of two ferroic materials (ferromagnetic CoFe 2 O 4 (CFO) and ferroelectric BaTiO 3 (BTO)) that possessed contrasting lattice parameters 28 . To maximize the ME coupling, we search -by employing the method of simulated annealing based in the Metropolis algorithm-for the optimal MS configuration that could enhance this effect especially utilizing the inherent anisotropy of the component phases in composite multiferroics 29,30 . Recent advances in manufacturing oriented single crystal multilayer films through heteroepitaxial growth on top of a crystalline substrate offer grounds for such a study as this could be realized into potential device applications 1,6 . As a guide to the experimental realization of the ME composite, we explore the texture information through the optimal pole figure analysis.
Problems inherent to bulk ME composites such as destruction of piezoelectric charges by leakage current can be avoided in horizontal layered geometries 1,15 . Since the lattice mismatch between CFO and BTO is large (~0.52), it can either give rise to an heteroepitaxial strain in ME thinfilms 31,32 or enhanced interface stress in other horizontal layered media 4,28 . The ME coupling coefficient tensor α ij corresponds to induction of electric polarization by a magnetic field or of magnetization by an electric field and is designated as the linear ME effect 7 . α of a laminar ME composite is inextricably related to the thickness and orientation (or design variables) of FE phase. As α could not be expressed explicitly as function of the design variables, the optimization problem should be treated as a combinatorial optimization with a generalized Monte Carlo scheme 30 . A three-dimensional (3D) generic model is developed to compute the homogenized ME property tensor of a multiferroic possessing the least crystallographic symmetry (i.e., triclinic) to be interfaced with the optimization program 33 . See Supplementary Material for further details.

Optimization of magnetoelectric composite
The general homogenization method applied to magnetoelectric composite is based upon assumptions of periodic boundary conditions on the microstructure and the separation of the microstructure scale through asymptotic expansion. The asymptotic analysis of multiferroic ME material leads to the product property of the homogenized ME coupling 33 Here R and Q are characteristic electric and magnetic displacements. Γ, and Ψ are characteristic coupled functions satisfying a set of microscopic equations 33 . These characteristic functions appear in the microscopic displacement, electric and magnetic field perturbations, viz., u 1 (x, y), ϕ 1 (x, y) and ψ 1 (x, y) respectively, due to the microstructure heterogeneity. The spacial derivatives of the macroscopic displacement, electric and magnetic fields viz., u 0 (x, y), ϕ 0 (x, y) and ψ 0 (x, y) can bridge the macroscale and the microscale through the characteristic functions mentioned above. e pkl and κ pj H  are the piezoelectric coefficient and dielectric permittivity at constant strain  and magnetic field vector H respectively. =  i j k l , , , , 1, 2, 3 are the 3-dimensional coordinate indices, and δ ij are the Kronecker delta symbol. Here we assume Einstein convention on summation about repeated indices. (Here the ~ over Greek or Latin letters represents homogenized values). We consider the ME composite as a heterogeneous body obtained by the translation of microstructure of size Y. For the homogenized ME composite laminate, the physical properties do not depend on x, the global frame of reference. Instead, if the material is heterogeneous, the magneto-electro-mechanical properties effectively depend on x such that H H etc. The material properties κ H  , e, α, etc as well as the characteristic functions R, Q, Γ, and Ψ are Y-periodic functions depending on the microscopic coordinates y.
Since the homogenized or effective α  is the measure of the efficiency of the control of the electrical voltage on the magnetization state of the ME composite, its maximization would be the key objective sought in potential applications 34 . It is possible to obtain other relevant and experimentally measured ME parameters such as the effective ME voltage coefficient α of the composite. Here, E, P and κ 0 are the electric filed, polarization and permittivity of free space respectively. κ  is the effective relative permittivity of the composite.
Single crystal BTO-ceramic CFO. The homogenized α  33 and α  11 are set as the objective functions in this study. (Here the ~ over α represents homogenized values). α  33 and α  11 depends on the detailed configuration (anisotropy characterized by the orientation of the constituent phases and volume fraction or thickness of the constituent phases) of that system (here the ME composite laminate). Using the objective function α  k ( ) ij and defining configurations by a set of parameters {k}, it is straightforward with the Metropolis procedure to generate a set of configurations of a given optimization problem at some effective temperature (a control parameter) 30 . See Supplementary Material for further details. In single crystalline BTO-ceramic CFO laminate the optimization problem is to find vector {k} of design variables to minimize the objective function α  ij , i.e., α ij f f Two problems viz., the optimization of α  33 and α  11 are summarised in the above Eq. 2. Here, the orientation of the single crystal BTO layer of ME composite laminate possessing a volume fraction v f is characterized by the Euler angles (φ, θ, ψ). We first optimize a laminate with the objective of maximizing the α | |  33 . The minimization program depicted in Eq. (2) should manifest the real optimization problem where the objective functions α  33 or α  11 conventionally assume negative values 33 . Multiferroic brick elements possessing five degrees of freedom (DOF) viz., three for displacements, one each for the electric and magnetic potentials and a total of 13, 720 DOF are used for the computation of homogenized α  ij . See the Supplementary material for further details. The results of two trials of the normalized α  33 and α  11 are shown in Fig. 1(a,b). (Here the "step length" or the number of iterations in each step is set to be four.) Here the α  ij of the optimized system is compared against the bulk α  ijS values previously computed for isotropic composite. The inset to Fig. 1 depicts the MS with constituent volumes alone, while the optimal orientations of the FE phase are given in Table 1.
The solutions obtained show that the c-axis of the single crystal BTO lies parallel (θ = 0) to the y 3 axis of the MS reference frame (y 1 , y 2 , y 3 ) ( Table 1). The y 3 axis coincides with the spontaneous polarization direction. The optimal values for both α  11 and α  33 are found to be independent of the Euler angle φ which characterizes the ab-plane rotation of the BTO about the y 3 axis (that coincides with the c-axis of the crystal before the Euler rotation is taking place). In this rotated geometry of the composite, the ME coupling attains manifold enhancements from the corresponding bulk composite values. Here we see the α  11 attains more than 7 times the bulk value while α  33 displays nearly 6-fold enhancement due to rotation of the FE phase. The optimal volume fraction v f = 0.69 of single crystal BTO is in agreement with the experiments (Table 1) in BTO-CFO thinfilm 35,36 . Heteroepitaxial films of CoFe 2 O 4 and BaTiO 3 grown on (001)-SrTiO 3 substrate exhibits the optimal c-domain orientation of BaTiO 3 layer 35 . i.e., the c-axis of the BTO unit cell is parallel to the substrate normal or the Euler angle θ subtended by the c-axis would be θ = 0 as in our study. Yang et al. 37 found that the α of the composite in a configuration with PFN-PT oriented with its [001] direction out of plane of the laminate is ~4 times as that with its [001] lying along the plane of the lamina. Our optimal θ value (θ = 0) in both instances of optimization (Table 1) shows that BTO oriented with its [001] axis out of plane of the lamina at maximum α. ME effect in a composite is a coupled electrical and magnetic phenomena via elastic interaction 8 . i.e., for the magnetoelectric effect, when a magnetic field is applied to a composite, the magnetic phase changes its shape magnetostrictively. The strain is then passed along the piezoelectric phase, resulting in an electric polarization. Thus the strain-mediated ME effect in composites is extrinsic, depending on the composite microstructure and coupling interaction across magnetic-piezoelectric interface 38 . The property ME coupling α which is absent in the constituent phases thus originates in the composite. In order to explore the the enhancement of ME coupling at the optimal orientation, several factors were analyzed including interface strain field of the laminate. The laminate plane is essentially the xy-plane of the composite and the lamination is along the z-axis direction. Hence the interface is ideally spread along the xy-plane. The local distortion caused by the rotation of the BTO layer would impart additional strains at the interface besides the ones existing due to the lattice mismatch. We observe an additional component of elastic stiffness viz., C 16 = 9.1 × 10 9 N/m 2 in the optimal composite in such a way that C 16 = −C 26 as shown in the elastic constant (C μν ) matrix given below; This property is non-existent in either of the BTO and CFO phases which were individually treated to be belonging to the 4 mm tetragonal class. Nevertheless, the C μν matrix of the composite reveals that the overall symmetry of the optimal composite changed to one of the tetragonal classes m 4, 4, 4/ 39 . This symmetry change might be a consequence of the strain effect at the interface and this structural transition can bring about an enhancement of ME coupling α as it is seen here 40 .
The interface strain field developed as a result of the composite geometry and the rotation of the ferroelectric phase plays a crucial role in the magnetostriction in the ferromagnetic phase 41 . To further explore the interfacial strain effect on the ME coupling α we study the effective compliance  S 66 which is equivalent to  σ ≡   S / 1212 12 12 in tensor form of the ME laminate at various orientations. Here 12  and σ 12 are the in-plane shear strain and stress respectively. Figure 2(a) depicts the effective compliance  S 66 of the ME laminate and the corresponding out-of-plane ME coupling at various Euler angles θ which corresponds to the orientation of the c-axis of the BaTiO 3 unit cell. It reveals that at optimal angle θ the  S 66 attains the maximum (and so is the ME coupling α | |  33 ) and decreases monotonically as we increase θ. The decreasing homogenized piezoelectric constant  d 15 shown in Fig. 2(b) with θ reinforces the role of interfacial shear strain since ≡   d d 15 113 in tensor form. (d 113 is essentially the piezostrain  13 per unit electric field E 1 acting along the x-direction of the crystal.). Thus the shear interfacial strains which get bigger at the optimal point in the solution space obviously contribute to the enhancement of ME coupling. Ceramic BTO-CFO composite. Next is the optimization problem of ceramic BTO-ceramic CFO laminate to find the vector k of design variables to minimize the objective function α  ij (i.e., either α  33 or α  11 ), i.e., α σ μ σ μ σ μ σ σ σ  Table 2 carries the optimal values obtained for three (σ χ , μ χ ) pairs. Figure 3(a,b) show the results of this optimization. The MSs shown inset provide only the volume fraction of the constituent BTO and CFO phases of the initial guesses and the solution. In bulk BTO-CFO type eutectic composite samples, the volume fraction at which the maximum ME coupling is realized comes below the value obtained in our study (Table 2), except the one reported in ref. 42 where it is congruent with our results. For instance, in BaTiO 3 -Ni(Co, Mn)Fe 2 O 4 , the volume fraction of BTO was obtained as v f = 0.4 (ref. 43 ) while in another BTO-CFO eutectic composite it was found to be v f = 0.6 (ref. 44 ). Analytical results on α show an increase in the value until v f of BTO turns 50% and then remains constant 4 . The ME coupling relies on many factors, such as the stoichiometry, lack of texture, disparity of the character of the composite due to chemical modifications, porosity etc 2,7 . The stereographic projections, the (hkl) pole figures of BTO essentially plots the poles projected by the [hkl] directions on the reference sphere constructed with the local coordinate system as its radii 45 . Figure 4 gives the pole figures of the planes {100}, {001} and {001} of the grains of BTO used in the starting configuration of the α  33 optimization. The orientations of BTO were generated based on the values μ σ μ σ μ σ = φ φ θ θ ψ ψ     ( , ) ( , ) ( , ) (0, 1) (0, 1) (0, 1) with 1183 grains for the starting point of the optimization algorithm. Figure 4 represents some texture as we can see  Here the number of grains in MS become 1690 at optimum. This signifies the fact that majority of [001] axes, the spontaneous polarization direction of BTO, are oriented parallel to the y 3 direction of the lamina. This implies that in order to achieve the enhanced out of plane coupling, one needs to have a texture which favours the alignment of maximum number of grains (or rather its spontaneous polarizations) towards the y 3 axis. Moreover, it is seen that majority of poles in the (100) and (010) pole figures in Fig. 5 are concentrated along the circumference. This indicates that in the optimum configuration where the ME coupling α  33 is maximized, the ab-planes of the majority of BTO crystallites lies parallel to the laminar plane while obviously the c-axes lie out of plane of the laminar plane. i.e., majority of the BTO domains constituting the FE phase are c-domains where the c-axes of the unit cells lie parallel to the laminate normal. Nanocomposite ME films were shown to exhibit strong ME coupling consequent to preferential orientation of the FE (PZT) layer favoring the formation of c-domains is shown in the present study 10,46 . The pole plots for the optimization of α  11 are shown in Figs 6 and 7. The initial guesses used for pole figures (Fig. 6) are μ π σ μ π σ μ π σ Here we can conclude from the third picture of Fig. 7 as the texture is [001]-oriented. Thus, if we introduce proper texture you can get enhancements of ME coupling to the tune of α ∼  5 11 times and α ∼  2 33 times that of the bulk composite (Table 2). ME voltage coefficient α E . Another parameter of importance is the overall ME voltage coefficient of the laminate. Generally, α  E is determined when the sample is subjected to a bias field H and an AC field δH by measuring the electric field δE 12 . The strong ME coupling is expected in a layered structure primarily due to low leakage currents and ease of poling to align the electric dipoles and thereby strengthen the piezoelectric effect of the ferroelectric phase of the composite 27 . For the out-of-plane ME voltage coefficient α α κ =    / E33 33 33 at the optimal configuration 8 , we get a value of α ≈  193 E33 V/cmOe for single crystal BTO-ceramic CFO laminate. Here, α  33 = −1.546 × 10 −9 Ns/VC and κ  33 = −0.637 × 10 −11 F/m at the optimum obtained in the study. While for ceramic BTO-CFO laminate the optimal value of α  E33 obtained is α ≈  143 E33 V/cmOe. (Here the α  33 = −0.574 × 10 −9 Ns/VC and κ  33 = −0.319 × 10 −11 F/m at the optimum. Epitaxially grown CoFe 2 O 4 -BaTiO 3 heterostructure on the 001-SrTiO 3 substrate via pulsed laser deposition displayed an out-of-plane ME coefficient of 104 mV/cmOe, which is orders of magnitude lower than the present value, however 35 . Ren et al. 47 obtained value up to 2.54 V/cmOe at 160 kHz, for a particulate composite of BTO-CFO synthesized by a one-pot process.
The experimental values are generally lower in the case of particulate composites compared to laminates, due to elimination of leakage problem 10 . We assume a perfect interface between the laminae which is not the ideal case in experiments 8 . In the case of heteroepitaxial thin film configurations such as in ref. 35 , the clamping due to stiff substrate would essentially restrict the ME coupling response 10 . Besides this, the discrepancy between the theory and the experiment is due to a host of other factors such as; defects of the sample chiefly due to interphase diffusion of the constituent atoms and the consequent deterioration of piezoelectricity and magnetostriction of the constituent phases resulted in poor strain transfer across phases, porosity, large thermal expansion mismatch and the subsequent formation of microcracks, demagnetization, leakage currents etc. 7,8 .
Some multiferroic laminate composites of the transversely poled Poly(vinylidene fluoride-cohexafluoropropylene-P(VDF-HFP) and the iron-based Metglas fabricated by Jin et al. 24 shows α E33 values as high as 320 V/cmOe. (001)-oriented PMN-PT crystal laminated with (211)-grain oriented Terfenol-D was found to yield ME voltage coefficient α E33 = 30.8 V/cmOe at a resonance frequency of ~78 kHz 25 . The ME response of a five layer Metglas/Terfenol-D/PMN-PZT single crystal/Terfenol-D/Metglas laminate at 1 kHz was found to be 5 V/cmOe 48 . These values fall in the range of and resembles the values obtained for the ME laminates obtained in the present study. Strong ME effects in a variety of ferrite based composites, including thin film piezoelectrics on nickel zinc ferrite, dense films of nickel ferrite and a piezoelectric deposited by aerosol deposition, in multiferroic magnetostrictive/piezoelectric composite bimorph structures are reported recently 9,[25][26][27] . Recent advances in manufacturing single crystal oriented multilayer films through heteroepitaxial growth and multilayering is found to enhance the magnetoelectric coupling 6 . In epitaxial growth of thin films, an oriented crystalline film is deposited onto a crystalline substrate, while in this study a polycrystalline ferromagnetic layer bonded to a crystalline ferroelectric layer. A recent report 49 , however presents an approach namely combinatorial substrate epitaxy (CSE) relies on the creation of polycrystalline substrates and provide access to a wide range of materials and numerous orientations in a single pass. Studies in this direction would plausibly promote experiments aimed at realizing optimal multiferroic material geometries where the ME coupling is maximized as is manifested in the present study.
In summary, a computational framework for the material design and pole figure analysis of ME composites identified a range of grain orientation/distribution of the FE phase in ME composite wherein the ME coupling can be enhanced manifold is identified. The single crystal phase of BTO is found to be ideal for the maximum ME coupling configuration compared to the polycrystalline phase. A systematic increase in interfacial shear strains is identified with the alignment of the c-axis of the BaTiO 3 single crystal phase of the composite along the global z-axis. Yet, the data pertaining to the ceramic BTO-CFO demonstrates unique solutions for the extreme coupling comparable to the single crystalline phase which can provide the manufacturers with more degrees of freedom. The insight obtained from the optimization has the potential to advance the design and upscaling of complex ME composite configurations with superior coupling. Further studies on laminates consisting of relaxor ferroelectrics, which possesses larger piezoelectricity, and ferromagnets with high piezomagnetic coupling such as ferrites, transition metals and rare-earth alloys could inaugurate new possibilities in technological applications demanding larger ME coupling.
See Supplementary material for the details on the multiferroic homogenization, convergence analysis for representative volume element and the optimization algorithm.