High frequency meta-ferroelectrics by inverse design

: Composites with subwavelength features exhibit effective properties that depend on microstructure morphology and materials, which can be adjusted to obtain enhanced characteris-tics. We detail the systematic design of electromagnetic metamaterials composed of dielectric inclusions in a ferroelectric matrix that, under an applied voltage, present an optimized effective tunability higher than the bulk due to a nonlinear local electric field enhancement. The effect of volume fraction, losses, and biasing field on homogenized properties is investigated and the analysis of the photonic band diagram is carried out, providing the frequency dependence of the anisotropic effective index and tunability. Such metaceramics can be used in microwave antennas and components with higher reconfigurability and reduced power consumption.


Introduction
Multifunctional and reconfigurable systems are increasingly demanded for a wide range of applications in Electromagnetics and Photonics [1][2][3]. In parallel, the development of artificial media and metamaterials [4,5], allowing unprecedented control of electromagnetic (EM) waves by carefully engineered subwavelength structures, has stimulated the research interest in tunable materials. Indeed, the ability to introduce tunable components has been established as a straightforward route to EM and metamaterials reconfigurability, albeit the added fabrication complexity in comparison with static devices. As a matter of fact, it is envisioned that reconfigurable metamaterial structures have the potential to mitigate some of the most significant limitations of static metamaterials such as narrow bandwidth operations, high losses, and tolerance sensitivity. Often, in addition to the necessary control logic, adding reconfigurability to a metamaterial structure requires the introduction of extra components and materials to a standard design, hence potentially affecting the system's performance, power consumption and weight. Consequently, there are still significant theoretical and technical obstacles to the production and deployment of reconfigurable metamaterials within commercial products, particularly regarding versatility, operational bandwidth, performance, complexity of realization, robustness, speed of tuning, and costs. Those engineering challenges are driving efforts of academia and industry towards new technological solutions and experimental testing of reconfigurable metamaterials, which is currently a very active area of research. Experimental demonstrations across frequency ranges have been reported and include microwave devices [6,7], lenses [8,9] tunable band diagrams [10,11], reconfigurable orbital angular momentum [12], topological insulators [13], infrared emitters [14] negative effective index media [15][16][17], hyperbolic metamaterials [18] and optically controlled metasurfaces [19] using different stimuli and tuning mechanisms such as mechanically adjustable elements [20][21][22][23][24], microfluidics [25,26], graphene [27] or liquid crystals [28,29].
In microwave engineering, ferroelectric materials play a crucial role in applications needing reconfigurability, such as antenna beam steering, phase shifters, filters, and tunable power splitters, [30]. Both thin films and bulk ceramics are used for frequency agile components [31,32] and metamaterials [15,33]. Ferroelectric materials possess a strong nonlinear dependence of their permittivity ε on an applied electric field E, with key requirements for antenna and microwave applications being large tunability and low losses. The permittivity values are usually high even at microwave frequencies, and therefore it has been considered to mix ferroelectric materials with low-index and low-loss non-tunable dielectrics or to synthesize porous ceramics in order to reduce both permittivity values and losses. The effective parameters of those composites have been investigated [34][35][36][37] using analytical effective medium approaches for low filling fraction of dielectric and have been successfully compared with numerical simulations and experiments. It has been shown that the permittivity can be reduced while losses are much less sensitive to the dielectric phase addition, and in some situations lead to a small increase of the tunability of the mixtures. Recently, refined theoretical and numerical models have been proposed [38][39][40], and the predicted behaviour compared well to experimental data [41]. Those methods account for the nonlinear coupling of electrostatic biasing and field dependent permittivity. It was shown that the permittivity can be reduced while maintaining high tunability and low losses, due to local field enhancement in the ferroelectric phase.
We propose here to design the microstructure of metamaterials composed of a ferroelectric and low index non tunable dielectric (see Fig. 1) through topology optimization [42]. Our aim is to find structures with desired effective properties: high tunability and minimal losses. Using a nonlinear coupled model that takes into consideration the subwavelength field enhancement and ferroelectric response, we systematically obtain the material distribution in the unit cell and study the effect of volume fraction and biasing field strength. Furthermore, we conduct a multi-objective optimization for maximal tunability and minimal losses and present the Pareto front, revealing a trade-off between the two characteristics. We finally study the photonic bands of a particular example of optimized geometry, showing an adjustable bandgap, derive the frequency dependence of the effective permittivity and compare the results with the homogenization approach. The structure is periodic along x and y, and is composed of a dielectric inclusion with permittivity ε d in a ferroelectric background with a field dependent permittivity ε f (E). A uniform electric field E B is applied to the metamaterial along the x direction.

