Simulating the non-monotonic strain response of nanoporous multiferroic composites under electric field control

Multiferroic composites have attracted a considerable amount of research in recent years as functional materials that exhibit the magnetoelectric (ME) effect. These materials transfer the strain induced by an external electric field in the ferroelectric component to the ferromagnetic component through elastic coupling. The application of multiferroic materials can drastically reduce power consumption in small-scale devices, such as sensors, switches, and data storage elements, compared to standard electromagnetic devices.^1,2

However, multiferroic devices composed of dissimilar materials arranged in alternating geometries^3–10 are often limited by poor mechanical coupling, thermal expansion coefficient mismatch, or irregular interfacial bonding. In many instances, these issues appear during synthesis and can all contribute to a weakened ME effect.^11–13 In recent years, research on nanoscale multiferroic heterostructures,^14,15 organic multiferroics,¹⁶ or thin films¹⁷ has demonstrated the feasibility of nanosynthesis techniques. As well, theory and simulation are being used to establish the suitability of downscaling different multiferroic designs as well as assessing their potential response.^18,19

In this work, we present finite element simulations of a multiferroic structure generated by infiltrating a ferroelectric phase, lead zirconate titanate, onward referred to as “PZT” (PbZr_xTi_1-xO₃, where, in this work, x = 0.52) via atomic layer deposition (ALD) into a nanoporous ferromagnetic template made of cobalt ferrite (CoFe₂O₄, referred to as “CFO”). Such structures have been seen to yield an enhanced ME response and be free of many of the issues limiting other multiferroic devices.^20,21 We focus on the dependence of the deformation behavior of the device on the thickness of the ALD film coating the interior of the pores. In particular, we model the non-monotonic strain response of the structure as a function of PZT thickness and analyze the physical causes behind such response.

We begin with a brief description of the synthesis procedure and magnetic measurements of the nanostructure. Nanostructured PZT/CFO composites were created in two parts: first making the porous CFO architecture, then conformally coating the inside of the pores using ALD. First, porous CFO structures were synthesized by block copolymer templating of sol gel films, as reported elsewhere.^22,23 In brief, metal precursors were dissolved and, over time, condensation reactions form a metal oxide polymer known as a “gel.” An amphiphilic di-block copolymer is added, which forms micelles in solution. As the solution is deposited onto the conductive silicon substrate, micelles self-assemble into ordered, periodic structures. Once pyrolyzed, the polymer is removed, leaving behind a porous, metal oxide framework. ALD was then used to conformally coat the inside of the pores with PZT, the details for which are discussed elsewhere.²⁴ Since ALD is a gas phase deposition process, the ALD precursors uniformly bind to available surface sites throughout the porous network. The thickness of the PZT coating is precisely controlled by the ALD process, and the ALD layers were then crystallized through rapid thermal annealing, or RTA, at 700 °C for 60 s under oxygen. Figure 1 shows the bare CFO template prior to ALD. To observe the magnetoelectric coupling of the PZT/CFO composite structure, samples were electrically poled ex situ with the electric field being applied perpendicular to the sample. Samples were cut down to be approximately 5 × 5 mm² and then sandwiched between aluminum electrodes with a 13 μm thick dielectric spacer above the film side to prevent shorting. The electric fields utilized in this work ranged between 0 and 1.42 MV/m. Once poled, high resolution diffraction measurements were performed at Stanford Synchrotron Light Source (SSRL), experimental station 7–2, with an x-ray energy of 0.9919 and 1.0332 Å. The d-spacing of the composites was determined by peak fitting the CFO (311) diffraction peaks to a Voigt function (through PeakFit v4.11²⁵). Measurements shown in this work were performed out of the plane of the sample, showing out-of-plane tension.

For their part, finite element (FE) simulations of a nanoporous network of CFO (consistent with the structure shown in Fig. 1) internally coated with PZT were carried out using COMSOL Multiphysics.²⁷ The CFO network functions as the ferromagnetic substructure while the PZT filling provides the ferroelectric response. Designed in this way, the simulated configuration represents a piezoelectric template with an infiltrated ferromagnetic component and not a multiferroic structure per se. However, the computational model acts as a viable surrogate for the actual multiferroic material in that it captures its physical response as if there was a true magnetic coupling. The pores of the CFO substructure are arranged into a close-packed three dimensional lattice (i.e., a face-centered cubic lattice, or FCC). The pores are coated with a PZT layer of varying thickness (3, 6, and 10 nm). The CFO and the PZT are assumed to be ideally bonded (no gaps or cracks) with the displacements being continuous across their interface (known as Dirichlet boundary condition). The size of the ‘unit’ FCC cell for the CFO structure was $≈$42 nm, while the pore diameters were in the 15–20-nm range. The radius of the ligaments connecting the pores (only along the first nearest-neighbor distance of the underlying FCC cell) was 8 nm. Figure 2(a) shows images of the unfilled CFO substructure, and the coated structures with 3, 6, and 10 nm-thick PZT [Figs. 2(b)–2(d), respectively]. With 10-nm PZT layers, the ligaments become fully infiltrated and the pores are no longer interconnected. The porosity of the entire structure for each ALD-coating thickness is 15.3, 6.6, and 0.03%, respectively.

