We investigate numerically the homogenized permittivities of composites made of low-index dielectric inclusions in a ferroelectric matrix under a static electric field. A refined model is used to take into account the coupling between the electrostatic problem and the electric field dependent permittivity of the ferroelectric material, leading to a local field enhancement and permittivity change in the ferroelectric. Periodic and pseudorandom structures in two dimensions are investigated, and we compute the effective permittivity, losses, electrically induced anisotropy, and tunability of these metamaterials. We show that the tunability of such composites might be substantially enhanced in the periodic case, whereas introducing disorder in the microstructure weakens the effect of the enhanced local permittivity change. Our results may be useful to guide the synthesis of novel composite ceramics with improved characteristics for controllable microwave devices.

## I. INTRODUCTION

Ferroelectric materials play a crucial role in reconfigurable microwave devices with typical applications including antenna beam steering, phase shifters, tunable power splitters, filters, voltage controlled oscillators, and matching networks.^{1} Both bulk ceramics and thin films have been employed to design frequency agile components^{2–4} and metamaterials.^{5,6} The main reason for using ferroelectric materials is their strong dependence of their permittivity $\epsilon $ on an applied electric field $E$, which is measured by their tunability defined as $n=\epsilon (0)/\epsilon (E)$, along with a nonhysteresis behavior when used in their paraelectric state. The key requirements for antenna and microwave applications are large tunability and low losses. These two characteristics are correlated and one has to find a trade-off for optimal device performance, which can be quantified by the so-called commutation quality factor $K=(n\u22121)2/(ntan\u2061\delta (0)tan\u2061\delta (E))$, where $tan\u2061\delta $ is the loss tangent. These materials have usually high permittivity values even at microwave frequencies, often leading to slow response time and impedance mismatch, which can be an issue in some practical applications. Thus, it has been considered to mix ferroelectric ceramics with low-index and low-loss nontunable dielectrics in order to reduce both permittivity values and losses or to use porous ceramics to achieve the same goals without unwanted chemical reactions at the boundaries between dissimilar materials. In particular, the addition of magnesium oxide in barium strontium titanate (BST) ceramics has been shown to decrease the losses while keeping good tunability.^{7,8} Ceramics such as $Pb(Zr,Ti)O3$ (PZT) and $BaTiO3$ (BT) have been used as fillers in polymer-based composites with high dielectric constant.^{9} Other mixtures include metal–polymer composites^{10} and electroactive polymers such as poly(vinylidene fluoride) (PVDF) with high index dielectric inclusions.^{11}

The effective parameters of those composites have been investigated^{12–15} and it has been found that the permittivity can be greatly reduced while losses are much less sensitive to the dielectric phase addition and, in some situations, can lead to a small increase of the tunability of the mixtures. Analytical models based on the Bruggeman effective medium approach for low concentration of dielectrics were derived for different configurations (columnar, layered, and spherical inclusions models) and have been successfully compared with numerical simulations and experiments.^{12} In the context of porous ferroelectrics, the homogenized properties strongly depend on the size and morphology of the pores.^{16,17} Recently, the concept of tailoring the nonlinear properties of ferroelectric and dielectric structures by local field engineering has been introduced.^{18–20} It was shown through finite element calculation including the nonlinear coupling that, by employing composite materials made of linear dielectric inclusions into a ferroelectric matrix, one can lower the permittivity while maintaining high tunability, due to the local field in the ferroelectric phase that is tuned by the linear dielectric phase. Moreover, the effect of grain sizes in ferroelectric ceramics was studied using a model taking the field enhancement into account at the grain boundaries, and the predicted behavior was successfully compared to experimental data.^{19} Generally, there is a need for refined theoretical and numerical models to explain and design tunable materials and composites with tailored nonlinear properties. Note that the general method followed by our coupled model could be applied to other types of tunable system where local field enhancement and amplification is relevant, including ferromagnetic metamaterials,^{21} liquid crystals based devices,^{22} or field-enhanced carrier dynamics in doped semiconductors at other frequency ranges, particularly in the terahertz and near-infrared.^{23,24}

