Dynamics and wave propagation in nonlinear piezoelectric metastructures

This paper reports dynamical effects in one-dimensional locally resonant piezoelectric metastructures leveraged by nonlinear electrical attachments featuring either combined quadratic and quartic, or essentially quartic potentials. The nonlinear electromechanical unit cell is built upon a linear host oscillator coupled to a nonlinear electrical circuit via piezoelectricity. Semi-analytical harmonic balance (HB)-based dispersion relations are derived to predict the location and edges of the nonlinear attenuation band. Numerical responses show that weakly and moderately nonlinear piezoelectric metastructures (NPMSs) promote a class of nonlinear attenuation band where a bandgap and a wave supratransmission band coexist, while also imparting nonlinear attenuation at the resonances around the underlying linear bandgap. Besides, strongly nonlinear regimes are shown to elicit broadband chaotic attenuation. Negative capacitance (NC)-based essentially cubic piezoelectric attachments are found to expand the aforementioned effects over a broader bandwidth. Excellent agreement is found between the predictions of the HB-based dispersion relations and the nonlinear transmissibility functions of undamped and weakly damped NPMSs at weakly and moderately nonlinear regimes, even in the presence of NC circuits. This research is expected to pave the way toward fully tunable smart periodic metastructures for vibration control via nonlinear piezoelectric attachments.