Application of an external electric field, $E$⁠, results in a piezoelectric coupling with the ferroelectric component of the device, i.e., the PZT coating. The PZT deforms, transferring part of the deformation to the ferromagnetic component CFO, which is what would then trigger a magnetic response. The governing equations for the piezoelectric effect can be expressed in the stress-charge form as

or in the strain-charge form

where T and S are the stress and strain tensors, respectively; E is the electric field (vector); D is the electric displacement vector; $ℂE=SE−1$ is the elasticity matrix; $e$ and $d$ are the converse and direct piezoelectric coefficient matrices; $ε0$ is the permittivity in free space; and $εS$ and $εT$ are the material permittivity matrices at constant strain and stress, respectively. $ℂE$ (⁠$SE$⁠), $e$ (⁠$d$⁠), and $εS$ (⁠$εT$⁠) are tensors of rank four, three, and two, respectively. These tensors, however, are highly symmetric for physical reasons, and they can generally be represented as sparse matrices using Voigt notation, which makes their algebraic treatment more convenient. COMSOL includes a standard piezoelectric constitutive model, such as this one in its Multiphysics release.²⁸ The values used in the present work for PZT are listed in the supplementary material document attached to this paper. CFO is simply treated as a compatible material with a cubic structure. Its elastic properties are given in the supplementary material as well.

Closure of the system of Eqs. (3) and (4) is achieved by adding the equilibrium conditions $∇·D=0, ∇×E=0, ∇·T=0$⁠, and the definition of the strain tensor from the displacement vector u: $S=1/2(∇u+u∇)$⁠. The system is assumed to be constrained by the device geometry along the x and y directions, and thus, no displacements are allowed on those boundaries, i.e., $u=(0,0,uz)$⁠. Along the z direction, displacements are uniform across the entire top surface, with u_z equal to the displacement of the center of mass of that surface. By way of example, the mesh for the 3.0-nm-thick PZT layer consists of over 100 000 elements and is shown in Fig. 2(b).

Next, we simulate the strain response of the structures shown in Fig. 2 as a function of the applied external electric field assuming perfectly linear piezoelectric behavior and ideal CFO/PZT interface bonding. Demonstrative examples are given in Fig. 3 for the displacement [Fig. 3(a) and stress Fig. 3(b)] fields for the 3-nm PZT layer case when $Ez=1.2$ MV m^–1. The results as a function of E_z are shown in Fig. 4, where the experimental measurements have also been included. While based on these results, it may appear that the relationship between PZT-coating thickness and applied voltage is monotonically decreasing, next we map the response of thicknesses from 1 nm to full infiltration (zero porosity, around 11 nm) to get a clearer picture of the dependence on PZT layer thickness. A surface plot showing the calculations of the strain as a function of these two variables is shown in Fig. 5. As the data show, the peak piezoelectric effect is reached for a coating thickness of 3 nm, increasing sharply from zero, and decreasing smoothly thereafter. To understand this behavior and also to explain the differences with the experimental data, next we explore the effect of non-ideal conditions on device response.

We start by considering the effect of nonlinear piezoelectricity on the current design response. In general, the polarization of a crystal under the application of an electric field or a mechanical stress consists of intrinsic and extrinsic contributions.^29,30 The intrinsic contribution is associated with the atomic structure of the crystal lattice, and it is strictly speaking the only one that can be considered linearly coupled to the electric field. The extrinsic response results principally from $90°$ domain wall movements,³⁰ which naturally exist in ferroelectric materials, like PZT.³¹ Poling of PZT prior to electrical stimulation enhances this nonlinear effect by aligning like-domains and enhancing the polarization response. Indeed, domain wall movement increases with the intensity of the electric field and, thus, the extrinsic part of the polarization is seen to couple more strongly to E. The standard model that captures this nonlinearity is the Rayleigh law,^32–34 which changes the axial component of $d$⁠, d₃₃ (the only one subjected to extrinsic polarization) to

where $dint=593$ fC N^–1 (cf. supplementary material) represents the intrinsic part of the coefficients and α is the Rayleigh constant. For tetragonal PZT, $α=6.02×10−18$ m² C N^–2 (Ref. 35). Consistent with ex situ poling of our multiferroic structure, we carry out computer simulations with the updated Rayleigh model. The results are shown in Fig. 6, with the calculated data points showing a larger deviation of the device strain with PZT-coating thickness as the electric field increases. As can be seen, this results in a better agreement with the experimental results.