This study investigates numerically the effective permittivity of composites made of dielectric inclusions in a ferroelectric matrix by using a two-scale convergence method.^{25,26} The originality lies in the fact that a fully coupling model is employed to calculate the electrostatic field distribution when a uniform biasing field is applied on the structures, which will result in a local modification of the permittivity in the ferroelectric phase due to the microstructure. As compared to a simple uncoupled model where the ferroelectric phase is only modified through the biasing field, the resulting effective permittivity, dielectric losses, tunability, and anisotropy significantly differ. In contrast with earlier studies in the literature,^{18,19} we account for the nonlinear coupling beyond the first iteration and use a two-scale convergence homogenization analysis to obtain the effective parameters at higher frequencies, instead of a capacitance-based model valid in the static regime. This is an important point as, contrary to most homogenization procedures that are based on a quasistatic approximation, the two-scale convergence method fixes the frequency and lets the characteristic size of the system (the periodicity of the composites) tend to zero.^{26} This asymptotic analysis allows one to study the frequency dependence of the effective parameters. In addition, analytical models for the effective permittivity routinely employed in the literature, such as Maxwell-Garnett or Bruggeman theories, are limited to a few canonical shapes of the inclusions and cannot handle arbitrary geometries and media with spatially varying properties. This last point is of particular importance in the context of this study since we have to account for the field induced local permittivity change.

The model we developed has been implemented with the finite element method (FEM), and we realize a systematic computational study of ferroelecric-dielectric mixtures. First, we consider metamaterials consisting of a square array of parallel dielectric rods with circular cross section in a ferroelectric host and then investigate the effect of random distribution of those rods within the unit cell.

## II. THEORY AND NUMERICAL MODEL

We consider a composite made of a ferroelectric material with anisotropic permittivity $\epsilon f(E)$ that is dependent on an applied electric field $E$ and a nontunable dielectric of permittivity $\epsilon d$, which are both nonmagnetic. The structures under study are invariant along the $z$ direction, which leads to the standard decomposition of the wave equation in the transverse electric case (TE, electric field parallel to the direction of invariance) and the transverse magnetic case (TM, magnetic field parallel to the direction of invariance). A uniform biasing field is applied in order to be able to tune the effective permittivity. Modeling homogenized properties of this type of mixtures can be done by assuming that the electric field distribution is uniform throughout the sample so that the study of the tunability is essentially achieved by changing the value of the properties in the ferroelectric phase and computing the effective permittivity of the composite. We refer to this approach as the uncoupled model in the following. However, a more accurate description is to take into account the change of the electric field by the microstructure, if any. We, therefore, need to solve an electrostatic equation to find the field distribution within the material, but its solution depends on the permittivities of both materials, and the permittivity in the ferroelectric phase depends on this induced electric field: this leads to a strongly coupled problem.

### A. Permittivity model

We use barium strontium titanate (BST) as our ferroelectric material.$BaxSr1\u2212xTiO3$ samples were fabricated using the conventional sintering method with a barium ratio of $x=0.6$ to obtain a dielectrically tunable material as reported in the literature.^{11,27} The tunability was measured using an impedance analyzer up to 100 MHz and at 3.8 GHz using a loaded microstrip split ring resonator.^{28} The measured tunability of the in-house BST samples of 27% under 1 kV/mm DC bias was in agreement with those reported elsewhere.^{11,27} The method presented is, however, general and relies only on the gradient of the dielectric tunability vs electric field and could be applied to any tunable host material. The normalized permittivity value as a function of biasing field is reported in Fig. 1.

Case . | $\epsilon f(0)$ . | $\alpha $ ($\mu m2/V2$) . | $\beta $ ($\mu m4/V4$) . |
---|---|---|---|

Static | 3050 | 0.120 | 0.024 |

$f=3.8GHz$ | 165 | 0.240 | 0.079 |

Case . | $\epsilon f(0)$ . | $\alpha $ ($\mu m2/V2$) . | $\beta $ ($\mu m4/V4$) . |
---|---|---|---|

Static | 3050 | 0.120 | 0.024 |

$f=3.8GHz$ | 165 | 0.240 | 0.079 |

