Brillouin-zone characterization of piezoelectric material intrinsic energy-harvesting availability

Vibration energy harvesting is an emerging technology that enables electric power generation using piezoelectric devices. The prevailing approach for characterization of the energy-harvesting capacity in these devices is to consider a finite structure operating under forced vibration conditions. Here, we present an alternative framework whereby the intrinsic energy-harvesting characteristics are formally quantified independent of the forcing and the structure size. In doing so, we consider the notion of a piezoelectric material rather than a finite piezoelectric structure. As an example, we consider a suspended piezoelectric phononic crystal to which we apply Bloch’s theorem and formally quantify the energy-harvesting characteristics within the span of the unit cell’s Brillouin zone (BZ). In the absence of shunted piezoelectric circuits, the wavenumber-dependent dissipation of the phononic crystal is calculated and shown to increase, as expected, with the level of prescribed damping. With the inclusion of the piezoelectric elements, the wavenumber-dependent dissipation rises by an amount proportional to the energy available for harvest which upon integration over the BZ and summing over all branches yields a quantity representative of the net available energy for harvesting. We investigate both monoatomic and diatomic phononic crystals and piezoelectric elements with and without an inductor. The paper concludes with a parametric design study yielding optimal piezoelectric element properties in terms of the proposed intrinsic energy-harvesting availability measure.