While the magnitude of the results is different between the results in Figs. 4 and 6, the S_zz-E_z trend remains. Next, we carry out extra analysis to understand the reasons behind this. As indicated above, the FE simulations are carried out using Dirichlet boundary conditions (continuity in the displacement field). This leads to a discontinuous stress distribution across a dissimilar material interface, as illustrated in Fig. 1(b). Due to the axial nature of the deformation (along the z direction) and the spherical nature of the internal cavity, the stresses that develop on elements of PZT near the equatorial plane of the pore are of tangential character, while those on elements aligned with z (near the spherical caps of the pore) are of tensile character. The relative magnitude of these two types of stresses depends on the value of E_z and on the geometry changes in the material due to deformation. It is, thus, important to understand the effect that these stresses have on the device response. To that end, we perform a series of simulated tests aimed at isolating the stress component in each case and analyze its impact on the material.

First we setup a simple test where a layer of PZT of varying thickness is attached to a CFO solid block 10-nm in size in the manner shown in Fig. 7(a). Then, a voltage differential of 0.05 V is applied to the PZT layer so as to induce piezoelectric expansion and the deformation and stress of the system are obtained. This resembles the conditions in the region of the PZT layer under tangential stress. However, because CFO is not directly mechanically deformed, the elements adjacent or near the PZT/CFO interface are under shear stresses, as shown in Fig. 7(b). Interestingly, while the region of CFO near the interface deforms in the same manner as the PZT, at a distance away from it the trend inverses and the CFO contracts along the direction of the applied field. Thus, shear stresses in CFO created by tangential stresses in PZT can partially induce an inverse displacement-voltage correlation, which is more pronounced the farther away from the interface.

The results as a function of PZT-layer thickness (at a fixed value of $Ez=5.0$ MV m^–1) are shown in Fig. 8, where we plot the strain S_zz on the edge of CFO block farthest from the interface. As the graph illustrates, this strain is negative under shear conditions at the interface, with a maximum value of $−0.08 %$ for a PZT thickness of 6 nm. By contrast, application of a constant electric field to the geometry as shown in Fig. 9(a) results in positive monotonically increasing strains as the PZT layer thickness increases, as seen in Fig. 8. These computational tests, thus, reveal a mixed picture whereby (i) a diminishing strain response of the system with PZT thickness to an external electric field is observed when the CFO/PZT interface is oriented tangentially to the loading direction, and (ii) an increasing strain response when oriented along the field direction. Note that these tests in Figs. 7 and 9 only capture the effect of the interface coupling, but not the effect of porosity. Nevertheless, we have carried out simulations of a simplified geometry consisting of a simple cubic cell with a single central internal pore (see the supplementary material) and have found the effect of porosity to be marginal.

In summary, we have developed a computational finite-element model of a nanoporous multiferroic composite consisting a CFO template coated in the interior of the pores with PZT of various thicknesses. The scale of the model mimics that of the experimentally synthesized structures, with the model subjected to nominally identical boundary conditions as in experimental tests. We have considered both linear and nonlinear piezoelectric responses to infer the degree of poling of PZT prior to the application of an electric field. The nonlinear response can be attributed to extrinsic poling of PZT, resulting a much improved agreement with the measurements. We find that the observed non-monotonic response stems from two competing processes. First, increased porosity works toward increasing the strain due to a reduced system stiffness. Second, decreased porosity involves a larger mass fraction of PZT, which drives the electro-mechanical response of the structure, thus leading to a larger strain. The balance between these two driving forces is controlled by the shear coupling at the CFO/PZT interface and the effective PZT cross section along the direction of the applied electric field.

See the supplementary material file contains all the numerical information needed to parameterize Eqs. (2)–(4), given in the form of coefficient matrices of measured values. As well, we have included extra analysis to support our conclusion that porosity only has a secondary effect in the strain coupling of PZT and CFO.

The authors acknowledge the support from UCLA's Office of the Vice Chancellor for Research under seed grant for developing interdisciplinary research groups. Part of this work was supported by the National Science Foundation Nanosystems Engineering Research Center for Translational Applications of Nanoscale Multiferroic Systems (TANMS) under Cooperative Agreement Award No. EEC-1160504. Additionally, author S.K.P is supported by the National Science Foundation Graduate Research Fellowship under Grant Nos. DGE-1650604 and DGE-2034835. This manuscript contains data collected at the Stanford Synchrotron Radiation Lightsource (SSRL), experimental station 7-2. SSRL and the SLAC National Accelerator Laboratory are supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences under Contract No. DE-AC02-76SF00515.