To describe the permittivity, we make use of the Landau potential given by $F(P,E)=F0+aP2/2+bP4/4+cP6/6\u2212EP$, where $E$ is the applied electric field and $P$ is the polarization.^{29,30} Variations of the permittivity with the temperature can be taken into account through the coefficients $a$, $b$, and $c$, but we assume that we are working at a constant room temperature. We further assume that the material is not subject to any stress so that the variation of permittivity due to mechanical constraints is irrelevant. The equation of state

gives the dependence of the polarization on the applied electric field, with $P0$ being the equilibrium polarization. Along the direction of a uniform applied electric field, the relative permittivity is given by

where $\epsilon f(0)=1/a$, $\alpha =3b/a$, and $\beta =5c/a$. The fitting parameters are given in Table I. As the norm of the field increases, the permittivity decreases with a characteristic bell curve typical for a ferroelectric material in its paraelectric state. Furthermore, assuming that the crystalline principal axes of the ferroelectric material are oriented in the coordinate directions and that the diagonal components of the permittivity tensor are only function of the corresponding bias electric field components,^{31} we have

where each of the diagonal components has the functional form given by Eq. (1). Note that we will use the static values of permittivity for the electrostatic modeling, while we are interested in the homogenized values of permittivity at microwaves.

### B. Electrostatic model

The composites under study are made of two materials, thus their permittivity is represented by a piecewise defined tensor $\epsilon (r,E)$, which is equal to $\epsilon f(E(r))$ in the ferroelectric phase and $diag(\epsilon d)$ in the dielectric phase. In the following, we consider two different cases for the biasing field. Because of the form (2) assumed for the ferroelectric permittivity tensor, $\epsilon zz$ will not be changing for a field in the plane orthogonal to the $z$ axis. This is the only component being relevant for TE polarization, so we consider in this case a uniform biasing electric field applied along the direction of invariance $E0=E0ez$. On the other hand, the in-plane components of $\epsilon f$ are tuned by $Ex$ and $Ey$, therefore, without loss of generality, we consider a uniform applied electric field directed along the $x$ axis $E0=E0ex$ for the TM polarization case. To calculate the total electric field in the material, one has to solve for the potential $V$ satisfying Gauss’ law

Note that for the TE case, the solution is trivial since the structures are invariant along $z$ so that the electric field is equal to the uniform biasing field, and we will thus not study it in the following. However, in the TM case, the situation is much more complex: this is a coupled problem since the electric field $E=\u2212\u2207V$ derived from the solution of Eq. (3) depends on the permittivity distribution, which itself depends on the electric field. The coupled system formed of Eqs. (2) and (3) is solved iteratively until there is convergence on the norm of the electric field. Here, we would like to emphasize that the permittivity in the ferroelectric material, although uniform initially, is spatially varying due to the nonuniform distribution of the total electric field.

### C. Homogenization

When the period of the composite metamaterial is much smaller than the wavelength, one can describe the properties of the composite by a bulk medium with homogenized parameters. The effective permittivity for TM polarization is calculated using a two-scale convergence homogenization technique.^{25,26} For this purpose, one has to find the solutions $\psi j$ of two annex problems $Pj$, $j={1,2}$,

where $r=(x,y)T$ is the position vector in the $xy$ plane and $\xi =\epsilon T/det(\epsilon )$. The homogenized tensor $\xi ~$ is obtained with

where $\u27e8\u22c5\u27e9$ denotes the mean value over the unit cell. The elements of the matrix $\varphi $ represent correction terms and are given by $\varphi ij=\u27e8\xi \u2207\psi i\u27e9j$. Finally, the effective permittivity tensor can be calculated using $\epsilon ~=\xi ~T/det(\xi ~)$.

Note that the TE case, which we shall not study here as no coupling happens, is trivial since the homogenized permittivity is simply the average of the permittivity in the unit cell: $\epsilon ~=\u27e8\epsilon \u27e9$.

## III. NUMERICAL RESULTS