Methods
The study is restricted to two-dimensional metamaterials with square unit cell Ω of size d<<λ and periodic boundary conditions, where λ is the wavelength under consideration. Since we target the microwave frequency range (f = 3.8 GHz, λ ≃ 8 cm), the maximal value for the periodicity is roughly λ/10 = 8 mm. On the other hand, we consider a bulk ferroelectric so that nano-scale effects such as domain walls motion [43], and micro-scale grains (ℓ ≃ 10 µm, [44]) are not resolved in the present model but are incorporated in the dielectric permittivity. This places a minimal length scale for the validity of the study at around 10ℓ ≃ 100 µm. The materials we employ are barium strontium titanate (BST) and a low index (ε d = 3) dielectric such as common thermoplastic polymer acrylonitrile butadiene styrene (ABS).

Coupling nonlinear ferroelectric permittivity and electrostatics
The tunability of barium strontium titanate Ba x Sr 1−x TiO 3 with a barium ratio of x = 0.6 samples fabricated using a conventional sintering method was measured in the static regime (DC) and at microwave frequencies (MW) [44] and fitted to a Vendik model [45] to represent the electric field dependence of the permittivity: with ξ = E/E N , ε f (0) = 3050 and E N = 1.65 kV/mm at DC, ε f (0) = 165 and E N = 1.12 kV/mm at f = 3.8 GHz. Effects of temperature and mechanical stress on the permittivity are neglected to simplify the model and focus on the electromagnetic properties of the composites. While multiple physical processes should be included since ferroelectric properties depends strongly on temperature, this is out of the scope of this paper.
To calculate the total electric field in the metamaterial, one has to solve for the potential V satisfying Gauss's law for a given permittivity distribution ε: biased by a constant uniform electric field E B = E B x. The electric field E = −∇V derived from the solution of Eq. (2) depends on the permittivity distribution, which itself depends on the electric field. The coupled system formed of Eqs. (1) and (2) is solved iteratively using a Picard method until convergence on the norm of the electric field ||E i+1 − E i ||<10 −5 . Although initially constant, the permittivity in the ferroelectric phase is spatially varying due to the non-uniform distribution of the total electric field.

Homogenized properties
The aim of homogenization approaches is to retrieve a simplified, but almost equally precise, description of the material response by averaging out material properties at the subwavelength level. Research in this area has been deeply linked to and influenced by the development of metamaterials and their application in physics, particularly for electromagnetic phenomena. Analytical models for the effective permittivity previously 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, which has to be accounted for in our model because of the field induced local permittivity change.
In this study, the effective permittivity of the metamaterials is calculated using a two scale convergence homogenization technique [46,47]. We will focus on 2D geometries and TM polarization (since the TE polarization effective permittivity is trivially the arithmetic mean of the basic components permittivity), resulting in a scalar wave equation for the z component of the magnetic field. The main idea is to select two scales in the study: a microscopic one (the size of the unit cell) and a mesoscopic one (the size of the bulk), controlled by a real parameter ν>0. Denoting r = (x, y) T the position vector, the approach consists in introducing the ansatz H (r) = H 0 (r, r/ν) + νH 1 (r, r/ν) + ν 2 H 2 (r, r/ν) + · · · for the magnetic field H solution of time harmonic Maxwell's equation and performing an asymptotic analysis as ν → 0. One then needs to find the solutions ψ j of two annex problems P j , j = x, y: The homogenized permittivityε is then obtained as: where ⟨.⟩ denotes the mean value over the unit cell, I is the 2 × 2 identity matrix and φ ij = ⟨ε∇ψ i ⟩ j are correction terms. This analysis allows us to obtain the effective parameters at higher frequencies, in contrast with capacitance-based models valid in the electrostatic regime. Contrarily to most homogenization procedures that are based on a quasi-static approximation, in the two scale convergence method, the frequency is fixed and the characteristic size of the system (the periodicity of the composites) tends to zero, making the study of the frequency dependence of the effective parameters possible.

