{'fraction_non_alphanumeric': 0.032017075773746, 'fraction_numerical': 0.011970348149643774, 'mean_word_length': 4.709815546772068, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Maximizing the electromomentum coupling in piezoelectric laminates

Asymmetric piezoelectric composites exhibit coupling between their macroscopic linear momentum and electric field, a coupling that does not appear at the microscopic scale. This electromomentum coupling constitutes an additional knob to tailor the dynamic response of the medium, analogously to the Willis coupling in elastic composites. Here, we employ topology- and free material optimization approaches to maximize the electromomentum coupling of periodic piezoelectric laminates in the low frequency, long-wavelength limit. We find that the coupling can be enhanced by orders of magnitude, depending on the degrees of freedom in the optimization process. The optimal compositions that we find provide guidelines for the design of metamaterials with maximum electromomentum coupling, paving the way for their integration in wave control applications.


Introduction
Metamaterials are artificial composites whose effective behavior is fundamentally different from the behavior of their constituents [1][2][3][4][5][6][7][8][9][10]. Of particular relevance to this work are Willis metamaterials , whose linear momentum and strain is coupled, and so are their stress and velocity: these couplings do not appear in homogeneous media. They are captured by the so-called Willis tensor that enters the effective constitutive relations of the composite, reflecting a designable degree of freedom to tailor its dynamic response. The Willis tensor is nonlocal in space and time; its optimization problems, we restrict attention to periodic cells that are made of three layers. Since this optimization is based on a small number of discrete variables, we use a genetic algorithm to solve the optimization problem. We find that the electromomentum coefficient of the optimal laminate is two orders of magnitude larger than the laminate that was first considered by Pernas-Salomón and Shmuel [34], which was chosen arbitrarily. Henceforth, we refer to the latter laminate as the reference laminate. We observe that the optimization yields a choice of materials that tends to maximize the contrast between their electromechanical properties, in a unit cell that comprises one thin layer in-between two thick layers. In Section 4 we consider a design space with continuous variables, whose values are bounded in-between the extreme values of the materials considered earlier. We solve different optimization problems using gradient-based algorithms that are implemented together with a sensitivity analysis of the objective function. We first consider optimization problems whose design variables are pairs that consist the mass density and one of the electromechanical properties, i.e., the dielectric-, piezoelectric-and stiffness coefficients, while all the rest of the properties are set to the properties of the reference laminate. We find that each one of the optimal pairs enlarge the electromomentum coefficient by another order of magnitude, in comparison with the optimal laminate that is obtained from the DMO. We then proceed and increase the number of design variables, and find that when all material properties are designable, the electromomentum coefficient is enlarged by another order of magnitude. We conclude this paper in Section 5, with a summary of our main results.

The homogenization process
The homogenization process that we use to extract the electromomentum coefficient is founded on two principles. (i) The macroscopic fields are defined by the product of the Bloch envelope of the microscopic fields and the volume average of their periodic part. Defined in this way, the effective fields satisfy identically macroscopic governing equations that are of the same form as the microscopic equations, see Refs. [34,46,47]. (ii) We account for several driving sources in the formulation [34,46,[48][49][50]; collectively, these sources allow us to obtain a unique set of effective properties [17,[51][52][53], which satisfies necessary physical laws [16,37,52,54,55].
Our averaging process uses the plane wave expansion method [18,23,[56][57][58], and the resultant scheme is essentially a one-dimensional reduction of the scheme that was developed by Muhafra et al. [37]. In order to make this manuscript self-contained and introduce the quantities that are used in the optimization process, the complete one-dimensional scheme is provided next.
Consider a periodic repetition of three piezoelectric layers, the poling direction of which is in the lamination direction, say x. The laminate is driven by a body force density f , axial inelastic strain η, and free charge density q, of the form where {s 0 } are constant. These sources generate a longitudinal motion, u (x,t), governed by the equations here, σ is the Cauchy stress, p is the linear momentum, and D is the electric displacement. The remaining Faraday equation for electric field E is satisfied by setting E = −φ ,x , where φ is termed the electric potential. The microscopic kinetic and kinematic fields in the laminate are related via the constitutive relations where C, B, A and ρ are the stiffness, piezoelectric coefficient, dielectric coefficient and the mass density, respectively. These properties are periodic, such that where n ∈ Z; l T = ∑ 3 m=1 l m and l m is the length of layer m, as shown in Figure 1. This periodicity implies that the kinetic and the kinematic fields are of the Bloch form The objective is to determine the effective constitutive relations, relating the macroscopic fields, which have the form To this end, we first consider the Bloch form of Eq. (3), namely, To relate the mean parts of the kinetic and kinematic fields, we expand the periodic functions in Eq. (7) in Fourier series, such that where {ψ G } are the Fourier coefficients ofψ (x); G = 2π l T m and m ∈ Z. Substituting Eq. (8) into Eq. (7) yields where G is defined in the same manner as G and Eq. (9) contains every combination of G and G. Since Eq. (9) holds for any value of x and t, the sum in the square brackets must be zero. We multiply those sums by e −iG x , and integrate the result over the unit cell. Due to the orthogonality property of the Fourier series, the only non-vanishing terms are those satisfying the condition G = G + G . Hence, substituting G = G into the resultant equation yields where (•) GG denotes the Fourier coefficient of (•) along the basis function e i(G−G )x . In the same manner, Eq. (2) reads σ ,x + iκσ +f e i(κx−ωt) = −iωpe i(κx−ωt) , Once again, we expand the periodic parts to A truncation to a finite number of plane waves N is carried out for implementational purposes, say by −s ≤ m ≤ s 1 , such that N = 2s + 1. We multiply those sums by e −iG x , and integrate the result over the unit cell. Due to the orthogonality property of the Fourier series, the only non-vanishing terms are those satisfying the condition G = G. We can write for each G where We assemble all the equations of each G into a single matrix system, namely, where f A is a column vector of 2N components that assembles all the Fourier coefficients off G andq G ; h A is a column vector of 3N components, which assembles all the Fourier coefficients of σ G ,Ď G andp G ; and D A is a 3N × 2N matrix that is composed from 3 diagonal and 3 zero matrices. For further details the reader is referred to Appendix B. In the same manner, h A is expressed by assembling the matrix equations of the constitutive relation in Eq. (10), namely, where L A is a 3N × 3N matrix that contains the Fourier components of the material properties. Similarly to D A , J A is a 3N × 2N matrix that is composed from 3 diagonal and 3 zero matrices; w A is a column vector of 2N coefficients that assembles the Fourier components of the electric potential the displacement; and m A is a column vector of 3N components that contains the Fourier coefficients of the eigenstrain. Extracting the average fields -that are associated with G = 0 -from Eq. (16) yieldsh where J 0 and L 0 are the parts of the matrices that multiply the average fields and (•) s denotes the reduced matrix or vector without the G = 0 terms 2 . Our next step is to express the terms w s usingw andm. To this end, we first write Eq. (15) in the form Next, we rewrite Eq. (18) by separating the equations that do not include G = 0 and obtain where Q s is a square matrix whose entries multiply the average fields, and T contains the terms of the matrix D T L which multiply the average fields. Next, the fluctuating terms of w can be expressed using its average value and the material properties, such that Substituting Eq. (20) into Eq. (17) yields We can now identify the effective properties that act on the macroscopic kinematic fields, i.e., on The effective properties (22) are of the form reported in Refs. [34,35,37]: they are nonlocal in space and time, i.e., functions of κ and ω; and include two pairs of adjoint properties-the Willis couplingsS,S † and the electromomentum couplingsW ,W † -which are absent in the microscopic relations. In this work, we focus on maximizing the electromomentum coupling at the long wavelength limit, namely at κ = 0, in which caseS =S † andW =W † . We further restrict attention to the low frequency range, where the electromomentum coefficient is expected to depend linearly on the frequency [34,36]; Accordingly, we optimizeW at 2 kHz. (Indeed, our calculations show that in the range 0 -10 kHz,W grows linearly with the frequency.) Henceforth, the matrix in Eq. (22) will be denoted by K.
In what follows, we analyze K, with the understanding that the component that is relevant to this work is K 23 .

Discrete material choice optimization
We begin by presenting a discrete material choice approach, where the base materials that compose the periodic structure are chosen from a predefined set of available materials given in Table 1. We focus on a laminate made of three constituents, optimize it over the choice of the materials, and evaluate simultaneously the length of each layer that leads to the optimal design, while the total length of the periodic cell is fixed to 3 mm. Thanks to the use of standard base materials, the results of this design parameterization are directly applicable and can contribute significantly to practical implementation. As the optimization is based on discrete variables and the number of design variables is small, a genetic algorithm (GA) is used for solving the optimization problem. The design variables are assigned to the choice of material in each of the three material layers, from a library of candidate materials, and to the length of each material phase. This corresponds to the maximization problem maximize y j ,l 1 ,l 2 |W |, subject to y j ∈ [1, 2, ...10], j = 1, 2, 3, where y j denotes the pointer to a discrete material that is a member of Table 1. For example, if y 2 = 3 then the second layer is Al 2 O 3 . The length of the third layer is computed explicitly as the complementary of the total length, hence, it is not considered as a design variable. Note that for each minimal (negative) value ofW , there is an equivalent maximal (positive) value ofW that can be designed by flipping the order of the structure, since the sign ofW depends on the coordinate system [36]. Other operations can also yield equivalent designs with the same absolute value ofW , as will be discussed in the next section. Hence, our optimization implementation was to minimize the objective without referring to the absolute value, namely, to obtain the most negative electromomentum coupling.  In general, when using optimization methods, the optimization process might end up with a local minimum instead of a global one, and the outcome can strongly depend on the initial starting point of the whole process. Hence, some form of verification is needed. We define a confidence parameter to determine the threshold for a satisfactory solution. The parameter is computed as follows where multiple runs have been considered; R best is the number of runs that resulted in the optimal solution; R total is the total number of runs; E avg and E best are the average and optimal electromomentum coupling out of all the outcomes, respectively. The minimum confidence required for ending the search for the optimal design depends on the number of runs, as presented in Table 2, together with the calculated confidence of the process.
Our results show that the optimal structure is 2.052 mm PMN-PT Single Crystal X2B-0.106 mm PMMA-0.842 mm Soft PZT HK1HD, (Figure 2), yieldingW = 6.312 · 10 −7 Cs/m 3 , which is one order of magnitude larger than the electromomentum coefficient of the reference laminate. We observe that the optimizer provided a laminate with high contrast between the properties of its constituents: PMMA has the smallest value of each one of the electromechanical properties; PMN-PT Single Crystal X2B has the largest piezoelectric coefficient and mass density; PZT HK1HD has the largest dielectric coefficient and second-largest mass density. The asymmetry is not only in the material properties but also in the geometry: the middle layer is very thin with respect to the other two layers. The above features reoccur in subsequent sections.

Gradient-based optimization
The efficiency of genetic algorithms strongly depends on the number of design variables, exhibiting a sharp decline when the number of design variables increases. Even though the one-dimensional problem that we address does not require a large number of design variables, the capability to optimize with a large number of variables is desirable in general; for example, in extensions of the current work to higher-dimensional problems. In such cases, gradient-based optimization is more suitable, and this section is dedicated to optimizations using continuous variables that is solved by such gradient-based algorithms. In its most general form, the optimization problem includes design variables that govern the material properties and the geometry of each of the three material phases. An essential building block is the derivation of the analytical sensitivity analysis of the objective function with respect to the design variables, which we carry out in the sequel.
The general problem consists of finding the optimal material properties and the optimal length of each layer. This leads to the following maximization problem maximize z mn ,l 1 ,l 2 |W |, subject to l 1 + l 2 ≤ l T ,  where z mn denotes the material property m, i.e., B, A, C, ρ, in this order, of layer n. For example, z 23 is the dielectric coefficient of the third layer. The notation z mL and z mU represents the lower and upper limits of property m. Each material property is optimized independently, thus relaxing the constraint on the relation between the various material properties. The purpose of this FMO [41] is twofold: first, it delivers a theoretical upper bound on the electromomentum coupling; second, it provides guidelines on how to compose a laminate that maximizesW .
In order to obtain results with some correlation to realistic material properties, the design variables are limited in-between the extreme properties among the real materials considered in Section 3, as presented in Table 3. The material properties were normalized with respect to the same limits. The normalized properties are computed as follows where χ m and χ phys are the normalized and the physical property, respectively; and χ L and χ U are the lower and upper limits, respectively.

Sensitivity analysis
We carry out a differentiation of the objective function with respect to each one of the design variables, i.e., the lengths of the first two layers and the material properties of each layer. Considering first the length variables, the derivative of Eq. (22) using the product rule yields Since the PWE expands the governing-and constitutive equations in Fourier series, the entries in each matrix listed in Eq. (27) are basically Fourier coefficients, as given in Eq. (8). Thus, differentiatingξ G and following the same procedure yields the derivatives that we are looking for, namely ∂ L s ∂ l n , ∂ Q s ∂ l n and ∂ T ∂ l n . Since we are analyzing trilayer laminates, the integral is divided into three parts, namely,ξ Then, following the product rule, the derivative is expressed as The chain rule yields and we recall that total length of the periodic cell is fixed, hence where (·) denotes the matrices L s , Q s and T. Now that we have the derivatives above, we can follow the same process as described in Section 2, and obtain the desired derivatives in the same manner as in computing L s , Q s and T. Note that according to Eq. (27), one must compute the derivative of Q −1 s , and this is carried out using the relation The sensitivities with respect to material properties are derived similarly to the procedure outlined above for the length variables. Analogously to Eq. (27), the sensitivity with respect to any material property is expressed as Again, the entries of the matrices are Fourier coefficients, therefore the derivative of Eq. (28) with respect to the material properties is Henceforth, a procedure analogous to the PWE is carried out to obtain the required sensitivities, namely ∂ L s ∂ z mn , ∂ Q s ∂ z mn and ∂ T ∂ z mn . Note that we have compared the evaluation of the analytical expressions for the derivatives with numerical derivatives that were obtained by finite differences. This comparison yielded an excellent agreement, thereby verifying the correctness of our sensitivity analysis.

Test case: optimization over the lengths
We begin our numerical investigation with a benchmark problem of length optimization; this simple problem allows us to graphically verify our optimization method.
We consider two triplets of materials. The first triplet consists of the materials PZT4, BaTiO 3 and Al 2 O 3 . The graphical solution is presented in Figure 3a, where we observe that the optimal design is located near the diagonal. This means that the third layer should be thin in order to maximize the absolute value ofW , which is found to be 8.462·10 −9 Cs/m 3 . Solving the problem by a numerical optimization scheme using the method of moving asymptotes [59], we obtain the same result with the iterative convergence as plotted in Figure 3b. The optimal solution is characterized by the material distribution of l 1 = 1.38 mm, l 2 = 1.28 mm and l 3 = 0.3979 mm.
The second triplet consists of the materials Al 2 O 3 , BaTiO 3 and PMMA, which corresponds to the reference laminate when all layers have the same length. Figure 3c shows the graphical solution, where we observe that the optimal design is located near the diagonal. The gradient-based optimization achieves the same result,W = 9.52 · 10 −9 Cs/m 3 , as indicated in Figure 3d, which shows the iterative convergence of the gradient-based optimization. This result reflects a fourfold improvement with respect to the reference laminate, achieved solely by changing its geometry. Specifically, this optimal design is the lengths l 1 = 0.69 mm, l 2 = 2.13 mm and l 3 = 0.19 mm, for layers one to three, respectively. Both results are summarized in Table 4. We observe that the optimal designs share a common feature: a thin layer with the lowest mass density and piezoelectric coefficient. This is a repeating feature throughout the following optimization problems.
The graphical solution in Figure 3c shows that the optimization problem with this set of materials exhibits several local minima. Using a gradient-based method, this might lead to results that  Table 4: Results of length optimization using gradient-based optimization. The optimal lengths of the layers in the first composition are 1.38 mm, 1.28 mm and 0.3979 mm, and for the second composition they are 0.69 mm, 2.13 mm and 0.19 mm. Each color represents a different material. Both optimal designs exhibit a structure with one relatively thin layer.
depend heavily on the initial structure. In order to overcome this problem, the optimization process has been divided into two steps. The first step seeks the optimal region by calculating the objective function at 20 random points inside the feasible region. Once the best starting point is found, in the second step we use the gradient-based optimization to steer the design towards a nearby minimum point.
As inferred from Eq. (25), the length of the third layer is not a design variable. This may hamper convergence because the maximum change in the length of the third layer is related to the design variables, which are the lengths of the two other layers. For example, in cases where the optimal solution is located near the diagonal, the length of the third layer is small compared to the first and the second layers. Then a design change that is suitable for the design variables, might be too liberal for the third layer and divergence may occur. Therefore, conservative values of a move limit on the change of the design variables were imposed. In all examples, a maximum move limit of 0.03 mm was enforced.

Optimization over different material properties
Having verified our gradient-based optimization using a test case, we proceed to investigate optimization problems whose design variables are chosen material properties, where the rest of the properties are set to the properties of Navy type II PZT-5H. (d) Convergence of the gradient-based optimization. In both contours the initial (optimal) design is marked with red cross (black dot), and the objective function is multiplied by 10 9 .
We begin by considering the design space of the piezoelectric coefficient and mass density of each layer, yielding six design variables. For this case, the optimal electromomentum coefficient is 5.687 · 10 −6 Cs/m 3 , which is three orders of magnitude greater than than that of the reference laminate. This value is achieved by four equivalent designs, as displayed in Figure 4; the properties were color-scaled as yellow (blue) color denotes the upper (lower) limit value. The same color scaling will be used in the subsequent results. The optimal design exhibits a large contrast between the material properties, as 5 out of 6 properties attain the limit values, and one remaining piezoelectric coefficient attains the intermediate value 96.25 C/m 2 . This result qualitatively agrees with the insights of Pernas-Salomón et al. [36], who carried out heuristic homogenization of a single trilayer piezoelectric element, based on its scattering properties. Specifically, Eq. (41f) therein suggests that the electromomentum coefficient depends on the contrast of B/A between the layers, with some weight that depends on the rest of the properties. Accordingly, it appears that the optimal laminates maximize the contrast between the two layers which have the largest weight function. The four optimal designs are similar to each other: the mapping from one optimal design to another is obtained by swapping the first and the third layer, and/or flipping the contrast of the mass density. This means that layers with minimal mass density are replaced with maximal mass density and vice versa; layers with maximal mass density are replaced with minimal mass density. We refer to this transformation as the flipping of ρ. Each of the two transformation steps alone flips the sign ofW .
We demonstrate these two transformations using the laminates in Figure 4a and in Figure 4d, that share the same arrangement of the piezoelectric coefficient, but flipped values of the mass density, leading to the same magnitude with opposite sign of the electromomentum coefficient. The laminates in Figure 4a and Figure Figure 5: Free material optimization over the dielectric and the mass density properties. Panels (a) and (b) show two equivalent optimal designs, where the properties were color-scaled as yellow (blue) color denotes the upper (lower) limit value, and the optimal |W | is 4.516 · 10 −6 Cs/m 3 .
values of ρ, leading to the same value ofW . These observations correspond to the direction-dependent response of the medium [36], captured by the anisotropic nature ofW . As noted by Pernas-Salomon et al. [36],W flips sign when the coordinate system is reversed, which is equivalent to the inversion of the unit cell.
Our second case study optimizes over the mass density and dielectric coefficient. We find four designs, two of which are displayed in Figure 5, which are related through the transformations mentioned earlier. The optimal coefficient is 4.516 · 10 −6 , which is smaller than the coefficient that was obtained when optimizing over ρ and B, yet it is still three order of magnitude larger than that of the reference laminate. Interestingly, here all six properties attain limiting values, not only five as in the optimization over ρ and B, i.e., the dependency ofW on A is different from its dependency on B.
The last optimization problem considers all these three material properties, i.e., A, B and ρ, as design variables. Figure 6 presents two optimal and equivalent designs of the design variables. Again, the mass density and dielectric coefficient of all the layers attain the limiting value, as well as the piezoelectric coefficient of two out of the three layers. The resultant electromomentum coefficient is larger by an additional order of magnitude, in comparison with the coefficient that was obtained when only pairs of material properties were optimized.
Finally, we report that additional computations (not shown here) with ρ or B/A as prescribed constant values in the unit cell, lead to zeroW , regardless of the design variables chosen. The fact that constant B/A yields zeroW is in agreement with the results of Pernas-Salomón et al. [36]. However, Eq. (41f) in the latter reference, which we recall is based on the homogenization of a single scatterer, implies thatW can be nonzero with constant ρ, a result that is different from our current calculations, which are based on homogenization of periodic media.  Figure 6: Free material optimization over the piezoelectric, the dielectric and the mass density properties. Panels (a) and (b) show two equivalent optimal designs, where the properties were color-scaled as yellow (blue) color denotes the upper (lower) limit value, and the optimal |W | is 1.0354 · 10 −5 Cs/m 3 . Figure 7a presents the hypothetical optimal design that is obtained by removing all design restrictions. Accordingly, the design variables are the four material properties and the length of the layers. The resultant electromomentum coefficient is 5.7·10 −5 Cs/m 3 , which is fourfold the coefficient that was obtained when the optimization was over the triplet ρ, A and B. All the material properties attain the limiting values, except the piezoelectric coefficient of one single layer. Interestingly, the optimal stiffness is constant in the cell, the value of which is the lower limit, corresponding to a soft laminate. The optimal lengths yield a thin layer (volume fraction of 11%), in-between two thick layers (volume fractions of 47% and 42%). This result resembles the one that we obtained in the length optimization problem. Note that, again, equivalent designs can be obtained by mirroring the periodic cell and/or flipping the mass density values. Such equivalent design is shown in Figure 7b. Interestingly, the optimization process not only increasesW , but also increases the relative contribution of both the electromomentum and the Willis effects to the effective linear momentum. To show this, we compare in Figure 7c the momentum ratiosS u ,x ρ u ,t (solid lines) andW φ ,x ρ u ,t (dashed lines), in the reference laminate (blue lines) and the FMO design-based laminate (black lines), as functions of the frequency. (The reader is referred to Appendix C in the Ref. [34] for details in the calculations of u ,t and φ ,x .) We also observe that whileW φ ,x in the reference laminate is negligible with respect toS u ,t , it surpassesS u ,t in the optimized laminate.

Complete free material optimization
To quantify the contribution of the geometrical variables to the optimization, we optimize next only over the four material properties, keeping the length of each layer fixed and equal to 1 mm. The resultant electromomentum coefficient is 4.192 · 10 −5 Cs/m 3 . This value is about 71% of the coefficient that was obtained when the geometrical and material parameters were included in the design space. Two of the equivalent, optimal designs are displayed in Figure 8. These laminates are different from those that were obtained so far, as none of the layers exhibits a maximal piezoelectric coefficient. show two equivalent optimal designs, where the properties were color-scaled as yellow (blue) color denotes the upper (lower) limit value, and the optimal |W | is 5.7·10 −5 Cs/m 3 . Panel (c) displays the momentum ratios of the off-diagonal couplings: Solid blue (black) curve represents the Willis momentum ratio, namely,S u ,x ρ u ,t for the reference laminate (FMO optimal design), and dashed blue (black) curve denotes the electromomentum coupling momentum ratio, namely,W φ ,x ρ u ,t for the reference laminate (FMO optimal design).

Optimization over the designable mass density and stiffness
The results of the previous section provide a hypothetical upper bound for the electromomentum coefficient, which is accessible only if all material properties can be engineered independently. The significant progress on the design of engineered stiffness and mass density [60,61], motivates us to term the optimal design of C and ρ as the applicable design. Figure 9 shows the solution of this optimization over C and ρ, when the rest of the properties are chosen according to the DMO optimal result (Figure 2). The resultant electromomentum coefficient is 1.282 · 10 −6 Cs/m 3 , which is only one order of magnitude lower than the hypothetical upper bound. Here again, the optimization process delivers a constant minimal stiffness, and maximizes the contrast of ρ in the cell, as it attains the lowest limit in layers two and three, and the upper limit in layer one.
A summary of the optimal electromomentum coefficients, together with the reference laminate and the improvement relatively to it, is given in Table 5.

Summary
A prominent challenge in engineering today is to control elastic waves; metamaterials exhibiting the electromomentum effect [34] offer a way to tune the phase velocity and generate asymmetry in the phase angle of such waves [34,36]. To make these phenomena more pronounced, it is required to design these metamaterials with large electromomentum coefficient. In this work, we have utilized discrete material-, topology-and free material optimization methods to maximize the electromomentum coupling of periodic piezoelectric trilayer laminates in the low frequency, longwavelength limit. These optimization methods are essentially inverse homogenization methods, since we have formulated the objective function, namely, the electromomentum coefficient, as the end result of a rigorous homogenization process [34,37]. The discrete material optimization, which uses a predefined set of real materials as its design space, has yielded an electromomentum coefficient that is two orders of magnitude larger than the arbitrary laminate that was first analyzed by Pernas-Salomón and Shmuel [34]. We observed that the optimization yields a choice of materials that tends to maximize the contrast between their electromechanical properties, in a unit cell that comprises one thin layer in-between two thick  and PMMA layers); the optimal laminates that were reported in this work, and the relative improvement with respect to the reference laminate coefficient.
layers. The free material optimization uses a design space with continuous variables, the values of which are bounded in-between the extreme values of the materials considered earlier. We have employed gradient-based algorithms that were implemented together with a sensitivity analysis of the objective function. Using these algorithms, we have solved optimization problems whose design variables are pairs of the mass density and one of the electromechanical properties: the dielectric-, piezoelectric-or stiffness coefficients. The rest of the properties were set to the properties of the reference laminate. We have found that the optimal values of each one of these pairs enlarge the electromomentum coefficient by another order of magnitude, in comparison with the optimal laminate that the discrete material optimization delivered. We observed that the pair that demonstrated the best improvement was with the piezoelectric coefficient. We further found that when all material properties are designable, the electromomentum coefficient is enlarged by another order of magnitude.
The optimal laminates reveal how the electromomentum depends on the contrast in the electromechanical properties between the layers. Specifically, we observed that the optimization tends to maximize the electromechanical properties in certain layers, and minimize the same properties in the remaining layers. Two exceptions are the piezoelectric coefficient, which, while tending to the two extreme values in two layers, tends to an intermediate value at the remaining layer; and the stiffness, which the optimization minimizes throughout the whole cell. The optimal laminates also show that the geometry that maximizes the electromomentum coefficient is of a thin layer, in-between two thick layers.
Collectively, our results provide guidelines for the design of metamaterials with maximum electromomentum coefficient, paving the way for future work on their integration in wave control ap-plications.

Appendices A The convergence of the PWE method
To examine the convergence of the method, the electromomentum coupling was evaluated for a combination of three fixed materials. Results show convergence at 30 Fourier waves, as shown in Figure 10. Moreover, to examine the rate of convergence, three different divisions were considered; one is equally divided, one with one short layer compared to the other two layers and the third is the average of the previous two. As expected, the higher the contrast between the length of each layer, the slower the rate of convergence of the method.

B Detailed expressions for the matrices in the PWE method
The matrices can be divided into sub-matrices following and as has been mentioned before, D A is composed of three diagonal matrices and three zero matrices, namely D The description of the vectors is similar, where h A is and w A w (1) and f A definitely and for m A m (1) =η G , m (2) = 0,m (3) = 0.