In the following numerical results, the dielectric phase is supposed to be lossless and nondispersive with $\epsilon d=3$, while the ferroelectric material follows the permittivity described in Sec. II A and has a constant loss tangent, $tan\u2061\delta f=10\u22122$. Equations (3) and (4) are solved with a finite element method using the open source packages Gmsh^{32} and GetDP.^{33} In both cases, we use a square unit cell $\Omega $ of length $d$ with periodic boundary conditions along $x$ and $y$. Second order Lagrange elements are used, and the solution is computed with a direct solver (MUMPS^{34}).

### A. Two dimensional periodic metamaterial

Let us now consider a periodic square array of infinitely long dielectric rods of a circular cross section of radius $r$ embedded in a ferroelectric matrix.

We first study the convergence of the coupled problem on the particular case with dielectric filling fraction $f=\pi r2/d2=0.5$ and $E0=2MV/m$. Figures 2(a) and 2(b) show the convergence of the real part and loss tangent of the components of the homogenized permittivity tensor, respectively. The $yy$ components converge quickly and are almost unaffected by the coupling process, whereas the $xx$ components change substantially from the initial conditions. This is due to the effect of the redistribution of the electrostatic field within the unit cell [see Figs. 2(c) and 2(d)], where the $x$ component of the electric field is still much stronger than the $y$ component, even if spatially varying in the ferroelectric medium. At equilibrium, the electric field is concentrated close to the $y$ axis in between two neighboring rods. This, in turn, affects the permittivity distribution [see Figs. 2(e) and 2(f)] and the homogenized properties of the composite.

We computed the effective parameters of these metamaterial structures for different radii of the rods and studied their behavior when subjected to an external electrostatic field (see Fig. 3). The results of our coupled model differ significantly from the uncoupled one. Increasing the dielectric fraction lowers the effective permittivity, while the losses are slightly reduced but much less sensitive. Due to the inhomogeneous redistribution of the permittivity over the ferroelectric domain, the overall tunability changes. For the periodic composites studied in this section, taking into account the coupling leads to an effective tunability increase with higher dielectric concentration, and that is larger than the tunability of bulk ferroelectric. This can be seen in Fig. 3(c), where we plot the tunability of the composites along the $x$ axis, $n~(E)=\epsilon ~xx(E)/\epsilon ~xx(0)$, normalized to the tunability of the bulk ferroelectric $n(E)=\epsilon xxf(E)/\epsilon xxf(0)$. Two concurrent effects are at stake here: on the one hand, the dilution of ferroelectric makes the composite less tunable, but on the other hand, the rearrangement of the electrostatic field surrounding the inclusion and its concentration in some region will cause a higher permittivity change locally. The relative strength of those phenomena is governed by the shape of the inclusion and its permittivity, and so, it is envisioned that the performance of the composites might be enhanced by engineering their microstructure. Those observations are consistent with previously published numerical and experimental results^{18} where the local field enhancement in porous ferroelectrics has been shown to possibly increase tunability with reducing permittivity for small porosity levels. Our approach also agrees with an analytical spherical inclusion model predicting an increase of the tunability with the dilution of the ferroelectric.^{12}

The geometry of the unit cell is symmetric so the homogenized material is isotropic when no field is applied. But when the sample is biased, the permittivity distribution becomes asymmetric due to the inhomogeneity of the electric field, thus making the effective material properties anisotropic. This geometric effect is added to the anisotropy arising from the material properties of the ferroelectric phase itself and depending on the topology and permittivity of the rods, one effect would be predominant. In the case studied here, the equilibrium permittivity distribution varies strongly along the bias direction and much less orthogonally to it, which adds anisotropy by diminishing the effective permittivity in the $x$ direction. This local field induced effect is what makes the anisotropy stronger in our coupled model compared to the uncoupled one [cf. Fig. 3(d), where we plot the anisotropy factor $\alpha =\epsilon xx/\epsilon yy$]. Those subtle phenomena can only be rigorously taken into account by employing a coupling formalism and are responsible for the difference observed when compared to a simple uncoupled model.

### B. Pseudorandom case