Topology optimization
Inverse design, where a specific target is searched according to engineering constraints, is an actively researched topic in the nanophotonics and metamaterial community [48][49][50]. Topology optimization has led to the design of a large range of devices such as invisibility cloaks [51,52], illusion devices [53], metasurfaces [54,55] and metamaterials with tailored properties [56][57][58][59], showing non intuitive material distributions, allowing to manipulate light matter interaction with unprecedented capability. In the topology optimization procedure, the permittivity distribution is parametrized by a continuous scalar density function ρ ∈ [0, 1]. To avoid small features and pathological "chessboard" patterns, we use a filtered density obtained by solving the following Helmholtz equation [60] with periodic boundary conditions on Ω: with R a real positive parameter characterizing the filter radius (R = 0.05d in the rest of the paper). In order to enforce a binary design, we apply a threshold projection [61,62]: with ν = 1/2 and β a real positive parameter characterizing the strength of the projection and is increased during the course of the optimization process. Finally, the permittivity is given by the solid isotropic material interpolation (SIMP) method [63] as: where q is a penalty factor (q = 1 here), ε min = ε d = 3 (the permittivity of the dielectric phase) and ε max = ε f (E) (ferroelectric phase). The gradient based optimization is initialized with a density ρ 0 and performed for 20 iterations or until convergence on the objective or design variables. This step is then repeated for n global iterations setting β = 2 n , where n is an integer between 0 and 7, and initializing the algorithm with the optimized density obtained at the previous step.

Results
Numerical results are obtained using Python and open source libraries: Eqs. (1), (2), (3) and (5) are solved by the finite element method with the fenics library [64], the sensitivity of the objective g with respect to the design variables p are computed with the adjoint method [65] using automatic differentiation through dolfin-adjoint [66] and the method of moving asymptotes (MMA, [67]) is employed for the gradient based optimization with the nlopt package [68].

Maximizing tunability
The first objective is to maximize the tunability of the composites along the direction of the applied electric field defined asη(E) =ε xx (0)/ε xx (E): where σ is the tunability enhancement compared to the bulk BST tunability We enforce a constraint on the volume fraction f = ⟨ρ⟩ of ferroelectric (i.e. the proportion of BST in the metamaterial). Indeed, there is an infinity of possible material distribution with the same volume fraction, but they would lead to different effective material properties. Our aim is to find the distribution that maximize our optimization objective. Those constraints are imposed as an inequality because the MMA algorithm we use cannot enforce equality constraints. To target approximately a given volume fraction f cons , we set f min = f cons − δ f and f max = f cons + δ f and δ f is decreased linearly from 0.1 to 0.01 during the course of the optimization. The optimization is initialized with a density ρ 0 = 1 − exp(−(x 2 + y 2 )/r 2 f ), with r f =

√︁
(1 − f )/π log 2). Additionally, our numerical experiments have shown that field enhancement is suppressed when the inclusion are connected along the y-axis, we thus constrain further the material to be ferroelectric (i.e. ρ = 1) at the top and bottom of the unit cell in a strip of height 0.05d.
Results are displayed on Fig. (2(a)), where we plot the effective tunability enhancement of structures optimized for different volume constraints and a biasing field of 1 kV/mm. Tunability higher than bulk BST is achieved in each case except for the smallest volume fraction where it is 4% weaker. A maximum enhancement of 34% is achieved for a structure containing 83% of ferroelectric, which indicates a trade-off between composites containing more tunable material and the volume of dielectric phase which induces a field concentration within the unit cell. Note that the obtained effective properties are anisotropic due to to the asymmetry of material distribution, but we are only focusing on the properties parallel to the applied electric field. As f increases, both components of the effective permittivity tensor increase (see Table (1)), with strong anisotropy for lower values of filling fraction, which correspond to structures connected along the x-axis. As the target volume fraction increases, the topology evolves from layered structure along y to cross like shapes, a "bow tie" looking dielectric structure for the highest enhancement (f = 83%) to a needle like inclusion with long side along y for f = 83%. These non-trivial topologies show the generality and versatility of the approach which does not require any geometrical parametrization, generating freeform microstructures. The field maps of the electric field and corresponding permittivity in subfigures (2(b-e)) show that the field enhancement is maximal between the dielectric parts of the microstructure, which leads to a localized tuning in the ceramic phase. The previous optimization was performed for a fixed biasing field of 1 kV/mm. In order to assess the response of the metamaterials, we computed the effective tunability as a function of E (see Fig. (3)). The maximal values of tunability enhancement are obtained as expected around 1 kV/mm. We note that the permittivity model has an inflexion point around E = 1.18 kV/mm so that a larger local field variation leads to a stronger change in permittivity. The improved tunability is maintained, albeit smaller, for different biasing intensities except at higher fields. Further optimization studies not reported here with lower and higher field strength led to similar material topologies, indicating the robustness of the design process.