Introduction
Chains of oscillators featuring nonlinearities exhibit unusual dynamical phenomena that enable remarkable ways for manipulating mechanical waves [16,48,52,56]. For instance, the acoustic switches, acoustic diodes, and acoustic rectifiers [9,28,30,31,[39][40][41]43,73], which enable one-way mechanical wave transmission through breaking down the linear time-invariant Maxwell-Betti reciprocity theorem 1 , are among the most notable examples of such mechanical devices for wave manipulation. Lattices built upon uniform anharmonic potentials also exhibit discrete breathers that open possibilities for novel energy harvesting approaches [10,12,21,47,69,76]. Since the seminal studies of the Fermi-Pasta-Ulam-Tsingou nonlinear model [27], from which wave supratransmission [28,44,45,67,70,73] and the existence of solitons [21,55,71] were demonstrated, a significant body of research has focused on the transmission features of such model and representative setups [47]. Most research focuses on nano-and micro-scale nonlinear lattices where the bandgap is generated by wave reflection (phononic crystals). A few efforts have investigated the interplay between Bragg scattering and nonlinear local resonances in weakly nonlinear lattices [12,35,36].
Periodic structures with linear resonating attachments are able to provide subwavelength wave filtering [42]. Since the forbidden band onset is given by the natural frequency of the attachments, wave manipulation is possible at wavelengths much larger than the lattice size. This is opposed to the bandgap formation mechanism in phononic crystals, which does depend on the host geometry and constituent material properties, thus demanding large lattices to tackle long wavelengths. The linear resonance of undamped mechanical oscillators is characterized by a large response at its natural frequency, which steeply decreases as the excitation frequency becomes mistuned. This narrowband feature can be enhanced through the addition of nonlinearities such as the Duffing oscillator, whose dynamical behaviors and applications are currently well understood [11,17,23,34,72]. However, only a few studies are attempted to harness its broadband and high-level responses in establishing periodic lattices for vibration attenuation or isolation purposes.
Lazarov and Jensen [38] numerically investigated the wave filtering properties of periodic linear lattices combined with local Duffing-type oscillators. An analytical equation based on the harmonic balance method (HBM) accurately predicts the onset and bandwidth of the nonlinear attenuation band, which is shifted as a function of both the input magnitude and the nonlinear coefficient of the Duffing attachment. Rothos and Vakakis [58] numerically examined the wave propagation features in linear chains with local essentially nonlinear attachments. Chaotic motion and subharmonic periodic orbits were unveiled by means of a nonlinear integro-differential equation that captures the full dynamics of one single nonlinear unit cell and subsequent periodic lattices. Manimala and Sun [46] investigated the filtering properties of linear chains of oscilla-tors with Duffing-type resonating attachments, where the propagation and attenuation properties of chains featuring cubic hardening and cubic softening (bilinear softening) springs were determined through interrogating the mathematical models with mono-, bi-harmonic, and wave packet-type loads, at low and high amplitude levels. They reported that the attenuation bandwidth depends on both the input amplitude level and the level and sort of nonlinearity.
Fang et al. [26] explored nonlinear mechanisms that give rise to the chaotic band formation in chains of linear oscillators with local Duffing-type (with positive linear term) oscillators. They observed that in undamped chains with several unit cells undergoing strong nonlinear regimes, there exists a relationship among the multiple saddle-node bifurcations at the acoustic branch and the multiple period-doubling bifurcations at the optical branch, with the formation of a chaotic band and a broadband nonlinear attenuation zone. They also discussed the dispersion properties of such chains of oscillators relying on the homotopy method. Cveticanin et al. [18] generalized the response of a nonlinear mass-in-mass unit cell to periodic forces of any sort by assuming Ateb-type periodic solutions for the set of equations of motion. Moore et al. [49] scrutinized a nonlinear unit cell composed by a linear large-scale oscillator coupled to a small-scale oscillator in an essential (or nearly essential) strongly nonlinear manner. They showed theoretically and experimentally that irreversible nonlinear energy transfer happens from the large to the small scale when the unit cell undergoes impulsive loads, thus enabling local energy nonreciprocity. Campana et al. [13] analyzed the dispersion relationship of linear chains of oscillators with nonlinear resonators by using a perturbation approach and compared the results with those reported by Lazarov and Jensen [38].
Piezoelectric periodic structures have also attracted attention due to their flexibility in tuning the resonators parameters via analog or digital electronics [3,7,24,[61][62][63][64]. While the most research in piezoelectric metastructures demonstrates the effects of linear resonating attachments on the bandgap formation, only a few examples exploit the development of NPMSs. Zheng et al. [73] proposed a nonlinear piezoelectric beam to realize one-dimensional propagation based on wave supratransmission. The nonlinear setup consisted of a homogeneous slender beam featuring several independent piezoelements bonded onto the middle portion of the clamped-clamped beam. Some piezoelements were connected to linear resonating circuits, others to nonlinear bistable circuits, and the remaining ones were left at short-circuit (SC) condition, thus providing an asymmetric feature to the nonlinear metastructure. They proved numerically and experimentally that nonlinear wave supratransmission leads to one-directional wave propagation for harmonic excitations whose frequencies fall into the underlying linear bandgap, i.e., existing bandgap at linear or quasi-linear oscillating regimes. Other few studies in NPMS have focused on vibration reduction [1] and energy harvesting techniques [15], by using lattices with nonlinear, nonresonant dedicated electronics.
This paper investigates the dynamics and wave propagation features of one-dimensional NPMS by means of a lumped-parameter model under low-amplitude axial harmonic excitations. The mechanical host chain is modeled by using linear restoring and dissipative laws, whereas the local nonlinear attachments are first modeled in the form of electrical Duffing oscillators [57,60,74]. The effect of negative capacitance (NC) circuits on the dynamics and wave propagation in NPMSs is afterward approached, in an attempt to address the effects of (nearly) essentially cubic-type piezoelectric attachments [60,74,75]. The dynamic behaviors of the nonlinear piezoelectric unit cell (NPUC) are investigated by using frequency-and time-domain numerical methods. The wave propagation features of NPMSs are derived through semi-analytical HB-based dispersion relations that account for damping in both the nonlinear attachments and the host structure. The validation of the equations is given through direct comparisons against the numerical nonlinear transmissibility functions (NTRFs) of a chain of 20 NPUCs. This paper is organized as follows: Sect. 2 derives the NPUC equations of motion. Section 3 discusses the steady-state dynamics and bifurcations of the NPUC by means of the single-harmonic HBM and time-domain methods. In Sect. 4, the analytical equations that describe the band structure of conservative and nonconservative NPMSs are derived, based on complexdomain single-harmonic HB approximations, and the linear piezoelectric metastructure (LPMS) band structure is presented as the reference for the nonlinear cases. The effects of NCs and the mass-inductor relationships are illustrated, in the presence of damping in both the attachment and the main chain. Section 5 shows validations of the analytical equations by means of numerical (a) (b) Fig. 1 a Lumped-parameter model of a nonlinear piezoelectric unit cell (NPUC): mechanical host oscillator piezoelectrically coupled with a nonlinear resonating attachment (shunt circuit); b free-free chain of NPUCs RK integration of the set of equations for several chains of NPUCs with corresponding discussions. Section 6 includes concluding remarks and future directions in this research. Figure 1b shows an electromechanically coupled mechanical oscillator that is characterized by the mass m s , the viscous damping d s , and the stiffness k s . The parameters of the piezoelectric element are given by its mass m p , damping d p , and stiffness k sc , where sc stands for short circuit. The electrical behavior of the piezoelectric element is modeled as a controlled voltage source with output v(t), in parallel with the inherent piezoelectric capacitance, C p . The operation mode of the piezoelectric material is the 33-one 2 as the mechanical host mass is subjected to longitudinal harmonic excitations.

Equations of motion
The shunt circuit shown in Fig. 1b allows for different unit-cell behaviors depending on its setup: while a resonant linear piezoelectric unit cell (LPUC) can be obtained by setting N (ω) with L c and R c , nonlinear piezoelectric unit cells (NPUCs) can be obtained by setting N (ω) with L c , R c , and α, the latter one characterizing an electrical device with cubic response. This NPUC can be further rendered an essentially cubic one through the inclusion of NCs, characterized by parameter C n , as it will be discussed next.

Duffing-type piezoelectric attachments
Using a Newtonian approach for modeling the NPUC, the governing equation can be obtained as follows: where m = m s + m p is the lumped system inertia; d = d s + d p is the lumped system viscous damping; F(t) = Fe −iωt is an external harmonic forcing term with F being the peak amplitude, ω the angular frequency, t is the time; i ∈ C = √ −1 is the imaginary unit; the over-dot(s) denote(s) time derivative(s); and k = k s + k p is the lumped system stiffness to which the piezoelectric contribution is to be added next through the piezoelectric constitutive relationships. In practice, slender piezoelectric elements contribute little to the total system's inertia and damping; hence, m p and d p can be neglected without loss of generality.
The piezoelectric effect on the structural system is accounted for in the form of an internal force, f p , which is developed through the linear piezoelectric constitutive relationships as follows: where ε S is the piezoelectric material's permittivity, and: are, respectively, the mechanical stress, T p ; the modulus of elasticity measured at constant electric field, c E p ; the piezoelectric stiffness at short-circuit (SC) condition, k sc ; the distance between electrodes, l p ; the piezoelectric's surface area, A p ; the piezoelectric constant, e; the charge density per unity stress at constant electric field, d 33 ; the mechanical strain, S p ; the electrical displacement, D; the electric field, E; and the electrical voltage, v.
By substituting Eq. 3 into Eq. 2, and after some manipulations, the piezoelectric constitutive relationships can be recast into the following form: where θ = k sc d 33 is the electromechanical coupling term, and C p = ε S A p /l p is the inherent piezoelectric capacitance. By isolating v from Eq. 4, and by substituting the results into Eq. 1, the following set of equations of motion is obtained: where k oc = k sc + θ 2 C −1 p is the open-circuit (OC) piezoelectric stiffness. Equations 5 and 6 govern the behavior of the NPUC at OC condition.
The nonlinear electrical oscillator is next modeled by using the physical laws that govern the voltage in an inductor, a resistor, and in the electrical device with cubic response, as follows: where R c , L c , and α = 1/α are the resistance, inductance, and nonlinear electrical constant, respectively, and C = 0 is a null integration constant provided that the initial condition in the nonlinear cubic device is zero. The realization of the nonlinear electrical device for piezoelectric shunt applications is currently accomplished through either ferroelectric capacitors [74] or analog/digital circuits [57,60]. The equation that governs the voltage in the nonlinear shunt circuit shown in Fig. 1 is obtained by applying Kirchhoff's voltage law to Eq. 7 and by substituting the results into Eq. 6, which writes: Equation 8 can be combined with Eq. 5 and recast into matrix form to give: where: is known as the influence coefficient matrix [8].
The dimensionless version of Eq. 9 is: with the following dimensionless parameters: where τ = ω sc t is the dimensionless time and ω sc is the SC resonance frequency; x c = h P and q c = (C p F)/θ are the system's characteristic units with h being a prescribed displacement (static distance per unit force), so that u = x −1 c x and r = q −1 c q are the dimensionless dependent variables; and the prime (resp. double prime) denotes first (resp. double) derivative(s) w.r.t. τ . The term K 2 in Eq. 11 is the squared generalized electromechanical coupling coefficient, which can also be calculated by using the "modal" electromechanical coupling coefficient given by [53]: where ω 2 oc = (k s + k oc ) /m is the OC resonance frequency.
After applying F D nl in Eq. 11, the dimensionless electrical stiffness function, F nl k , writes: where K δ sc couples the mechanical and electrical domains, whereas δ 2 sc and α N are the linear and cubic electric terms, respectively. The inherent piezoelectric capacitance, given by δ sc (cf. Eq. 12), is the term that enables the realization of the piezoelectric tuned vibration damper [32], by setting the inductance in accordance to the following tuning rule: where ω t = ω e (cf. Eq. 12) is the oscillator's target (natural) frequency. When adding the nonlinear cubic term to the stiffness function, as in Eq. 14, a Duffingtype electrical oscillator becomes realized that features a linear positive restoring force term, which can be used for diverse piezoelectric shunt applications, e.g., energy harvesting [17,23,72].

Essentially cubic-type piezoelectric attachments via NC circuits
Negative capacitance (NC) circuits are known in piezoelectric applications for improving the electromechanical coupling. Moreover, the realization of essentially nonlinear piezoelectric attachments requires the use of NC circuits, as demonstrated recently [57,60,74,75]. As soon as F nl k in Eq. 14 becomes a nonlinearizable stiffness function in the presence of the cubic term, the nonlinear oscillator will have no preferential oscillating frequency anymore [14,74]. Otherwise, the attachment will exhibit a preferential oscillating frequency that will depend on the residual linear term, which is related to the resultant capacitance.
NC circuits should hold the same physical law that govern the voltage in a capacitor, even though NCs deliver a negative impedance. Hence, the voltage law: where C n denotes the impedance value delivered by the NC circuit, also holds for NCs. A convenient way of quantifying the extent to which the NC impedance matches the one of the piezoelectric capacitance is through defining the variable μ ∈ (−1, 0), such that: and C n = 0 when no NC circuit is present in the electrical oscillator [53].
The goal of canceling out the piezoelectric capacitance via NC circuits accounts to making μ → −1. However, several aspects regarding both the stability of the overall electromechanical system and the performance of the electrical components used to realize an NC circuit have to be taken into account when setting a NC according to Eq. 17. Further details about the interplay between NC circuits and structural dynamics can be verified in "Appendix A" and elsewhere [2,6,19,29,50,53,68].
By substituting Eq. 17 into 16, and by using Kirchhoff's voltage law with the resulting NC's voltage law and remaining voltage laws in Eq. 7, Eq. 8 can be rewritten as follows: where is the residual capacitance. A straightforward inspection of Eq. 19 reveals that μ = −1, otherwise, C r → ∞ which causes system instability.
Since the residual capacitance is modified by the addition of the NC circuit, the inductance can be recalculated by making obvious substitution in Eq. 15, so that a given target frequency previously calculated for the RL-α-shunted NPUC is kept for the RLCα-shunted (negative-C) NPUC. In this regard, it is expected that L cc << L c as C r >> C p , where L cc is the recalculated inductance of the electrical oscillator in the presence of a NC device.
With the addition of the NC circuit, the electrical oscillator's nonlinear stiffness function becomes accordingly modified as follows: which makes the set of equations in Eq. 9 to be recast in the following form: Accordingly, the set of dimensionless equations is obtained as follows: where γ = L cc /L c is the dimensionless inductance rate between the inductor tuned without NC circuit (using C p ) and the inductor tuned in the presence of a NC circuit (using C r ); and μ N = C p /C r is the dimensionless capacitance ratio between the piezoelectric capacitance and the residual capacitance. Equation 22 governs the motion of the NPUC in the presence of an NC circuit and the Duffing-type NPUC in Eq. 11 can be obtained from Eq. 22 by making γ = 1 and μ N = 1. From Eq. 22, it can be found that the dimensionless attachment's stiffness function becomes: from which the strength of nonlinearity, σ [11], can be calculated after comparing the cubic and linear coefficients, as follows: Equation 24 shows that the strength or degree of nonlinearity can be controlled by either the nonlinear coefficient or by the input amplitude level (cf. term α N in Eq. 12 which depends on q c = C p F/θ ), or by both, while inversely controlled by the relationship between the residual and inherent piezoelectric capacitances, μ N , and by δ sc that denotes the rate between the mechanical host's natural frequency and the one of the electrical attachment.

Dynamics of the NPUC under harmonic excitation
This section investigates the NPUC harmonic steadystate features. The equation of motion given in Eq. 22 is compactly rewritten as follows: where the over-bars denote coupled mass, stiffness, and damping matrices, F denotes the external force vector, and z = [u, r ] T denotes the generalized dimensionless coordinates of the electroelastic body. The NPUC harmonic steady-state response is studied next by assuming single-harmonic solutions of the form: where Ω = ω/ω sc is the dimensionless frequency with U c = Acosφ u , U s = Asinφ u , and R c = Bcosφ r , R s = Bsinφ r , being the amplitude components of the assumed host and attachment harmonic amplitudes, given by A and B, respectively, and relative phases given by φ u and φ r , respectively. The subscripts u and r are relative to the dimensionless generalized coordinates. By neglecting the terms that are not related to Eq. 26, the following implicit frequency-domain system of equations is obtained after balancing the harmonics: The codimension-1 bifurcations associated to the nonlinear system of Eq. 25 can be obtained by rearranging Eq. 27 as follows: and by solving for its discriminant set to zero [11].
In special, fold-type bifurcations of cycles, which are related to the jump phenomenon [11,65], are studied and tracked by means of a numerical scheme, aiming at establishing a relation between this phenomenon and the emergence of harmonic solutions within the underlying linear bandgap, as it will be shown later in detail.

Steady-state responses of the NPUC
Since no closed-form solutions for Z are available for Eq. 27, an arc-length continuation technique together with a HB numerical algorithm, as well as an alternating time-frequency (AFT) scheme to solve for the cubic forcing term [37,66], is used to investigate the NPUC steady-state response. The relevant codimension-1 bifurcations are calculated by solving Eq. 25 with a shooting method, from which the monodromy matrix is obtained and the corresponding Floquet multipliers calculated as a function of Ω [22]. The numerical stability of the branches is assessed by evaluating the complex-valued eigenvalues of the monodromy matrices at each excitation frequency [54]. Furthermore, Eq. 25 is solved using a Runge-Kutta (RK)-based time-domain integration scheme (MAT-LAB's function ode113) as a reference for the HBand shooting-based predictions. The time-domain integration is performed for 50 × 10 3 periods of the input excitation. The RMS amplitude value at steady state, divided by the input displacement h, is used to provide the transmissibility functions. Focus is given to the dynamical response around the in-phase and outof-phase modes, in the bandwidth given by Ω ∈ [0.94, 1.08] for the assumed parameters. The adopted parameters for all the numerical simulations presented onward are: h = 1.00 × 10 −3 ; λ 1 = 3.16 × 10 −3 ; λ 2 = 2.38 × 10 −5 ; K 2 33 = 4.36 × 10 −3 ; δ sc = 1.002. To assess the effect of the NC circuit, the   adopted μ N = 0.85, whereas μ N = 1.00 for the shunt circuit without NC. NTRFs are calculated for six values of α = 1.00 × 10 7 , 1.00 × 10 8 , 1.00 × 10 9 , 2.30 × 10 9 , 5.00 × 10 9 , and 1.00 × 10 10 , which, according to Eq. 13, yield into the following α N = 5.18 × 10 −6 , 5.18 × 10 −5 , 5.18 × 10 −4 , 1.19 × 10 −3 , 2.59 × 10 −3 , and 5.18 × 10 −3 . They are referred to as α N ,l onward, with l = 1, 2, ... 6. For the Duffing-type piezoelectric attachment (μ N = 1.00), it follows that γ = 1, whereas for the essentially cubic-type piezoelectric attachment (μ N = 0.85), γ = 0.85, which indicates that the inductance has been retuned by using the residual capacitance. By making use of the dimensional coefficients [11], the corresponding σ 's for each degree of nonlinearity with μ N = 1.  Figure 2 shows the harmonic steady-state responses of the Duffing-type (in plot (a)) and the essentially cubic-type (in plot (b)) NPUCs, for the set of the chosen nonlinear coefficients. Each plot illustrates the responses obtained from the HB and shooting methods, together with the results from the direct RK integration. The dotted blue and light blue lines indicate the LPUC (linear piezoelectric unit cell) responses at SC and OC conditions, respectively. Since the inductance in the shunt circuit has been tuned using the OC resonance frequency, a antiresonance arises at that frequency and shifts toward higher frequencies as the degree of nonlinearity increases. The two sharp peaks at both sides of the antiresonance represent the in-phase and the out-of-phase normal modes. The steady-state LPUC displacement response is calculated by making u(τ ) = U e −iΩτ and r (τ ) = Re −iΩτ in Eq. 22, which, after some manipulations, yields in: from which the LTRFs (linear transmissibility functions) are calculated by making T R = U (Ω)/P, whereas the steady-state attachment responses can be calculated by using: The response of the weakest NPUC in Fig. 2 (α N ,1 ) matches the behavior of the LPUC calculated with Eq. 29. For this nonlinear case, the results from both the frequency-domain HB and shooting methods (in dashdotted lines) match the ones from the RK integration (in continuous lines). Since neither bifurcations nor unstable branches were detected, the weakest NPUC can be regarded as a quasi-linear NPUC.
The next degree of nonlinearity, α N ,2 , exhibits noticeable distortions at the two normal modes, which cause attenuation of the linear resonances. The response of this weakly NPUC exhibits jump phenomenon and bifurcations (in markers), as well as unstable responses (in thin black lines connecting the markers) between bifurcations. Four fold-type bifurcations have been detected, wherein the first two ones are located around the in-phase normal mode, and the remaining ones around the out-of-phase normal mode. The steady-state responses are in excellent agreement among each other (HBM, SM, and RK methods).
The next three nonlinear cases represent moderately nonlinear unit cells, illustrated in Fig. 2 (α N ,3 , α N ,4 , and α N ,5 ). The results from HBM and SM (frequencydomain methods) present some discrepancies with respect to RK integration. While the steady-state amplitude values seem properly predicted by using either HB or shooting methods, in relation to the time-domain integration, the loci of the resonance jumps and bifurcations show small inaccuracies that increase as the degree of nonlinearity increases. Even so, the three numerical methods agree in showing large attenuations at the in-phase and the out-of-phase normal modes and also in estimating the shift of the antiresonance toward higher frequencies. Besides the occurrence of two fold-type bifurcations close to the in-phase normal mode, one Neimark-Sacker bifurcation has been detected close to the out-of-phase normal mode, the dynamics of which should be taken into account when exciting the NPUC with nonstationary loads, e.g., frequency sweeps.
The last nonlinear case, α N ,6 in blue lines, corresponds to a strongly NPUC that exhibits the largest amplitude attenuations around the normal modes. The attenuation trends are similar as the observed ones at the previous weakly and moderately nonlinear regimes, yet the predictions of the loci of the resonance jumps and bifurcations are increasingly inaccurate. This nonlinear case also exhibits the largest antiresonance shift toward high frequencies.
With regard to the essentially (cubic-type) piezoelectric attachment shown in Fig. 2b, the foremost effect is that of enlarging the bandwidth at which both the linear and nonlinear effects arise. The weakly essentially cubic-type NPUC exhibits its nonlinear normal modes further apart from each other, in comparison with the Duffing-type NPUC. Besides, similar dynamical behaviors observed for the Duffing-type piezoelectric attachment are verified for the essentially cubictype one in the steady-state plots: prominent nonlinear distortions, shifts of the jump resonances (bifurcations) toward high frequencies, and shift of the antiresonance toward high frequencies. Figure 2 also presents the long-term phase planes as insets, which evidence the boundness of both the host mass displacement and the attachment electrical charge, and also the relative phase relations between both sub-domains.

Steady-state responses of a 5-NPUC NPMS
Before to give more insights into the reliability of the single-harmonic HB approximation in estimating the NPMS' harmonic steady-state responses, Eq. 22 is generalized to multiple unit cells, as follows: where I ∈ 1 N ×N is the identity matrix, with N being the number of unit cells, are the normalized damping and short-circuit stiffness matrices, respectively; and their time derivatives are N × 1-sized vectors; P ∈ R N ×1 is the host input vector, and 0 ∈ 0 N ×1 represents the null electrical charge vector. The NTRFs of a 5-NPUC NPMS are calculated next by using the HB algorithms, and the results are compared to the RK schemes. The computational cost of one single shooting method's iteration renders the calculation of both the steady-state responses and bifurcations of the 5-NPUC NPMS unfeasible. The 5-NPUC NPMS is built upon the repetition of the previously analyzed NPUC along one single direction, cf. Fig. 1a. Harmonic displacements are prescribed for the NPUC mass at the left end of the system (with the same h = 1 × 10 −3 ), the responses are calculated at the host mass of the unit cell at the right end, and the NTRFs are calculated w.r.t. h.
The NTRFs in Fig. 3 show that the HB-based numerical continuation scheme is able to follow all the complicated branches of the 5-NPUC NPMS, which unveils bandwidths with several multiple solutions. It can then be expected that the complexity of the branches increases as the number of NPUCs increases. Moreover, the cost of assessing the stability of the frequencydomain solutions along the multiple branches becomes computationally comparable to the time-domain integration of the equations.
In contrast to the complicated HB-based frequencydomain branches illustrated in Fig. 3 (thin, dashed lines), the trends unveiled by the time-domain RK scheme (thick, continuous lines) resolve in a clear way the steady-state NPMS responses under the aforementioned input excitation. While the undamped LPMS (linear piezoelectric metastructure) exhibits a locally resonant bandgap around Ω = 1 (L or LC curves in Fig. 3), yet showing new resonances at the bandgap edges, the weakly NPMSs (in Fig. 3a, b) show noticeable nonlinear attenuations of the linear resonances at the bandgap edges. Moreover, the moderately and strongly NPMSs (in Fig. 3c, d) present a wide nonlinear attenuation band, characterized by strong attenuations at the linear resonances, while the moderately nonlinear metastructure (light-blue curves) still exhibits a small bandgap-like attenuation band located around the right linear bandgap edge. Figure 3 evidences that the nonlinear attenuations at the linear resonances around the bandgap become stronger as the nonlinear coefficients become increasingly larger. At its turn, this behavior causes the emergence of harmonic solutions within the forbidden band, cf. α N ,2 and α N ,3 cases in Fig. 3a, b, which start spanning from the left to the right bandgap edges, as a function of the nonlinear coefficient level as well. This effect becomes more prominent in case α N ,4 (Fig. 3c,  d), where the linear bandgap is mostly overlapped with such waves, whereas the case α N ,5 shows no bandgaplike region at all, but a wider attenuation zone that flattens the steady-state responses at and around the underlying linear bandgap. Similar behaviors were already unveiled in the NPUC's steady-state responses shown in Fig. 2, where the fold-type bifurcation that arises in the vicinity of the in-phase normal mode becomes increasingly shifted to high frequencies as the level of nonlinearity increases, causing large attenuations at the resonances, and small antiresonance shifts toward high frequencies. Then, the bifurcation shifts toward high frequencies that happens at each unit cell in the NPMS can be linked to the emergence of such harmonic solutions within the otherwise (linear) forbidden band. This essentially nonlinear phenomenon is characterized as wave supratransmission [70] that suddenly underpins wave propagation within the forbidden band, as a function of either the input amplitude level or nonlinear attachment coefficient.
More prominent effects due to the presence of the nonlinearities in the electrical attachments, e.g., wave supratransmission and localized modes will be discussed in Sect. 5 in the light of the dispersion relations presented in the following section.

Wave propagation in NPMSs
Wave transmission properties in NPMSs are investigated next by using the HB technique. The NPUC in Fig. 1 is periodically repeated an infinite number of times, as ideally shown in Fig. 1a, and the j-th NPUC is considered for the mathematical developments.
The Floquet-Bloch theorem that governs the propagation of plane waves in a periodic medium [33]: is used as a solution for the set of equations of motion, where x j = a j = j denotes the position of the j-th mass in the NPMS with a = 1 the dimensionless lattice constant; κ = aυ is the wave propagation constant of harmonic wave [38], with dimensionless frequency Ω, which is equal to the wavenumber,υ, multiplied by the lattice constant; and the spatial part of Eq. 32 is given by: where C is an arbitrary amplitude value. By assuming that the j-th NPUC is at the 0-th position in the NPMS, the displacements of the j = −1 and j = 1 masses become: Different from the cases in Sect. 3, where either a force or prescribed displacement is applied at the first host mass, the study of the infinite chain of oscillators is tackled by assuming that the j-th unit cell is an autonomous system, whose free response is to be approximated in the frequency domain using HB techniques.

Dispersion relation for LPMSs with undamped host chain
The following set of equations of motion governs the dynamics of the j-th undamped LPUC: The mechanical equation of Eq. 35 is rewritten by using Eq. 34, which yields in: The temporal part of the solution of Eq. 36 is sought in the form: which is substituted into the electrical equation of Eq. 35 to give: By substituting Eq. 38 into Eq. 36, equating to zero the coefficients of e −iΩτ and rearranging terms, the equation that governs the dispersion relation for the LPMS with undamped host chain is obtained: Figure 4 illustrates some band structures of the LPUC calculated through Eq. 39, for a set of μ N . For the conservative LPUC, no mechanical wave with excitation frequency equal to or close to the frequency at which the propagation branch becomes divided into acoustic and optical branches is expected to propagate along the host chain. The bandgap onset is located at the frequency at which Re(κ) = π . The bandgap end is located at the frequency at which Re(κ) = 0. In between the bandgap onset and end, only harmonic evanescent wave solutions exist, i.e., Re(κ) = 0 and Im(κ) > 0. The relative phase of the displacement of any two contiguous unit cells in the chain will increase until Ω = 2, at which Re(κ) = π and a Bragg-type bandgap starts.
As previously stated in Sect. 2.2, the inductance can be recalculated as a function of μ N . As μ N → 0, γ → 0, once that the residual capacitance is increasingly larger (cf. Eqs. 15, 17, and 19). Moreover, the coupling between the structural and electrical domains improves as μ N → 0, which explains the noticeable broadening of the bandgap limits in Fig. 4a as the inductance becomes smaller. On the other hand, if the inductance is not recalculated as a function of μ N , the target frequency becomes modified (ω t → 0) due to the changes in the residual capacitance as μ N → 0. In the latter case, the bandgap becomes shifted toward low frequencies and broader as μ N → 0, cf. Fig. 4b, which again evidences that NCs improve the coupling between the structural and electrical domains.

Dispersion relation for NPMSs with undamped host chain
The cubic term is now added to the electrical domain of Eq. 35, which writes: The temporal part of the solution of Eq. 35 is now sought in the form [38]: where is a dimensionless parameter that shows the order of the amplitude of the motion, superscript * denotes complex conjugate, R is the amplitude of the lth harmonic of the series with l = 1, 3, ..., and Ω is the fundamental dimensionless frequency of the assumed motion. Equation 41 is substituted into Eq. 40 by assuming that the contributions of the harmonics of order equal to and higher than l = 3 are neglected, which gives: By substituting Eq. 42 into Eq. 40, equating to zero the coefficients of e −iΩτ and rearranging terms, the dispersion relation for the NPMS with undamped host chain is obtained as: which is similar to Eq. 39, except from term −3α N R 2 1,0 that makes κ dependent on both the nonlinear coefficient and the steady-state attachment amplitude, given by Re(R) in Eq. 30. Figure 5 illustrates some band structures of the NPMS, for a set of α N values. Similar to the linear case, the resonating attachments split the longitudinal propagation branch into acoustic and optical branches. However, the bandgap edges shift toward high frequencies as the degree of nonlinearity increases. This is a well-known effect of cubic nonlinearities in mechanical metastructures featuring positive cubic restoring forces (hardening springs) [13,26,35,38]. The same shift effects are observed at both NPUCs with and without NCs in Fig. 5a, b, respectively. Although the bandgaps are also expected to shift toward high frequencies as the level of the input magnitude increases, the linear host model would not predict accurate responses in such a situation.
It is worth noticing in Fig. 5 that the strongly nonlinear case (in blue continuous lines) presents bandgap edges well ahead of the cut-on frequency of the LPMS. The remaining nonlinear cases present bandgap onsets that fall into the frequency range of the LPMS bandgap. By taking as a reference the onset of the LPMS bandgap, the nonlinear band structures suggest that wave propagation occur within the forbidden band of the LPMS. The band at which such phenomenon exists starts from the onset of the underlying LPMS and ends at the onset of the nonlinear bandgap. This phenomenon is known as wave supratransmission, which will be more clearly evidenced through the NTRFs in Sect. 5. Figure 6 shows the effects of damping in the attachments on the NPMS band structures. Two cases have been chosen, namely λ 2 = 2.38 × 10 −5 that can be regarded as an undamped NPMS, and λ 2 = 1.19 × 10 −4 , i.e., a weakly damped NPMS. The presence of dissipation mechanisms causes the acoustic and optical branches to merge into one single propagation branch, in which the truly forbidden band is no longer available, but an attenuation band that presents similar behaviors as the underlying bandgap. This is especially verified for the first weakly nonlinear case, in which the band structures in Fig. 6a, b are quite similar to those shown in Fig. 5a, b. Moreover, the case with weakly damped attachments shows noticeable changes in both the real and imaginary parts of κ. Although there are no evanescent wave solutions in the attenuation band (Re(κ) = 0, Im(κ) > 0), strongly attenuated oscillating harmonic wave solutions are present. While the attenuation band shifts toward high frequencies due to the nonlinearity, the small damping still has the capability of broadening the attenuation bands observed in the undamped linear and nonlinear band structures.

Dispersion relation for LPMSs with damped host chain
When linear viscous damping is added to the LPMS, the first line of Eq. 35 should be accordingly modified as follows: By substituting Eq. 38 into Eq. 44, equating to zero the coefficients in front of e −iΩτ , and rearranging terms, the following equation is reached for the dispersion relation of the LPMS with damped host chain: where and

Dispersion relation for NPMSs with damped host chain
The dispersion relation that involves the damped NPUC parameters is reached by substituting Eq. 42 into Eq. 44, which yields: where and Equation 48 is similar to Eq. 45, except from nonlinear terms that make κ dependent on both the nonlin-ear coefficient and the steady-state attachment amplitude, given by Re(R) in Eq. 30. While the results from Eq. 45 for the damped-host LPMS are interesting for applications in the linear regime, Eq. 48 constitutes the main result of this paper that links all the structural and nonlinear electric parameters of piezoelectric metastructures to the wave propagation constant of harmonic waves.

Bandgap edges
As observed in Figs. 4 and 5, the bandgap onset is given by Re(κ) = π , whereas the bandgap end is given by Re(κ) = 0. Hence, by taking the real part of Eq. 48, equating to −1 and solving for Ω, a set of four solutions can be found as follows: two of which are positive-valued and the remaining ones are negative-valued, i.e., meaningfulness. While one of the positive roots gives the resonant bandgap onset, the other one provides the onset of the Braggtype bandgap around Ω = 2z, cf. Fig. 4. Moreover, by taking the real part of Eq. 48, equating to 1 and solving for Ω, a set of four solutions is found as follows: two of which are real-valued solutions and the remaining ones are purely imaginary. One of the real roots is positive and gives the resonant bandgap end, whereas the another one is negative, so that it can be regarded as meaningfulness.
It is worth noticing that, according to Eq. 51, the bandgap onset is dependent on the simultaneous presence of damping at both the host and the attachment, whereas, according to Eq. 52, the bandgap end is completely independent of such mechanisms. The bandgap onset for the NPMS with undamped host (Sect. 4.2) can be obtained after making λ 1 = 0 in Eq. 51, while Eq. 52 still holds for the bandgap end. The bandgap edges for the LPMS (Sect. 4.1) can be obtained after making α N = 0 in Eqs. 51 and 52, which also cancels out their dependence on the amplitude. LPMSs with either undamped or damped host chains can be eventually obtained by keeping a sufficiently small λ 2 so that parametric analyses with λ 1 may be verified.

Numerical results
This section validates the theoretical dispersion relations derived in Sect. 4 for NPMSs through the results from time integration of the equations of motion. Sufficiently small damping coefficients i.e., λ 1 = 3.16 × 10 −3 and λ 2 = 2.38 × 10 −5 are considered so that the nonlinear forces are the dominant mechanisms in the dynamics and wave propagation of the investigated NPMSs. This also facilitates the attainment of the steady-state responses in affordable computation times. A second simulation scenario using λ 2 = 1.19 × 10 −4 is also investigated to determine the effects of weak attachment damping on the dynamics and wave propagation of a 20-NPUC NPMS. A Runge-Kutta scheme (MATLAB's function ode113 with stringent relative and absolute error criteria set at 1.0 × 10 −9 ) was used to perform the numerical integration for 50 × 10 3 periods of each harmonic prescribed displacement imposed at the first mass of the host. The frequency resolution used to obtain the steady-state responses within the bandwidth of interest, Ω ∈ [0.9875, 1.025], is 5.0 × 10 −5 ω SC , thus giving 7.71 × 10 2 different harmonic prescribed displacements. The dimensionless displacement, velocity, electrical charge, and electrical current are measured at the last unit cell, for each discrete frequency in the bandwidth of interest. These quantities are also recorded at all the NPUCs for some frequencies at which natural resonances of the underlying LPMS located around the linear attenuation band occur, to get more insights into the steadystate features along the chain. All the quantities are given as root-mean-square (RMS) values, calculated by using the last quarter of the corresponding time histories. Figure 7a, b shows, respectively, the steady-state responses of the 20-NPUC NPMS (cf. Figure 1a) shows dispersion relations considering μ N = 1.00 for a set of nonlinear coefficients. Figure 7c, d shows the same analyses considering μ N = 0.85 for the same set of nonlinear coefficients. Starting from the underlying linear attenuation band (in dashed lines), the figure evidences that the prediction of the relevant attenuation band features, such as its location, width, and limits, is in excellent agreement with those features directly observed in the NTRFs for low to moderate nonlinear coefficients. Table 1 presents the numerical predictions made with Eqs. 51 and 52 for the bandgap left and right edges, respectively, and corresponding bandwidths for cases α N = 0 (LPMS) and α N ,1...5 (weakly and moderately nonlinear metastructures) with μ N = 1.00. Figure 7c, d also evidences that the predictions of the dispersion relations for the attenuation band characteristics in the presence of NC circuits are as accurate as for the case without NC circuits. Table 2 presents the numerical predictions made with Eqs. 51 and 52 for the bandgap left and right edges, respectively, and corresponding bandwidths for cases α N = 0 (LPMS) and α N ,1...5 with μ N = 0.85. Moreover, the accuracy of the dispersion relations is not appropriate for the strongly nonlinear regime (α N ,6 ). This is due to the inherent limitation of the HBM in predicting quasiperiodic and chaotic steady-state behaviors that are leveraged by the NPMS at such a strong nonlinear regime, which will be unveiled and discussed in Sect. 5.2. The predictions for the bandgap edges at the strongest nonlinear regime are not shown in the Tables.

Dynamics and wave transmission in undamped NPMSs
The steady-state amplitude responses shown in Fig. 7a, c also evidence the existence of harmonic wave solutions that overlap a portion of, or entirely, the underlying linear attenuation band. As the degree of nonlinearity increases, a bandwidth composed of highamplitude harmonic waves emerges and starts overlapping the linear attenuation band from its left edge until the first frequency predicted by the dispersion relations at which Re {κ} → π (bandgap onset), for a given α N . The existence of waves within the frequency range of the underlying linear bandgap suggests that wave supratransmission is supported by the NPMS, notably at moderately and strongly nonlinear regimes. For the NPMS composed of a linear host chain and either combined quadratic and quartic or essentially quartic local potentials, wave supratransmission can be explained upon the nonlinear distortions that happen in the acoustic branch, which are more expressive than the ones in the nonlinear mechanical metamaterials studied by, e.g., Lazarov and Jensen [38] or Fang et al. [26]. Thus, nonlinear wave supratransmission in NPMSs has its underlying mechanism on the various saddle-node bifurcations at the acoustic branch, which depart toward high frequencies as a consequence of the increase of the nonlinearity. Figure 7a, c evidences that wave supratransmission at weakly and moderately nonlinear regimes does not completely overlap the underlying linear bandgap, but the amplitude function does suddenly drop off at the frequency for which Re {κ} ≈ π holds, still falling into the linear bandgap range. This frequency indicates the onset of the nonlinear bandgap for weakly and moderately nonlinear regimes. Hence, weakly and moderately NPMSs exhibit a frequency band in which both a wave supratransmission band and a bandgap coexist. Moreover, the insets in the NTRF plots in Fig. 7a, c evidence that the right edge of the nonlinear attenuation band shifts to the right with increasing the nonlinear coefficient, which brings large reductions for the linear resonances located at that edge. As the level of nonlinearity becomes high (α N ,6 ), the underlying linear bandgap becomes completely dominated by wave supratransmission. The bandwidth at which this phenomenon manifests is significantly broader than the one at which the underlying linear bandgap exists. The attenuation band in strongly NPMS is seen to be dominated by wave supratransmission, which provides broadband nonlinear attenuation such that all the linear resonances of the underlying LPMS become significantly reduced.
A localized mode for α N ,5 (moderately nonlinear regime) at Ω ≈ 1.004 has been detected, which can be better observed through the insets of Fig. 7a, c. This localized mode is evidenced for both cases of NC, i.e., for μ N = 1.00 and μ N = 0.85, in the form of a sharp peak within the forbidden band that exceeds the transmissibility reference value 3 (dotted light-blue curves). This localized mode may be leveraged by the relatively high degree of anharmonicity in the lattice [10,55,59,69,76]. Although out of the scope of this paper, a comprehensive study on the emergence, features, and applications of such discrete breathers is a relevant topic.

Dynamics and wave transmission in weakly
damped NPMSs Figure 8 compares the NTRFs of the NPMS endowed with weakly damped attachments with the results of the developed dispersion relations. The trends observed in this figure are qualitatively similar to those shown in Fig. 7, although the overall steady-state amplitude trends become smoother. The sudden amplitude dropoffs at the onset of the attenuation bands become smoother as well, although they happen at the same frequency as in the case with undamped attachments. These observations are mostly noticed at moderately and strongly nonlinear regimes. Wave supratransmission is still present in the NTRFs, and the degree of agreement between the dispersion relations and the NTRFs is still high. Moreover, the localized mode observed at Ω ≈ 1.004 when α N ,5 for both cases of NC setup is strongly attenuated due to damping in the attachment. Damping mechanisms in LPMS are regarded as a means to enhance the attenuation bandwidth, and this fact can be verified for the weakly and moderately NPMSs, as the width of the attenuation bands becomes slightly broader compared to those observed in Fig. 7. Tables 1 and 2 also evidence this fact along the BW columns of the second rows, wherein more significant digits are available, yet the improvement on the bandwidth enlargement is quite small.

Analysis of the unveiled behaviors
The dispersion relations developed for NPMSs are highly accurate in predicting the bandgap behavior at weakly and moderately nonlinear regimes. The onset of the nonlinear attenuation band exactly matches the fre- quency predicted by the dispersion relations at which Re {κ} → π . Taking as a reference the onset of the underlying linear bandgap and the onset of a given nonlinear bandgap, the supratransmission frequency range could also be predicted for weakly and moderately nonlinear regimes. Moreover, the prediction of the cut-on frequency of the nonlinear bandgap degrades as the degree of nonlinearity increases. Recall from Sect. 3 that the nonlinear phenomena that happen at the in-phase and the out-of-phase nonlinear normal modes produce two different sorts of dynamical behaviors. While quasi-periodic and/or chaotic behavior are expected for moderately and strongly nonlinear regimes at the out-of-phase nonlinear normal mode, nonlinear attenuation and shifts of the antiresonance toward high frequencies given by the successive foldtype bifurcation to the right are expected at the inphase nonlinear normal mode. Then, the loss of performance in predicting the cut-on frequency is given by the inability of the HBM in predicting quasi-periodic and/or chaotic behaviors, predominantly happening at the out-of-phase nonlinear normal mode. Figure 9 illustrates some time histories of the last NPUC displacement. The time lapse is equivalent to ten periods the fundamental prescribed harmonic displacement at four selected resonant frequencies (Ω = 0.99, 1.00, 1.005, and 1.02) of the underlying LPMS. Each time history is plotted together with its corre- Table 1 Bandgap edges and bandwidths for case   sponding long-term phase plane (displacement vs. electrical charge of the last attachment), which allows verifying the existence of long-term quasi-periodic and chaotic responses. While for Ω = 0.99 and Ω = 1.02 (linear resonances located on both sides of the linear bandgap) only nonlinear attenuation is verified, highamplitude broadband strongly chaotic responses are observed at the innermost frequencies Ω = 1.00 and Ω = 1.005. Moreover, while the outermost frequencies show well-defined host-attachment phase relations, i.e., in-phase time histories for Ω = 0.99 and out-of-phase time histories for Ω = 1.02, at all the nonlinear regimes, the time histories of the innermost frequencies evidence that the host-attachment relations become quite complicated, yet the displacements and electrical charges remain bounded during uptime. In special, the time-domain trends for Ω = 1.00 clearly reveal the prominent wave supratransmission feature, since the displacement dramatically increases with increasing degree of nonlinearity. Another intriguing observation is that, for the case α N ,5 (red continuous lines), the displacement seems to have suddenly changed its phase, relative to the observed one in the remaining nonlinear regimes. The time histories for Ω = 1.005 show similar trends as the observed ones at Ω = 1.02, regarding both the displacement attenuations as the nonlinear coefficient increases, and the out-of-phase host-attachment behavior. However, for the strongly nonlinear regime, the host-attachment relationship becomes chaotic. Similar observations can be extracted from the time histories of the last NPUC featuring the NC with μ N = 0.85, wherein the phase planes are quite similar to those of the last NPUC without NC circuits.
The quasi-periodic and chaotic behaviors are further confirmed through the spectra of the displacements, illustrated in Fig. 10. While the spectra of the outermost frequencies, Ω = 0.99 and Ω = 1.02, verify a harmonic response mainly centered at the excitation frequency, the observed spectra for the innermost frequencies, Ω = 1.00 and Ω = 1.005, show harmonic responses at several discrete frequencies, even for the weakly nonlinear regimes, and quite broadband responses for the moderately and strongly nonlinear regimes. The spectra also unveil the existence of suband superharmonic responses of very low amplitude (small insets in Fig. 10). It is worth remarking that, even when the moderately and strongly nonlinear regimes elucidate broadband responses with a quite simple har-monic input, the highest peak in the spectra is still localized around the excitation frequency. Then, the broadband feature of such responses does not invalidate the methodology used in this research of obtaining the frequency-domain NRTFs based on taking the RMS value of the last quarter of the time histories as the steady-state response for a prescribed harmonic displacement. Figures 11 and 12 illustrate the steady-state amplitudes of each mass along the undamped and the weakly damped 20-NPUC NPMSs, respectively, for a set of nonlinear coefficients. The excitation frequencies are the same as the ones visualized in the time-and frequency-domain plots of Figs. 9 and 10. While slight changes are seen at Ω = 0.99 and Ω = 1.02 along the host chain with increasing the nonlinear coefficient (plots at the first and last columns of Figs. 11 and 12), dramatic qualitative changes are observed for Ω = 1.00 and Ω = 1.005 cases. Focusing on the Ω = 1.00 case (second column of Figs. 11 and 12), while the weakly nonlinear regimes still cause the NPMS to show the typical well-behaved exponential decrease of the displacements along the chain (quasi-linear behavior), moderately nonlinear regimes start leveraging irregular steady-state behaviors along the chain that cause the sudden increase of the displacement of the last host mass. The strongly nonlinear regime promotes steady-state amplitudes slightly greater than the observed ones at moderately nonlinear regimes. This fact evidences in a different manner that both the undamped and the weakly damped NPMSs do support nonlinear wave supratransmission. Regarding the Ω = 1.005 case (third column of Figs. 11 and 12), well-defined steady-state exponential decreasing behavior is observed along the host, even for weakly and moderately nonlinear regimes (quasi-linear behavior). The steady-state amplitudes are reduced as the degree of nonlinearity increases, yet the strongest nonlinear regime produces dramatic changes similar to those observed for Ω = 1.00, wherein the steady-state displacements are no longer well behaved.

Conclusions
This paper investigated the dynamics and wave propagation in nonlinear piezoelectric metastructures (NPMSs). The state of NPMSs undergoing low-amplitude harmonic prescribed displacements was modeled through   chaotic responses that were not properly predicted by the HBM-based procedures. Moreover, the NTRFs unveiled wave supratransmission, which is due to the pronounced shifts of the saddle-node bifurcations at the acoustic branch that overlap and dominate the frequency range of the underlying linear bandgap, when at moderately and strongly nonlinear regimes. The presented formulation could also predict the frequency band at which wave supratransmission manifests, as a function of both the nonlinear coefficient and the onset of the underlying linear bandgap.
The numerical results revealed that NPMSs present four relevant effects leveraged by the cubic term in the electrical circuit of the piezoelectric attachments, which could be harnessed for vibration attenuation and mechanical isolation purposes, as follows: (1) ondemand reconfiguration of the attenuation band as its onset shifts to high frequencies with increasing the nonlinear coefficient, which could be set in real-world applications through either analog or digital electronics; (2) large nonlinear reductions provided at the resonances located at both sides of the linear bandgap, due to the nonlinear jump resonance phenomenon; (3) wave supratransmission due to the shift of the multiple saddle-node bifurcations happening at the acoustic branch toward high frequencies, so that the NPMS becomes underpinning stable harmonic solutions within the frequency range of the underlying linear bandgap; and (4) chaotic attenuations at the resonances located at the right side of the underlying linear bandgap, promoted by the strongest nonlinear regimes. The presence of NC circuits to realize nearly essentially nonlinear attachments was seen to extend the NPMS attenuation band and underlying linear bandgap, thus enhancing all the aforementioned positive effects for vibration attenuation and mechanical isolation purposes. A discrete breather was also detected when the undamped NPMS was subjected to a relatively high degree of anharmonicity, regardless of whether NC circuits are or are not present in the shunt circuit.
In summary, this paper numerically demonstrated that piezoelectric lattices coupled to well-designed Duffing-type and essentially cubic-type circuits leverage relevant dynamical effects and unusual wave propagation properties that bring substantial benefits for vibration reduction and isolation applications. Future directions in this research include the numerical treatment of realistic distributed-parameter, one-and twodimensional NPMSs, for which the present study establishes a foundation to the understanding of the vast plethora of nonlinear phenomena and unveiled wave propagation features.
The stability limit has been determined for the LPUC and several LPMSs, formed by 5, 10, 15, and 20 LPUCs, with the parameters reported in Sect bility. Yet, this does not imply in that the twofold NC effects in reducing the inherent piezoelectric capacitance while increasing the electromechanical coupling are degraded in some way, as it will be shown next. Figure 13 shows the LPUC's and the LPMS' transmissibilities (LPMSs formed by 5, 10, 15, and 20 LPUCs), wherein the main plot illustrates the denser amount of resonances that arise as the number of cells increases, whereas the inset evidences the location of the first resonance at almost Ω = 0, placed at that frequency by the effect of the previously found μ N values close to the linear stability bound. Any attempt to improve μ N above the obtained stability bounds renders the linear electroelastic system unstable.
Another form of gaining insights into the stability of the linear electroelastic system under the effect of the S-type NC circuit is by iteratively calculating the generalized electromechanical coupling coefficient, K 2 e f f , by using: where ω 1 refers to the first eigenvalue of the generalized eigenvalue problem that is formed with matrices M and K without NC (cf. Eq. 22 with μ N = 1 and γ = 1), and ω NC,1 refers to the first eigenvalue of the generalized eigenvalue problem that is formed with the same matrices, yet with 0 < μ N < 1. Equation 57 can be recursively calculated as μ → −1, which shows that, as ω NC → 0 + , K 2 e f f → ∞. Moreover, as ω NC becomes a complex-valued eigenvalue due to μ N values that overcome the stability limit given by either Eq. 55 or 56, the negative real part of ω NC makes K 2 e f f to take negative, meaningfulness values [53]. Figure 14 illustrates the trends of K 2 e f f as a function of μ N , for the LPUC and a set of LPMS formed with 5, 10, 15, and 20 LPUCs. The figure has been set so as to show K 2 e f f values between the K 2 e f f value that corresponds to the electromechanical coupling without NC and K 2 e f f = 2.5, but even larger values aretheoretically-possible. Figure 14 also shows that a certain K 2 e f f value can be attained by either the LPUC or a LPMS with several LPUCs, by properly fixing the μ N value of all the LPUCs in the LPMS. That is to say, the same generalized electromechanical coupling and inherent piezoelectric capacitance reduction effects can be attained by one or several sequentially connected LPUCs, by properly setting all the μ N 's to the same value, yet the numerical value should be increased toward zero as the number of LPUCs increases. This explains the small μ N value that was chosen for the LPMS with 20 unit cells shunted with negative capacitances, although it produces the same net effect as μ N = 0.01 for the LPUC case, i.e., K 2 e f f = 1.0. From this discussion, it becomes straightforward to see that the μ N value that makes K 2 e f f = 1.0 for the LPMS with 20 LPUCs will cause a rather small effect in the electroelastic system consisting of one single LPUC, cf. Fig. 2, but the obtainment of the LPMS effects shown in Fig. 7 follows from repeating 20 times the investigated LPUC with μ N = 0.85. Figure 14 evidences that all the simulations in this paper have considered linearly stable electroelastic systems, both for the NPUC and the 20-NPUC NPMS cases. The inclusion of α N with the values reported in Sect. 3.2 was not shown to change the system's stability. An analytical expression and bounds for α N that ensure the electroelastic system's stability are left for future works.

Appendix B
In uniform nonlinear metamaterials 4 with cubic restoring forces [25,35,38,51], the saddle-node bifurcation that emerges and shifts to the right side in the spectrum as a function of the input amplitude has been associated with the onset of transmission within the underlying linear bandgap [70]. It is through this essentially nonlinear mechanism that wave propagation becomes enabled within the otherwise forbidden band. As the periodic solutions lose their stability through the several saddlenode bifurcations of the unit cells in the metastructure, the essentially nonlinear wave supratransmission phenomenon becomes enabled in the host structure [70].
Wave supratransmission is the essentially nonlinear phenomenon that allows the metamaterial host to underpin waves within the forbidden band [69,70]. The most research has focused on nonlinear hosts that are able to respond to amplitude-level variations. On the one hand, the versatility of nonlinear mechanical host structures is found to be limited, from the point of view of designing, realizing, and controlling the degree of their nonlinearity, which makes the input amplitude variations to be a better approach in practical terms. On the other hand, piezoelectric metastructures based on the nonlinear electrical circuit shown in Fig. 1b can achieve a great level of versatility, since their coefficients are synthesized in the electrical domain with some analog circuits and passive components, whose parameters are comparatively easier to tailor. In this line, the input excitation given for the linear mechanical host can be kept at the same level, whereas the nonlinear coefficients are varied by adjusting the electrical circuit's parameters [57,60]. According to Eq. 24, this approach is also a means for exciting and controlling the overall system's level of nonlinearity, which has a direct relationship with enabling and enhancing the nonlinear wave supratransmission phenomenon across the NPMS.
From the results shown in Fig. 7 for a set of α N values, five discrete frequencies that fall into the forbidden band have been chosen to demonstrate the nonlinear wave supratransmission phenomenon, as shown in Fig. 15. In special, the frequencies located at Ω = 1.002, Ω = 1.003, and Ω = 1.004 that fall into the forbidden band show similar trends to those found in the specialized literature concerning amplitude-induced wave supratransmission, cf. right column of Fig. 15. The frequencies Ω = 1.001 and Ω = 1.004, which are located at the left and right bandgap edges, respectively, also show the same behavior of rising their amplitudes and reaching a stable value, which is similar to those observed at the remaining frequencies, yet the initial amplitude of the latter frequencies is lesser than the one of the former frequencies because they are located at the bandgap edges.
The plots at the left column of Fig. 15 also help visualize the behaviors inside the bandgap and at its edges, for the same set of α coefficients defined in Sect. 3.1.