Introduction
Energy harvesting [1,2] has been a topic of intense and fruitful research for the past few decades. The notion of capturing or harvesting different types of energies, such as kinetic (vibrational), thermal, and electromagnetic, in nature or from man-made constructs has led to a variety of modern techniques of power generation, and this endeavor continues to grow vigorously among the research continuity. The advent of low-power energy-harvesting systems that can be incorporated into devices that are usually battery operated has brought about new opportunities for the design of energy-efficient products. Among a diverse range of energy-harvesting techniques, vibration-based energy harvesting has received considerable attention from researchers across a variety of disciplines [3]. Energy can be harvested from sources of mechanical vibration by means of piezoelectrics [4][5][6], shape memory alloys [7], ionic polymer metal composites [8,9], magnetostrictive materials [10], and structural instabilities [11]. Piezoelectric energy harvesting in particular has shown practical promise for the generation of low-power electricity for a variety of applications, including micro electrical devices [12], wireless sensors [13], transducers [14], microelectromechanical systems (MEMS) portable devices [15], structural health monitoring devices [16][17][18], biomedical implants [19], among others. Sodano et al [20] and Anton et al [4] presented in-depth reviews of research on vibrational-power harvesting using piezoelectric material-based stuctures.
A considerable amount of theoretical and experimental research has been carried out on harvesting energy from the mechanical vibrations of cantilever beams [2,21]. In this context, different configurations have been considered such as a beam with a tip mass and a piezoelectric/piezoceramic substrate bonded to its surface. Adhikari et al [22], for example, performed a study on the design and analysis of beam-based harvesters under random excitation and based on three extreme response (peak voltage) traits, namely, fractional time spent above a specific level, level crossing, and response peaks above a specific level. Energy harvesting in beams with a tip mass have also been examined in the nonlinear vibrations regime [23]. An extensive review article on the "role of nonlinearities in vibratory energy harvesting" was published by Daqaq et al [24].
In parallel to energy-harvesting research, the study of wave propagation in artificially structured materials has been an explosive area of research over the past three decades [25][26][27][28]. In this domain, two key classes of engineered materials emerged, namely phononic crystals (PnCs) [29][30][31][32][33][34][35] and locally resonant acoustic/elastic metamaterials [36][37][38][39][40]. The band-gap phenomenon and wave propagation characteristics in reciprocal lattice space [41] have been among the key features of interest in both these classes of materials. An important trait in particular is the spatial and temporal attenuation properties. If the interest is in free waves of a damped model, the temporal attenuation is quantified by the wavenumberdependent damping ratio corresponding to each dispersion branch [42,43]. Regardless of the level of damping, and for both free and driven waves, spatial attenuation exists in phononic materials and is quantified by the imaginary part of the wavenumber spectrum [43].
The merger of the two aforementioned fields-phononic crystal/metamaterial-based energy harvesting-has emerged as a natural outcome. Incorporation of piezoelectricity into phononic crystals for active tuning or energy harvesting has been the focus of numerous studies [44][45][46][47][48][49][50][51][52][53][54]. Similarly, elastic metamaterials have been linked with energy harvesting as surveyed in the review by Chen et al [55]. An attractive feature in this context is the ability to use piezoelectricity to actively tune the local resonance properties for attenuation [56], waveguiding [57], or energy harvesting [58]. Simultaneous vibration suppression and energy harvesting may also be realized, as demonstrated by Hu et al [59,60]. A more recent study by Dwivedi et al [61] also explored simultaneous vibration attenuation and energy harvesting in an elastic metamaterial. The metamaterial employed exhibited negative stiffness and consisted of an energy-harvesting 'smart material' embedded within the resonating units. The results of a parametric investigation followed and suggested that a piezo-embedded negative stiffness metamaterial boasts enhanced performance in the context of vibration attenuation and energy harvesting.
The theoretical and experimental studies of harnessing energy and generating power efficiently from artificial structures based on phononic crystals and acoustic/elastic metamaterials continue to receive significant attention and carry an immense amount of scope and future potential. While some of the research mentioned above has examined the dispersion characteristics considering a combination of elastic and piezoelectric properties [1,2,4,20], when it comes to analysis of the actual energy-harvesting capacity this has been done at the structural level where the focus is on the structural dynamics of the systems presented. In this paper, we present a new, formal approach for the characterization of energy-harvesting capacity that is fundamentally at the material level, rather than the structural or device level. We consider a phononic crystal with integrated piezoelectric elements and model it as a damped medium within a Bloch wave propagation analysis framework [42,43]. We obtain the set of wavenumber-dependent damping ratio curveswhich is a rigorous measure of dissipation capacity-and compare it directly with the corresponding set for the same phononic crystal without the piezoelectric elements. The dissipation curves for the latter, the non-piezo medium, give a direct indication of the 'raw' dissipation which represents unutilized/lost energy. On the other hand, the difference in the dissipation curves between the two systemsthe phononic crystal with piezoelectric elements versus the phononic crystal without piezoelectric elements-provides a formal wavenumber-dependent representation of the amount of energy available for harvesting. This Brillouin zone (BZ)based approach for energy-harvesting characterization directly follows the characterization framework for the concept of 'metadamping' [40,62], where the wavenumber-dependent dissipation of two statically equivalent waveguides are directly compared to produce a formal measure of either enhanced [40,[62][63][64] or reduced [63,64] dissipation (or damping capacity).
The rest of the paper is divided as follows. In section 2, an example of a piezoelectric phononic crystal (piezo PnC) modeled from a materials perspective is presented conceptually, considering both monoatomic and diatomic configurations. In section 3, we present the application of Bloch's theorem to a monoatomic PnC and a diatomic PnC. The application of Bloch's theorem to piezo PnCs is then followed in section 4. Wave propagation analysis and the calculation of our representation of energy-harvesting availability is detailed in section 5. A non-dimensional parametric study is provided in section 6, followed by concluding remarks in section 7.

Examples of piezoelectric phononic crystals modeled as materials with intrinsic properties
In this section, we present a physical conceptualization of two configurations of piezoelectric phononic crystals featuring monoatomic and diatomic periodicities; see figures 1(a) and (b), respectively. These configurations are selected for their suitability for future experimental realization; however, they serve only as example configurations. The monoatomic configuration comprises a periodic arrangement of equally sized masses, each suspended by a rectangular beam. Each of the rectangular beams is equipped with a unimorph piezoelectric patch consisting of a shunt circuit with or without an inductor; this patch performs the energy-harvesting function. The piezoelectric patch shown in figure 1 can use polyvinylderylene fluoride (PVDF) films, which are known to have a high d31 coefficient [1]. Each mass has one degree of freedom. The masses are interconnected by springs and viscous damping elements (dashpots). The diatomic piezoelectric phononic crystal configuration comprises a periodic arrangement of alternating small and large masses, each also suspended by a rectangular beam. The rest of the set-up is similar to the monoatomic configuration. In analysing each case, we apply a generalized form of Bloch's theorem on a single representative unit cell to accommodate complex frequencies [42,43] and then obtain the wavenumber-dependent dispersion and damping ratio diagrams for a specific set of parameters. For simplicity, we will utilize lumped parameter mass-springdashpot models; see section 3 for details on the models and the application of Bloch's theorem.

Monoatomic phononic crystal
The monoatomic PnC without energy-harvesting patches is represented by the schematic shown in figure 2(a). In the figure, a is the distance between two consecutive unit cells, m denotes the mass, K (··· ) and C (··· ) denote spring stiffness and damping coefficient, respectively, and n is used to identify the unit cell under consideration, i.e. n + 1 and n − 1 identify the unit cells to the right and to the left, respectively. The time-dependent displacements of the masses (identified by the superscript index) are denoted by x n (t), x n + 1 (t), and x n − 1 (t). The governing equation of motion for mass m in the nth unit cell is written as: where a dot denotes differentiation with respect to the time.
We impose a plane-wave Bloch solution [65] as: wherex κ is the wave amplitude, r is the mass position, κ is the wavenumber, t is the time, j = √ −1, and n + g is used to identify a particular unit cell. For the unit cell under consideration, the index g = 0 and for the unit cells towards the left and right of the central unit cell, g =−1 and g =+1, respectively. Substituting equation (2) into equation (1) yields a homogeneous equation forx κ , which is written as: Equation (3) can be converted to a first-order problem through a state-space transformation [65][66][67] in the following manner: where the system matrices and the state vector are: , , For equation (4), we assume a solution of the formŶ κ = Y λ κ e λt , whereŶ λ κ is a complex amplitude state-space vector corresponding to eigenvalue λ. The damped frequency band structure can now be obtained by applying this solution form and solving the resulting eigenvalue problem given by: Solving equation (6) and using the identities e jκa + e −jκa = 2 cos κa and 1 − cos κa = 2 sin 2 κa 2 yields a second-order equation (dispersion relation) in terms of λ, which has two complex roots appearing as a complex conjugate pair. The complex solution is expressed as: where s refers to the mode or branch number indicating the number of complex conjugate pairs in the solution; for the monoatomic model, s = 1. From this expression, we can obtain the wavenumber-dependent resonant frequency ω rs (κ), which is defined as: the wavenumber-dependent damped frequency: as well as the wavenumber-dependent damping ratio: The phase velocity and the group velocity corresponding to a particular mode can also be obtained using: respectively. The second-order equation for this model in terms of λ can be non-dimensionalized by utilizing a specific set of parameters. For this purpose, we choose the lumped model parameters of the system. We define and substitute K 1 = K, K 2 = γ 1 K, C 1 = C, C 2 = γ 2 C, and λ = jω, where the branch index s has been dropped for brevity and ω is introduced to denote a complex frequency as the general solution of a characteristic equation corresponding to the state-space eigenvalue problem of equation (6). We now define the following quantities: and where Ω denotes non-dimensional damped frequency, ω n is a characteristic frequency defined in terms of K and m, and ζ n denotes damping ratio. Dividing by ω n 2 yields the characteristic equation Equation (13) has two complex roots appearing as a complex conjugate pair. The coefficients of this non-dimensional equation are found in the appendix. In section 5 where we present the main results, we will resort again to the dimensional form in defining the parameters; however, the nondimensional version provides a more general presentation of the characteristic equation. We will follow this interchangeable approach of use of either dimensional or non-dimensional notation where appropriate throughout the paper.

Diatomic phononic crystal
The diatomic configuration without energy-harvesting patches are similarly represented, using the schematic shown in figure 2(b). The governing equations of motion for masses m 1 and m 2 in the nth unit cell are written as: As in the monoatomic model, we impose a plane-wave Bloch solution: where l = 1, 2 is an index. Substituting equation (16) into equations (14), and (15). yields two homogeneous equations forx κ1 andx κ2 , which can be written in a matrix form as: where,X Converting equation (17) to a first-order problem through a state-space transformation results in equation (4), where the system matrices and the state vector are now: The damped frequency and damping ratio band structures can now be obtained by applying a solution of the form Y κ =Ŷ λ κ e λt and solving the associated eigenvalue problem given by equation (6). Solving equation (6) for the current model yields a fourth-order equation in terms of λ, which has four complex roots appearing as two complex conjugate pairs. Given that we have a diatomic unit-cell model, we obtain a physical root corresponding to the acoustic branch (s = 1) and a physical root corresponding to the optical branch (s = 2).
In order to non-dimensionalize the present model's fourthorder equation in terms of λ, we choose m 1 and the lumped model parameters representing the beam attached to it. We define and substitute m 1 = m, m 2 = γ 1 m, and λ = jω. Dividing by ω n 4 and using equation (12) yields the characteristic equation: Equation (20) has four complex roots appearing as two complex conjugate pairs. The coefficients of the non-dimensional equation are found in the appendix.

Piezoelectric phononic crystals
In this section, we now consider monoatomic and diatomic phononic crystals without and with an inductor and apply unitcell Bloch analysis to the displacement and voltage fields.

Monoatomic phononic crystals with piezoelectric elements
4.1.1. Piezo PnC without inductor. The monoatomic configuration consisting of shunted piezoelectric patches where the shunt circuit does not consist of an inductor is represented using the schematic shown in figure 3(a). In figure 3(a), R l denotes the resistance of the shunt circuit. The mechanical and electrical governing equations with electromechanical coupling [22], for the n th unit cell, are written as: where v is the voltage generated by the shunted piezoelectric patch, θ is the electromechanical coupling, and C p is the capacitance of the shunt circuit.
We impose a plane-wave Bloch solution for the voltage: whereν κ is the voltage wave amplitude. Substituting equations (2) and (23) into equations (21) and (22) yields two homogeneous equations forx κ andν κ , respectively, which can be written in matrix form in the following manner: Converting equation (24) to a first-order problem by a statespace transformation results in equation (4), for which the system matrices and the state vector are: Similar to the non-piezo monoatomic model, the damped frequency and damping ratio band structures are obtained by applying a solution of the formŶ κ =Ŷ λ κ e λt and solving the associated eigenvalue problem given by equation (6). From equation (6), we obtain a third-order equation in terms of λ, which has a real root and two complex roots appearing as a complex conjugate pair, hence this particular case exhibits only one branch.
In order to non-dimensionalize the present model's thirdorder equation in terms of λ, we choose the lumped model parameters of the beam and introduce certain non-dimensional parameters incorporating the parameters of the energy harvester. We define and substitute K 1 = K, K 2 = γ 1 K, C 1 = C, C 2 = γ 2 C, and λ = jω. Dividing by ω n 3 and using equation (12) and, and where α is a non-dimensional time constant and k 2 is an electromechanical coupling coefficient, yields the characteristic equation: Equation (27) has a real root and two complex roots appearing as a complex conjugate pair. The coefficients of the nondimensional equation are found in the appendix.

Piezo PnC with inductor.
The monoatomic configuration consisting of shunted piezoelectric patches where the shunt circuit consists of an inductor may also be represented using the schematic shown in figure 3(a). In figure 3(a), L denotes the inductance of the shunt circuit. The mechanical governing equation with electromechanical coupling remains the same as in the previous case, but the electrical governing equation with electromechanical coupling, for the nth unit cell now becomes: Substituting equations (2) and (23) into equations (21) and (28) yields two homogeneous equations forx κ andν κ , respectively, which can be written in matrix form in the following manner: where, Converting equation (30) to a first-order problem, again by a state-space transformation, results in equation (4) where the system matrices and the state vector are: The damped frequency and damping ratio diagrams are now obtained by considering a solution of the formŶ κ =Ŷ λ κ e λt and solving the associated eigenvalue problem given by equation (6). From equation (6), for this problem we obtain a fourth-order equation in terms of λ, which, depending on the choice of parameters, either has two complex roots appearing as a complex conjugate pair and two real roots or four complex roots appearing as two complex conjugate pairs. For the purpose of obtaining wave dispersion characteristics, the branch or mode that always gives a complex conjugate pair must be considered.
Non-dimensionalization of this model's fourth-order equation in terms of λ is done by using an approach similar to the previous cases. We first define and substitute K 1 = K, K 2 = γ 1 K, C 1 = C, C 2 = γ 2 C, and λ = jω. Dividing by ω n 4 , and using equations (12) and (26) and where β is a non-dimensional constant, yields the characteristic equation: Equation (33) has four complex roots appearing as two complex conjugate pairs, however one of them must be ignored as it gives non-dimensional damped frequencies equal to zero for κ = [0, π]. The coefficients of equation (33) are provided in the appendix.

Diatomic phononic crystals with piezoelectric elements
In this section, we consider diatomic piezolectric phononic crystal configurations.

Piezo PnC without inductor.
A diatomic configuration consisting of shunted piezoelectric patches where the shunt circuit does not consist of an inductor are represented using the schematic shown in figure 3(b). The mechanical and electrical governing equations with electromechanical coupling for the nth unit cell are written as: Following the same procedure as done earlier, we impose a plane-wave Bloch solution for the voltage, which is written as: where l = 1, 2 is an index as defined earlier. Substituting equations (16) and (38) into equations (34)-(37) yields four homogeneous equations forx κ1 ,x κ2 ,ν κ1 andν κ2 , which can be written in matrix form in the following manner: where, Using state-space transformation, we convert equations (39) and (40) to a first-order problem resulting in equation (4), where the system matrices and the state vector are: The damped frequency and damping ratio band structures are then obtained by applying a solution of formŶ κ =Ŷ λ κ e λt and solving the associated eigenvalue problem given by equation (6). Solving equation (6), for this particular case, yields a sixthorder equation in terms of λ, which has two real roots and four complex roots appearing as two complex conjugate pairs.
In order to non-dimensionalize the sixth-order equation of the present model in terms of λ, we choose m 1 , the lumped model parameters of the beam attached to it, and the electrical parameters of the shunted piezoelectric patch bonded to that beam. We define and substitute m 1 = m, m 2 = γ 1 m, 10 R l , and λ = jω. Dividing by ω n 6 and using equations (12) and (26) produces the characteristic equation: Equation (43) has two real roots and four complex roots appearing as two complex conjugate pairs. The coefficients of equation (43) are available in the appendix.

Piezo PnC with inductor.
Finally, the diatomic configuration consisting of shunted piezoelectric patches where the shunt circuit consists of an inductor may as well be represented using the schematic shown in figure 3(b). The mechanical governing equations with electromechanical coupling remain the same as in the previous case, but the electrical governing equations with electromechanical coupling for the nth unit cell now become: Substituting equations (16) and (38) into equations (34), (35), (44), and (45) yields four homogeneous equations for x κ1 ,x κ2 ,ν κ1 andν κ2 , which can be written as two equations in matrix form. The first equation is the same as equation (39) and the second equation is: where, and M, C, K, T, C p and R l are the same as in the previous case. Upon state-space transformation, equations (39) and (46) are converted to a first-order problem resulting in equation (4), where the system matrices and the state vector are: The damped frequency and damping ratio diagram are then obtained by applying the solution formŶ κ =Ŷ λ κ e λt and solving the associated eigenvalue problem given by equation (6). Upon solving equation (6), we obtain for this case an eighthorder equation in terms of λ, which, depending on the choice of parameters, either has four complex roots appearing as two complex conjugate pairs and four real roots or eight complex roots appearing as four complex conjugate pairs. For the purpose of obtaining wave dispersion characteristics, the branches or modes that always generate complex conjugate pairs must be considered.
Equation (49) has eight complex roots appearing as four complex conjugate pairs; however, two of them must be ignored since they are equal to zero for κ = [0, π].

Wave propagation, raw dissipation, and energy-harvesting availability of piezoelectric phononic crystals
In this section, we examine all configurations, namely the PnC, piezo PnC without inductor, and piezo PnC with inductor.
For each case, we first obtain the damped dispersion and damping ratio diagrams. We then subtract the wavenumberdependent damping ratio quantities of the non-piezo PnC from those of the corresponding piezo PnCs to quantify the energyharvesting availability for each model. All analysis is done for both the monoatomic and diatomic models.

Monoatomic configurations
In this section, we examine the monoatomic class of models and compare the frequency and damping ratio band structures for the PnC, piezo PnC without inductor, and piezo PnC with inductor, the parameters used to generate the plots are given in table 1. The parameters have been taken from a previous study [68] which investigates vibrational energy harvesting using a cantilever beam with a tip mass in a flow field. For the three cases, concerning the damping prescription, we consider only the special case where the two dashpots in each system have the same prescribed damping value i.e. η = C 1 = C 2 . We analyze three levels of damping, namely, no damping (η = 0), a low damping (η = 0.218), and a high damping (η = 0.436). Figure 4 depicts a comparison of the wavenumberdependent normalized damped frequency dispersion curve for the three monoatomic configurations, at the three selected damping levels. Figure 4(a) illustrates a near overlap of the frequencies, across the entire BZ, as η varies for the PnC. In figure 4(b), a modest drop in the frequencies is noticed for the piezo PnC without inductor across the entire BZ, as η increases. In figure 4(c), a more noticeable drop in the frequencies takes place for the piezo PnC with inductor across the entire BZ, particularly in the interval κ = [0, 1]. as η increases. Since the mechanical parameters are the same, the three systems are statically equivalent, i.e. the slope of the dispersion curve (group velocity, C g ) in the long-wave limit is the same for an isolated chain. Since in our examples. Since the chains are grounded, i.e. attached to a frame as depicted in figures 2(a) and 3(a), the slope of the curve near κ = 0 is zero for all three models. Figure 5 illustrates a comparison of the damping ratio ζ between the monoatomic PnC, piezo PnC without inductor, and piezo PnC with inductor, for the three prescribed damping levels. When η = 0, ζ = 0 for the PnC for all wavenumbers, but figure 5(a) shows non-zero dissipation in the other systems. This directly demonstrates that adding piezoelectric elements introduces dissipation by a system, which is in fact consistent with the very nature of the energy-harvesting aspect since dissipated energy provides the harvesting source. For the three prescribed damping levels, a significantly higher damping ratio is exhibited by the piezo PnC without inductor and more so by the piezo PnC with inductor, when compared to a non-piezo PnC. As η increases, the damping ratio of all three models increases, as expected. For the piezo PnC with inductor, a discernible change can be seen in the interval κ = [0, 0.5]. For the three prescribed damping levels, the area under the curves have been shaded with different intensities in order to distinguish between the various regions. The bottom shaded region for the damped (η ̸ = 0) cases is directly attributed to the 'raw system damping' since this region is not associated with the presence of piezo elements. The darker colored regions above are therefore, in turn, attributed to dissipation made 'available' for energy harvesting by the addition of the piezoelectric elements (without or with an inductor). These dissipation regions associated with the piezoelectric PnCs are intrinsic in nature, i.e. dependent only on the unit cell properties including the electromechanical coupling and the electrical parameters of the piezoelectric element and not on the dimensions and boundary conditions of a finite structure or any external forcing.
To provide a compact numerical quantification of the useful dissipation added to the PnC due to the utilization of piezoelectric elements, i.e. the level of energy-harvesting availability, we introduce a wavenumber-dependent metric which we refer to as the energy harvest (EH) availability metric, given by: where the subscript ' * ' indicates monoatomic piezo PnC (without inductor or with inductor). We also use its cumulative value over the BZ: and its total value Z tot = Z cum (π) to quantify the overall difference in the damping ratio ζ and hence arrive at a single number to represent the total energy-harvesting availability for a given model.    Figure 6 displays the emerging energy harvest availability metric Z(κ)| * (subscript ' * ' indicating piezo PnC without inductor or piezo PnC with inductor) and its BZ cumulative value, for the three damping levels. The value of Z tot (a single total representative quantity) is also provided. Higher values of Z tot are observed for the piezo PnC with inductor compared to without inductor at all damping levels; this quantifies the inductor's role in improving the energy-harvesting availability. A small increase in Z tot | * is noticed as the damping level increases. This, in turn, indicates that the prescribed damping level in a piezo phononic crystal affects the energy-harvesting availability.

Parameter
Value Unit

Diatomic configurations
In this section, we undergo similar calculations for the diatomic model, also considering the three cases of PnC, piezo PnC without inductor, and piezo PnC with inductor. The parameters used for the diatomic models are given in table 2. Similar to the monoatomic models investigation, we consider only the special case where the four dashpots in each system have the same prescribed damping value, i.e. η = C 1 = C 2 = C 3 = C 4 . We analyze the same three levels of damping as for the monoatomic case. Figure 7 depicts a comparison of the damped dispersion curves for the three diatomic configurations, at the three prescribed damping levels. Similar to the monoatomic case, it is seen in figure 7(a) that for the PnC, the frequencies associated with the acoustic and the optical branches nearly overlap at the three damping levels. In figure 7(b) for the piezo PnC without inductor, a noticeable uniform drop in both branches is observed as η increases. For the piezo PnC with inductor, figure 7(c) depicts even stronger drops especially for the acoustic branch at lower wavenumbers. Since the mechanical parameters are the same, the three systems are statically equivalent, which again provides an appropriate comparison between the cases. Similar to the monoatomic models, the acoustic branches in the long-wave limit are non-zero due to the chains grounding as depicted in figures 2(b) and 3(b). Figure 8 provides a comparison of the damping ratios ζ 1 (acoustic branch) and ζ 2 (optical branch) and their summation ζ sum , between the diatomic PnC, piezo PnC without inductor, and piezo PnC with inductor, at the three prescribed damping levels. When η = 0, ζ 1 , ζ 2 , and ζ sum are equal to zero for the PnC since it is not damped. In contrast, the piezo PnC models exhibit dissipation even in the absence of prescribed damping. As we explained for the monoatomic models, this is direct indication of the energy-harvesting character of these systems as this dissipation stems exclusively from the electromechanical coupling and electrical parameters. For the three configurations, a significant increase in ζ 1 , ζ 2 , and ζ sum is seen as the damping level, η, increases. In analogy to the monoatomic results, the bottom shaded region for the damped (η ̸ = 0) cases are directly attributed to the 'raw system damping' since this region is not associated with the presence of piezo elements. The darker colored regions above, on the other hand, are representative of the intrinsic energy-harvesting availability.
For the diatomic models, the wavenumber-dependent energy harvest availability metric is slightly modified from the monoatomic models to account for the presence of two branches; it is therefore written as: An alternative representation (not adopted here) is to replace the "sum" by an average, i.e., Z avg (κ)| * = ( 2 l=1 Z l (κ)| * )/2. The associated cumulative value, in turn, is given by: and the total value is Z tot l = Z cum l (π). Figures 9 depicts the wavenumber-dependent energy harvest availability metric, its cumulative values, and its total value based on the acoustic branch damping ratio ζ 1 , the optical branch damping ratio ζ 2 , and the summation of the quantities ζ sum , respectively. Similar to the monoatomic models, higher values of Z tot l (l = 1, 2, or sum) are noticed for the piezo PnC with inductor compared to the piezo PnC without inductor at the three damping levels. A relatively small increase is seen in Z tot l | * as the damping level increases, again indicating that the level of prescribed damping has a relatively modest, but noticeable, effect on piezoelectric performance. Upon comparing between the branches, it is observed that ζ 1 for the acoustic branch increases with wavenumber, similar to the monoatomic behaviour. However, in contrast, ζ 2 for the optical branch increases slightly with wavenumber. This bebaviour has implications for frequency-dependent energy harvesting since the optical branch corresponds to higher frequencies.

Parametric analysis of the energy-harvest availability metric
The energy harvest availability metric introduced in this work represents a novel intrinsic characterization tool. In this section, we examine the effects of the various piezoelectric design parameters on this metric. We analyze the behavior of the total energy harvest availability metric Z tot for varying values of the dimensionless parameters used to formulate the governing equations for each piezo PnC model. This enables us to explore the energy-harvesting design space and gain insights on optimal performance.

Monoatomic piezoelectric phononic crystals
In this section, we present a parametric analysis for our monoatomic piezo PnC models at the three selected prescribed damping levels. Figure 10 shows the variation of Z tot with changes in the prescribed damping level η while keeping all the other parameters given in table 1 fixed. As observed earlier, it is clear that increasing the level of prescribed damping improves the energy-harvesting availability, and the rate of this improvement with η is higher when the piezo PnC incorporates an     inductor. We ensure that the prescribed damping levels are chosen such that the piezoelectric PnCs do not become overdamped and remain in the regime of oscillatory vibrations. An η value that renders the piezoelectric PnCs overdamped is considered a cut-off prescribed damping threshold. We find that this threshold is η = 2.92 for the piezo PnC without inductor and η = 1.07 for the piezo PnC with inductor. This points to the ability of the no-inductor piezo PnC chain to endure more prescribed damping than the inductor piezo PnC chain while still sustaining a vibratory response character. However, as noted earlier, the presence of an inductor significantly increases the energy-harvesting availability. Figure 11 shows the variation of Z tot as a function of the dimensionless time constant α for the monoatomic piezoelectric PnCs. Since this quantity is defined as α = ω n C p R l , its value is varied by altering the value of the resistance R l while keeping the rest of the parameters in table 1 constant. Table 3 presents the optimal values of α for the three selected prescribed damping levels. It is noticed that for all three levels, the optimal α value for the piezo PnC with inductor is lower compared to the piezo PnC without inductor; and as η increases, the optimal α value decreases slightly for the two cases. As for the energy-harvesting availability, its optimal value appears to decrease but fairly modestly with increase in η. From α ≈ 0.6 onwards, Z tot | Piezo PnC without inductor begins to increase as η increases, and appears to be on a trajectory to surpass Z tot | Piezo PnC with inductor at higher prescribed damping levels. Figure 12 shows the variation of Z tot as the dimensionless inductor constant β varies for the monoatomic piezoelectric PnCs. Since β = ω 2 n C p L, the value of β is varied by changing the value of the inductance L while keeping the rest of the parameters in table 1 constant. It is seen that for the chosen parameter set, the optimal value of β is 0.3111 for the three selected prescribed damping levels. As η increases, there is no marked change in Z tot | Piezo PnC with inductor . The relationship between L and η is important due to the fact that coupling any inductance with a high damping level has a tendency to push the system towards critical damping and eventually make the system overdamped. The figures broach that for the case of monoatomic piezo PnC with inductor; favorable energy-harvesting availability can be achieved at lower β values and prescribed damping levels. Figure 13 shows the variation of Z tot as the electromechanical coupling θ varies for the monoatomic piezoelectric PnCs. We alter the value of θ, while keeping constant the rest of the parameters in table 1. Table 4 shows the optimal values of θ for the three selected prescribed damping levels. For both the piezoelectric PnCs, it is observed that as η increases, optimal Z tot decreases slightly and optimal θ decreases but insignificantly. For all the three prescribed damping levels, optimal θ is lower for the piezo PnC with inductor than the piezo PnC without inductor. Beyond the optimal value of θ for the piezo PnC without inductor, Z tot is higher for the piezo PnC without inductor, indicating higher energy-harvesting availability compared to the piezo PnC with inductor. From these results, we can deduce that it is favourable to use a piezo PnC with inductor at lower prescribed damping levels although the difference in performance is relatively small.

Diatomic piezoelectric phononic crystals
In this section, we analyze the effect of the previously considered dimensionless parameters on the behavior of Z tot sum for the diatomic piezo PnCs.    Figure 14 shows the variation of Z tot sum , as η changes, while the other parameters given in table 2 remain the same. We find that the cut-off value for the damping level is 2.63 for the piezo PnC without inductor and 0.98 for the piezo PnC with inductor, but it is again shown that higher dissipation is achieved by using a shunt circuit with an inductor at a relatively lower damping level (half of that required for the piezo PnC without inductor) damping level. Figure 15 shows the variation of Z tot sum as α varies for the diatomic piezoelectric PnCs. We utilize the same procedure as for the monoatomic case to obtain the optimal values. Table 5 displays the optimal values of α at the three damping levels. We observe a trend similar to the monoatomic case. For the three damping levels, the optimal α for the piezo PnC with inductor is lower compared to the piezo PnC without inductor; and as η increases, the optimal α decreases and there is no noticeable change in Z tot sum | Piezo PnC with inductor . However, from α ≈ 0.6 onwards, Z tot sum | Piezo PnC without inductor begins to increase and appears to be on a trajectory to surpass Z tot sum | Piezo PnC with inductor at higher η values. Figure 16 shows the variation of Z tot sum as β varies for the diatomic piezoelectric PnCs. We again follow the same procedure as for the monoatomic case to obtain the optimal values. Table 6 shows the optimal values of β at the three prescribed damping levels. The optimal β is the same for η = 0 and η = 0.218, but increases for η = 0.436, which shows that as the prescribed raw damping level increases considerably,   the optimal value of β will increase as well. As η increases, there is no marked change in Z tot sum | Piezo PnC with inductor . The figures indicate that favorable energy-harvesting availability can be achieved at lower β values and prescribed damping levels.
Finally, figure 17 presents the variation of Z tot sum as a function of θ varies for the diatomic piezoelectric PnCs. Table 7 lists the optimal values of θ at the three prescribed damping levels. For both the piezoelectric PnCs, it is observed that as η increases, there is no considerable change in Z tot sum , and the optimal value of θ decreases although insignificantly. At all the three prescribed damping levels, θ is lower for the piezo PnC with inductor compared to the piezo PnC without inductor. Beyond the optimal value of θ for the piezo PnC without inductor, Z tot sum is higher for the piezo PnC without inductor, indicating higher energy-harvesting availability compared to the piezo PnC with inductor. From the plots, we can infer that it is favourable to use a piezo PnC with inductor at lower prescribed damping levels although with very modest comparative gains.

Conclusions
Phononic crystals with piezoelectric patches are capable of harvesting energy from natural and artificial vibrations. This paper proposes an original approach for quantifying the intrinsic energy harvesting availability-a new measure-for a given phononic crystal unit-cell configuration and piezoelectric patch characteristics. This is done by analytically evaluating the difference between the wavenumber-dependent damping ratio for a piezo PnC and a corresponding staticallyequivalent non-piezo PnC. This represents a rigorous BZ characterization approach. The intrinsic energy-harvesting availability is a quantity representative of the intrinsic energyharvesting capacity. Unlike the prevailing approach that examines the extrinsic piezoelectric energy-harvesting capacity of forced finite structures, the proposed approach is based on a material's perspective-no forcing is considered, nor are the overall structure size and global boundary conditions. Four different cases have been investigated by way of demonstration of our new measure, namely, monoatomic piezoelectric PnC without inductor, monoatomic piezoelectric PnC with inductor, diatomic piezoelectric PnC without inductor, and diatomic piezoelectric PnC with inductor. All models are suspended, but the analysis is not restricted to suspended systems. The characteristic equation corresponding to each of these cases has been obtained in closed form. It was shown that the addition of piezoelectric elements (piezo patches) to the monoatomic and diatomic phononic crystals admitting free wave motion introduces 'useful' dissipation, which is essentially the measure mentioned above for the amount of energy available for harvest (i.e. the energy that otherwise would be dissipated into the environment). The amount of raw dissipation has also been quantified (i.e. the energy that gets dissipated into the environment in any case). The specific technical contributions are summarized as follows: • Using Bloch's theorem, a formal characterization of the intrinsic energy-harvesting availability-expressed using the Z(κ) metric-of piezoelectric PnCs under free wave motion has been presented. • Cut-off prescribed damping levels for the piezoelectric PnCs to maintain vibration oscillations have been quantified for a given set of parameters. • Introduction of an inductor in the PnC piezoelectric circuit was demonstrated to substantially increase the energyharvesting availability. • For maximum energy-harvesting availability, optimal values of key non-dimensional electromechanical and electrical design parameters widely employed in the piezoelectric energy-harvesting literature have been computed. This parameter search was driven by the new measure of energyharvesting availability. • For the parameters chosen, the diatomic piezo PnC exhibits generally higher intrinsic energy harvesting availability than the monoatomic piezo PnC when comparing the first damping ratio branch of the former with the sole damping ratio branch of the latter, particularly at high wavenumbers.
The proposed theory provides a fundamental, intrinsic characterization tool for the design and/or selection of piezoelectric energy-harvesting materials from which structure/devices may be formed. The underlying approach for quantifying intrinsic energy-harvesting availability will be useful in formally and generally categorising and designing future energy-harvesting phononic materials, including other types of phononic crystals and also elastic metamaterials. Different configurations from among the vast body of research investigating phononic crystals and elastic metamaterials [25][26][27][28] may be investigated for intrinsic energyharvesting availability. While the concept configurations presented appear large-scale and bulky, and only lumpedparameter models were investigated, future research will easily Adapt this new characterization method to (1) small-scale miniaturized piezoelectric material systems and (2) more complex models.

Data availability statement
All data is producible from the mathematical equations and formulae provided.