Bi-objective optimization: high tunability and low loss metamaterials
Our next aim is to find metamaterial structures that will at the same time enhance the effective tunability while reducing the effective losses. We now use a complex permittivity value ε C f with the same dependence on the electric field (1) but with constant value of the loss tangent tan δ f = 10 −3 , that is ε C f (E) = ε f (E) (1 − tan δ f ). We define the loss reduction factor as: where the homogenized loss tangent xx component is tanδ xx (E) = −Imε xx (E)/Reε xx (E). We then have to solve the following maximization problem: where the weighting parameter p ∈ [0, 1] is varied to give more importance to one objective or the other. Note that here the effective tunability is defined with the real parts of the homogenized permittivity:η(E) = Reε xx (0)/Reε xx (E). The applied electric field is 1KV/mm and the target volume fraction 50%. The Pareto front of this bi-objective problem is given in Fig. (4), where each point represent a different value for p. One can notice that there is a trade-off between those two concurrent goals and observe two regimes: the first with enhanced tunability (18% to 25%) and moderate loss reduction (18% to 24%), and the second where the tunability is weaker than bulk BST (-11% to -7%) but with a stronger decrease in losses (33% to 43%). Depending on the application requirements, one may choose between different topologies to mitigate the advantages and drawbacks of the two material properties. Note that the Pareto front is effectively discontinuous, and is jumping in between the two regimes of high tunability/moderate loss and low tunability/low losses. This is unrelated to the parameter p which is just a numerical coefficient used in the optimization procedure. We now look at the response of the optimized composites as the biasing field is varied. A metric for the performance of ceramics that takes into account tunability and losses is the commutation quality factor (CQF) defined as: For the structures associated with the second region of the Pareto set (p = 0.10 and p = 0.25), the tunability decreases significantly as the field increases (see Fig. (5(a))), while losses are strongly reduced, and this reduction actually increases with biasing strength (see Fig. (5(b))). We show in Fig. (5(c)) the normalized CQF defined as κ =K/K f , whereK and K f correspond to the CQF for the homogenized metamaterial and bulk BST respectively. For p = 0.10, the QCF is reduced while it is of the same order than than the bulk for the case p = 0.25, for all electric fields. This highlights the limits of such a metric where it is hard to distinguish independently the contributions of losses and tunability: indeed one can have, as the later case, a material with reduced tunability and losses with similar QCF. For the metaferroelectrics belonging to the other optimized region (p = 0.35 and p = 0.95), a maximum of tunability enhancement is observed as expected around 1KV/mm and decreasing for higher fields, similarly to the lossless case. Losses are reduced significantly and the stronger the bias, the smaller the loss tangent. Finally, in both cases the associated QCF is much larger than the bulk (around 4 times at 1KV/mm) and decreases as E increases. The enhanced performances reported here are clearly a trade-off in between those two conflicting objectives (losses and tunability), but other considerations may come into play in the choice of and adequate metamaterial such as the volume fraction or the effective anisotropy.