We finally study the effect of random distribution of the inclusions within the unit cell on the effective parameters of the composites. This is an important point as fabrication of randomly dispersed inclusions is much more easy from a technological perspective. For each filling fraction of the dielectric, we generated 21 numerical samples with inclusions of circular cross section of average radius $r=d/20$ that can vary by $\xb130%$. Their center is chosen randomly and the rods are allowed to overlap. An example of distribution for $f=0.5$ is given in Fig. 4. The effective material properties are plotted in Fig. 5. Similar to the periodic case, the permittivity decreases with increasing dilution of ferroelectric, but for identical filling fraction, the permittivity is lower as compared to the periodic array, and the smaller the dielectric concentration, the larger is the difference. Losses decrease as well and the reduction is substantially larger than the periodic case, with higher variation from sample to sample as $f$ increases. The effective tunability is on average smaller than that in the periodic case, and for low biasing fields and for some particular samples, it can be greater than the bulk tunability. However, at higher applied electric fields, normalized tunability becomes smaller than unity and is reduced as one adds more dielectric. For comparison, the homogenized parameters are plotted in Fig. 6 in the case where the coupling is neglected. One can see that the coupled and uncoupled models give similar results for the tunability, whereas the losses are still smaller for the coupled case at higher fields.

The redistribution of electric field, permittivity, and convergence of the effective parameters are displayed in Fig. 7. The effect of disorder plays an important role here: the electrostatic field gets concentrated in between neighboring inclusions and the smaller the gap, the higher the field, hence a greater local permittivity change. In addition, even if the distribution is random, one expects that the anisotropy due to geometry would cancel for a sufficiently large number of rods (which is the case as the mean anisotropy factor is close to $1$ when no bias field is applied). However, the anisotropy due to ferroelectric properties is important in this case as well, as both the $x$ and $y$ components of the electrostatic field are playing a role. Because of the relative positions of the rods, both $\epsilon xx$ and $\epsilon yy$ are affected by the coupling so that the anisotropy factor for higher fields is reduced as compared to the periodic case. However, even if there is a substantial variability from sample to sample, on average, the anisotropy factor decreases with increasing dielectric concentration.

## IV. CONCLUSION

We have studied the homogenized properties of dielectric/ferroelectric mixtures using a rigorous model that takes into account the coupling between the electrostatic field distribution and the field dependent ferroelectric permittivity tensor. After convergence of the coupled problem, the effective permittivity tensor is calculated using a two-scale convergence homogenization theory. The results obtained by this model differ significantly from a simple assumption that the permittivity of the ferroelectric responds just to the uniform biasing field. We have considered both periodic and random arrays of dielectric rods in a ferroelectric matrix in 2D and studied their effective properties for TM polarization as a function of dielectric concentration and bias field. Importantly, adding more low index and low loss dielectric allows to decrease the overall permittivity significantly and slightly lower the losses. For the periodic case, the tunability is higher than the bulk due to local field enhancement, whereas this effect is strongly suppressed when the disorder is introduced. The asymmetric redistribution of the permittivity induces an effective anisotropy that is added to the one arising purely from the ferroelectric material. The properties of the composites are affected by multiple factors: geometry and the spatially dependent electric field that will induce locally a tunable, anisotropic response in the ferroelectric phase depending on its amplitude and direction. This suggests that the performances of the composites may be enhanced by distributing the two phases in an optimal way to get high tunability and low losses. Further work in that direction is needed as well as extending this study to 3D media. Finally, because the permittivity of the dielectric is much smaller than the ferroelectric one, it would be of great interest to use high contrast homogenization theory^{35,36} to study these kinds of mixtures. This would reveal the frequency dependent artificial magnetism due to “microresonances” in the high index phase and potentially lead to composites with tunable effective permeability.

## ACKNOWLEDGMENTS

This work was funded by the Engineering and Physical Sciences Research Council (EPSRC), UK, under a grant (No. EP/P005578/1) “Adaptive Tools for Electromagnetics and Materials Modelling to Bridge the Gap between Design and Manufacturing (AOTOMAT).”

The authors would like to thank Henry Giddens for performing the measurements of ferroelectric permittivity used in this paper.

The codes necessary to reproduce the results in this article are freely available online at this address: https://www.github.com/benvial/ferromtm.

## REFERENCES

*14th IEEE International Symposium on Applications of Ferroelectrics, ISAF-04, 2004*(