Photonic bands properties
We now investigate the spectral properties of one example of optimized composite (Section (3.1) for f = 0.47). The eigenvalues k n = ω n /c, of Maxwell's operator for TM polarization, i.e the solutions of: with c the speed of light in a vacuum, ω n the eigenpulsations and H n the eigenmodes, subject to Bloch-periodic conditions with wavevector k = (k x , k y ) T . This spectral problem is solved using the open-source finite element package GetDP [69] with its interface to the SLEPc library [70]. The effective permittivity can be obtained from the band diagram according to: Figure (6(a)) shows the photonic band diagram associated with the optimized lossless composite for f = 0.47 obtained in Section (3.1) (see inset in Fig. (2) for the unit cell topology) and observe a tunable band gap between the first and second bands. When unbiased, the photonic band gap is centred at ω = 0.0829ω 0 with a width of w = 0.0173ω 0 , and when the electric field is applied, its centre is blueshifted to 0.0848ω 0 (∼2%) and slightly broadened (w = 0.0175ω 0 , ∼0.6%). Another much narrower bandgap between the second and third bands is present as well, with centre ω = 0.114ω 0 (resp. 0.118ω 0 ) and width 0.0017ω 0 (resp. 0.0019ω 0 ) for the unbiased and biased case respectively. Note that while in the first case the two band gap overlap when the electric field is on or off, in the later case the two forbidden propagation bands are disjoint since the bandwidth is smaller. The isofrequency contours for the first band are displayed on (6(b)). When no field is applied (full lines), we have elliptical contours with long axis along k y , meaning a diagonal anisotropic effective index with higher values along x, which is consistent with the results of the previous homogenization study. This is valid for ω → 0, and when the frequency increases so does the effective index. Applying a voltage on the photonic crystal deforms the isofrequency contours, increasing the semi-axis of the ellipse along k y but almost keeping the same along k x : in other words, the tunability is high in the x direction but practically negligible in the orthogonal direction, in agreement with the results from Section (3.1).
To corroborate those qualitative results, we extracted the effective permittivity tensor from the first band and plotted the results in Fig. (7(a)). Firstly, the values agree very well with those obtained with two scale homogenization (see the markers at ω = 0), and we observe an increase of the components of the effective permittivity. The permittivity along x is enhanced further as we get closer to the edge of the gap, and hence, due to the shift of the bandgap from applying a voltage, the tunability along x is increased, reaching a value of around 4 at ω = 0.041ω 0 (see Fig. (7(b))). This effect, even if it is narrowband, provides a way to even further enhance the tunability in this kind of metastructures. Note that the design was optimized to enhance the homogenized tunability, not to maximize bandgaps / tunability frequency dependence, which is out of the scope of this paper.

Conclusion
We presented an inverse design methodology to obtain microstructures of ferroelectric/dielectric composites with desired and enhanced properties. This inverse homogenization approach allows us to find the material distribution within the periodic unit cell that gives a maximized effective tunability for a given volume fraction of ceramic. Tunability enhancement as high as 34% compared to bulk BST for a filling fraction of 87% are obtained. Taking dissipation into account, we have shown that one has to mitigate two competing objectives, namely high tunability and low losses. In addition, we have studied the spectral properties of one optimized photonic crystal, showing a tunable bandgaps, frequency dependent enhanced tunability and confirming the homogenization approach. Extension to 3D metamaterial topologies, fabrication and experimental validation of the designed composites will be the subject of future studies. The method proposed here provide guidelines for the systematic design of metamaterials with desired properties, allowing further control of electromagnetic waves in antennas and radio frequency applications. The proposed designs can provide a range of permittivity values suitable for instance for the implementation of graded index lenses employed in beam steering antennas. At the same time, the optimized tunability would increase the reconfigurability of such devices, avoiding the need for complex feed network. On the other hand, a co-design of and integrated biasing network of electrodes will be needed. The general approach followed here can be applied to other type of tunable systems in different frequency ranges where local field enhancement and amplification is relevant, including for instance liquid crystals based devices, ferromagnetic metamaterials and field-enhanced carrier dynamics in doped semiconductors. Finally, since ferroelectric material behaviour is inherently governed by multiple physical processes, it would be of great interest to optimize in parallel the effective electro-mechanical and piezoelectric properties of metaferroelectrics, in order for example to mitigate the brittleness of ceramic materials and obtain composites with high tunability, low losses and strong mechanical compliance.