Arbitrary-Order Sensitivity Analysis in Phononic Metamaterials Using the Multicomplex Taylor Series Expansion Method Coupled With Bloch ’ s Theorem

Phononic metamaterials (PMs) exhibit frequency ranges at which elastic waves are atten-uated called band gaps. However, this phenomenon is highly sensitive to geometrical var- iations and the unit cell ’ s mechanical properties. It is useful to have accurate sensitivity information to identify the variables that produce the highest impact on band gaps and guide the design of PMs with a desired wave propagation behavior. Current methodologies for sensitivity analysis in PMs, such as the ﬁ nite difference method (FDM), are computationally inef ﬁ cient, subjected to subtraction cancelation errors, and their accuracy is highly dependent on the magnitude of the perturbation step size. In this study, we introduce a new computational methodology to perform parameter sensitivity in the dynamic behavior of PMs using the multicomplex Taylor series expansion (ZTSE) coupled with Bloch ’ s theorem. The methodology allows one to obtain arbitrary-order sensitivities with high accu- racy. In contrast to FDM, this methodology is computationally more ef ﬁ cient, eliminates the step size selection issue, and is not subjected to subtractive cancelation errors. Also, we show how the method can be applied using real algebra solvers. We limit our analysis to linear undamped PMs. The methodology using ZTSE with Bloch ’ s theorem is presented in numerical examples for the diatomic lattice and a 2D square lattice, where we compute up to third-order sensitivities. The results show a maximum normalized root-mean-squared deviation in the order of 10 − 9 for the diatomic lattice and in the order of 10 − 8 for the 2D square lattice when compared to the analytical solutions. [DOI: 10.1115/1.4052830]


Introduction
Phononic metamaterials (PMs) are a class of manmade periodic materials that exhibit unusual dynamic behavior. The mechanical properties of PMs are determined by the geometry, topology, and materials of its repeating unit cell, yielding to Bragg reflections of acoustic or elastic waves. Therefore, these materials exhibit ranges of frequencies at which waves are either allowed to propagate (pass bands) or blocked in one or more directions (stop bands). These particular dynamic behaviors observed in PMs have drawn the attention of researchers in recent years, allowing advances in applications such as waveguides and filters [1], imaging [2], focusing [3], sound collimation [4], energy harvesting [5], negative refraction index [6], negative elastic modulus [7,8], locally resonant sonic materials [9], negative bulk modulus [10], negative mass density [11], peculiar transmission peaks [12], structures for acoustic collimation [13], slow-wave effect [14], and room temperature THz modulators for ultrafast wireless communication [15].
The wave propagation characteristics of PMs are described by their dispersion relations. These relations can be used to identify the wave group velocity, the directionality of the propagating waves, and the bandgaps presented in the material [16,17]. It is widely reported that the dispersion relations for PMs are highly sensitive to mechanical and geometrical variations [18,19]. Therefore, performing sensitivity analyses for these relations is essential for understanding the wave propagation behavior of periodic materials and design goals. For instance, the calculation of sensitivities is a crucial step for designing PMs through topology optimization strategies. However, optimizing the dispersion diagrams of phononic metamaterials is considered a complicated task as it involves calculating sensitivities of multiple eigenvalues [20]. Issues related to the calculation of eigenvalue sensitivities with respect to system design variables have been reported as one of the core limiting factors in the design of PMs using topology optimization [21]. Inaccurate sensitivities result in slow convergence ratios in optimization algorithms or can even halt the optimization process [22]. Different numerical methods have been developed to compute sensitivities in eigenvalue problems. Most of the seminal works in this area have been dedicated to the development of semianalytical methods for the calculation of first-order and second-order sensitivities for the eigenvalues of structural systems [23][24][25][26][27][28][29]. In these works, the eigenvectors are calculated using numerical methods and then used with the sensitivities of the stiffness and mass matrices in algebraic equations to obtain the eigenvalues' sensitivities. It is important to highlight that these methods do not specify a procedure to obtain the sensitivities of the stiffness and mass matrices, and there are no expressions available for calculating higher-order sensitivities of the eigenvalues. Alternatively, one can use the finite difference method (FDM) for the calculation of sensitivities in eigenvalue problems [30]. In FDM, sensitivities are obtained by solving at least two eigenvalue problems: one at the nominal value of the input parameters and the other after applying a real-valued perturbation to the nominal value of the input parameters(s) of interest. The sensitivity of the eigenvalues with respect to the perturbed parameter(s) is then calculated by subtracting these two solutions and dividing by the magnitude of the perturbation used. The implementation of FDM is straightforward; however, the application of this method may lead to subtraction cancelation errors, and its accuracy is highly dependent on the magnitude of the applied perturbation [31][32][33][34]. Despite the existence of methods to compute sensitivities of the eigenvalues (i.e., the modal method), these methods require one to obtain the eigenvectors in a first instance and then use them in conjunction with the sensitivities of the mass and stiffness matrices, which traditionally are obtained with FDM [35]. However, this means that current methods to obtain sensitivities on eigenvalue problems suffer from the same disadvantages inherited from FDM.
The sensitivities of eigenvalue problems can be also calculated using the complex Taylor series expansion (CTSE) method [36,37], a first-order numerical differentiation technique similar in spirit and concept to finite differencing but with significant numerical advantages. CTSE uses the orthogonality of the real and imaginary axes of the complex plane to calculate sensitivities without difference operations; hence, the specification of the step size is not an issue as long as it is sufficiently small, e.g., less that 10 −10 . Unlike finite differencing, CTSE does not require the difference of two analyses. Instead, the sensitivities are computed using a small perturbation along the imaginary axis of the input parameters, solving the eigenvalue problem, and then extracting the imaginary part of the eigenvalues and dividing by the magnitude of the perturbation. In this method, the sensitivity is accurate because no differencing operations are needed. Therefore, the step size can be made arbitrarily small with no concern about round-off error and also that the computation of first-order modal sensitivities is not hampered by repeated eigenvalues [38]. However, CTSE has limited application in the analysis of PMs, as the dispersions relations for these materials are obtained via Bloch's theorem [39,40], which requires the solution of a complex value eigenvalue problem. Therefore, the imaginary axis required for the perturbation is already occupied with the numerical considerations resulting from the application of Bloch's theorem.
In this study, we extend the CTSE method to calculate sensitivities in eigenvalue problems with hypercomplex mathematics. The extension of CTSE to the hypercomplex mathematics method is known as the multicomplex Taylor series expansion (ZTSE) and allows one to compute high-order and mixed sensitivities up to any order [41]. In essence, ZTSE consists of multiple independent imaginary axes that are independent and orthogonal, e.g., i 1 , i 2 , and i 3 . Then, the calculation of high-order and mixed partial sensitivities is obtained by perturbing the parameter of interest along the various imaginary axes. Besides computing higher-order sensitivities, one of the advantages of using hypercomplex mathematics is that one of the complex variable axes can be reserved to apply Bloch's theorem. Therefore, this approach provides a new methodology to enable the sensitivity analysis of the dynamic behavior of linear undamped PMs by integrating multicomplex algebra with Bloch's theorem. Also, as most open-source eigenvalue solvers cannot handle complex algebra operations, we demonstrate this new methodology using real algebra by taking advantage of the Cauchy-Riemann (CR) matrix representation of a complex number [42]. Finally, it is worth highlighting that the ZTSE has been verified in a number of application areas such as elastic analysis [43], fatigue [38,44], thermoelasticity [45], fracture mechanics [46][47][48], plasticity [49], residual stresses [50], creep [51], heat transfer [52], variance estimates [53], and structural dynamics [32], among others; however, to the authors' best knowledge, this is the first demonstration of ZTSE for phenomena described by physics involving complex value mathematics.
The outline for the rest of the article is as follows: Sec. 2 provides a summary of the background information, including the analysis of wave propagation in PMs using Bloch's theorem, the calculation of first-order sensitivities using CTSE, the extension of CTSE to multicomplex mathematics, and the Cauchy-Riemann matrix representation of multicomplex numbers that eliminates the need of using complex variable solvers to obtain sensitivities using ZTSE. Section 3 introduces the methodology to calculate high-order sensitivities in PMs by integrating ZTSE with Bloch's theorem and includes application examples in 1D and 2D mass-spring lattices. Finally, in Sec. 4, we provide the conclusions and future work.

Background
2.1 Wave Propagation in Phononic Metamaterials. The wave propagation in PMs is traditionally studied by applying Bloch's theorem [39,40]. According to this theorem, the dynamic behavior of a PM can be described by analyzing the wave motion of its fundamental unit cell. For a wave propagating without attenuation through such material with spatial periodicity, Bloch's theorem states that the local (i.e., over a unit cell) change in wave amplitude does not depend on the specific location of the cell in the periodic ensemble. Then, the displacement field over a unit cell satisfies [54]: where κ is the wave vector, a is the lattice translation vector, and u(x) and u(x + a) describe the response at the field points x and x + a. The dispersion relation, ω(κ), of a material under free wave motion can be obtained by discretizing its equation of motion in a finite element framework. Considering a linear undamped system, the discretized equation of motion has the form: where ω is the frequency of harmonic wave propagation, K and M are the stiffness and mass matrices, respectively, and U is the nodal displacement vector for the unit cell. Applying the Bloch-Floquet Theorem to this relation results in the generalized eigenvalue problem described by where K R , M R , and u R are the reduced stiffness, mass matrices, and nodal displacement vector, respectively, resulting after removing the redundant equations from Eq. (2). Each assignment of κ yields a particular instance of the ω(κ). The relations ω(κ) provide essential information about the wave propagation behavior in a material, such as wave group velocity, wave directionality, and frequency band gaps. Therefore, obtaining sensitivities for ω(κ) is essential for understanding the behavior of periodic materials and for design goals.
2.2 First-Order Sensitivity Using CTSE. The CTSE is a first-order differentiation technique that uses the orthogonality of the real and imaginary axis of the complex plane to calculate sensitivities. Unlike FDM, CTSE does not require the difference between the two analyses. Instead, the sensitivities are computed using a small (i.e., h ≤ 10 −10 ) perturbation along the imaginary axis [36]. That is, variable X = x 0 is perturbed to X = x 0 + hi 1 , where i denotes an imaginary number and h denotes the perturbation step size. For reference, schematics of the perturbation schemes for FDM and CTSE are shown in Figs. 1(a) and 1(b), respectively. By using the complex Taylor series expansion, a complex function f(x + hi 1 ) can be expanded as follows: where d n f (x) dx n is the nth sensitivity of f(x) with respect to x. Then, taking the imaginary part of both sides of Eq. (4) and ignoring the higher-order terms yields an estimate of the first sensitivity: As there is no subtraction involved, the round-off error is eliminated. On the other hand, the truncation error is reduced by reducing h, e.g., using h ≤ 10 −10 times, the smallest problem parameter will produce a truncation error below the machine precision.

Extension of CTSE to Multicomplex Mathematics:
High-Order Sensitivities. Extending CTSE to multicomplex mathematics allows computing higher-order and mixed sensitivities up to any order [41]. Multicomplex numbers expand the concept of complex numbers by adding multiple imaginary axes to the real quantity (e.g., i 1 , i 2 , i 3 , …, i n ). These imaginary axes are independent and orthogonal [55]. For more details about the properties of multicomplex mathematics, the reader is directed to Ref. [42]. The calculation of high-order and mixed partial sensitivities using the ZTSE is achieved by perturbing the parameter of interest along the imaginary axes, and the number of imaginary directions perturbed determines the maximum order of the sensitivity to be computed. For instance, a second-order sensitivity of f (x) will require a bicomplex number (e.g., two independent imaginary axes) with a small perturbation along the two imaginary axes as shown in Fig. 1(c). By performing the Taylor series expansion of the bicomplex function f (x + h(i 1 + i 2 )), we have The error of the approximation O() depends on the number of terms preserved from the complex Taylor series expansion before the truncation. Considering that multicomplex multiplication is commutative [56], we have that Therefore, Eq. (7) can be simplified to Second-order sensitivities can be obtained by taking the imaginary part corresponding to i 1 i 2 and ignoring the higher-order terms: Similarly, first-order sensitivities are also obtained by taking the imaginary part corresponding to i 1 and ignoring the higher-order terms or by taking the imaginary part corresponding to i 2 and ignoring the higher-order terms: In the case of cross sensitivities, each independent variable is perturbed along different imaginary axes, as shown in Fig. 1(d). Then, the sensitivities of the function f (x, y) are given by (assuming x is perturbed along i 1 and y along i 2 ), Extending this procedure to obtain nth-order sensitivities of a function f (x): It is important to note that the numerical considerations of Bloch's theorem in PMs are already expressed in complex variables (see Eq. (1)). Therefore, it is necessary to reserve one imaginary axis for these periodic boundary conditions and additional imaginary axes for the sensitivities. Then, to compute first-order sensitivities, it is necessary to represent all the input variables at least as bicomplex numbers (e.g., two imaginary axes). For second-order sensitivities, all the variables should be expressed at least as tricomplex numbers (e.g., three imaginary axes). Generalizing, if n imaginary axes are available, it is possible to calculate the sensitivities from order 1 to n−1, e.g., tricomplex numbers allow calculating firstand second-order sensitivities in PMs.

Cauchy-Riemann Matrix Notation of Multicomplex
Numbers. As most numerical packages do not support multicomplex variable operations, it is possible to take advantage of the multicomplex algebra isomorphism to represent any multicomplex number as a matrix with all real numbers. This representation of multicomplex numbers is known as CR matrix notation and allows the use of matrix functions and arithmetic operations as homologous of the operations between multicomplex numbers [42]. For reference, for a multicomplex number with n imaginary axes z η = z η−1 are recursively CR representations of multicomplex numbers with η − 1 imaginary axes: For example, in the case when η = 1, the CR representation of the complex number z η=1 = x 1 + x 2 i 1 , where z η=1 ∈ ℂ 1 and x 1 , x 2 ∈ ℝ are shown in Eq. (14).
In the case of a bicomplex number with η = 2 imaginary axes, The CR representation is given by Eq. (15): For a tricomplex number with η = 3 imaginary axes, The CR is given as follows: 3 Sensitivity Analysis in Phononic Materials Using ZTSE We propose a new methodology that combines ZTSE with Bloch's theorem described in Eqs. (1) and (3) to obtain sensitivity information on the dispersion diagrams that describe the dynamic behavior of PMs. The methodology is summarized in Fig. 2(a), and it requires converting all the parameters involved in the problem (e.g., mass, stiffness, geometry, and nodal coordinates) to multicomplex variables represented as real-value matrices using the CR notation. This conversion facilitates the operations between multicomplex quantities by allowing the use of matrix operations, eliminating the necessity of supporting multicomplex algebra by the numerical package used. When converting all the parameters to the CR notation, it is important to recall that one imaginary axis should be reserved to apply periodic boundary conditions using Bloch's theorem. Therefore, n+1 imaginary axes should be considered to obtain sensitivities of order n. Then, to get sensitivities, imaginary perturbations are added to the parameters of interest. These perturbations must be incorporated using their CR notation. To minimize truncation errors, the perturbation magnitudes must be a low number (e.g., h ≤ 10 −10 ) times the nominal value of the variable of interest.
To exemplify this procedure, consider the discrete system with two masses shown in Fig. 2(b), representing the unit cell of an infinitely periodic monoatomic lattice. Each point mass is represented with a different shape: circles for m 1 , squares for m 2 , and the stiffness of the springs, c, is represented with parallelograms. Considering first-order sensitivities for any input parameter, it is necessary to use bicomplex numbers. Each axis representing the bicomplex numbers is identified with different colors: orange for the real axis, blue for i 1 , red for i 2 , and green for the cross-terms i 12 .
In this example, we are reserving i 1 for sensitivities with respect to any input parameter denoted from now on with i s , and i 2 is reserved for applying Bloch's theorem and denoted with i B . In Fig. 2(c), we present the assembled generalized eigenvalue problem corresponding to Eq. (2). As shown, each position of the mass matrix is replaced with a CR representation for each bicomplex mass. The displacement vectors contain all the real and imaginary degrees-of-freedom, and the mass and stiffness matrices are constructed using the finite element formulation of spring elements as follows: After applying boundary conditions using Bloch's theorem, the dispersion diagram for the diatomic system can be obtained by solving the generalized eigenvalue problem described in Fig. 3, where a is the size of the unit cell and κ is the wavenumber.
The resulting system after applying Bloch's theorem is nonself-adjoint. Therefore, the QR algorithm is used to solve the eigenvalue problem [57][58][59]. The QR algorithm is recognized as one of the most important algorithms for eigenvalue computations due to its robustness [60][61][62]. The most traditional numerical algorithms to solve eigenvalue problems achieve the stability for selfadjoint, Hermitian, or tridiagonal Hermitian eigenvalue problems [63]. This includes the Davidson method [64], the implicit restarted Lanczos [65], and the locally optimal block preconditioned conjugate gradient [66]. In this study, we aim to solve a nonself-adjoint eigenvalue problem. Nonself-adjoint eigenvalue problems are more difficult and computationally demanding to solve than selfadjoint problems [67]. Generally, self-adjoint eigenvalue problems are well conditioned, while nonself-adjoint problems are not [68]. Schur decomposition-based methods to solve eigenvalue problems, such as the QR algorithm, have been standard methods for solving small-to medium-sized eigenvalue problems [69,70]. The reader is referred to Ref. [63] for a thorough comparison of the available Open-source libraries to solve nonself-adjoint eigenvalue problems of large-scale power systems.
The QR algorithm requires converting the generalized eigenvalue problem into a standard eigenvalue problem. To this end, both sides of Eq. (3) are premultiplied by the negative of the inverse of the reduced mass matrix to obtain The eigenvalues of Eq. (18) If the size of the matrix [S] 0 is N 2 , the QR algorithm requires N 3 iterations to find the eigenvalues. In a traditional analysis, the QR algorithm would produce an upper triangular matrix with the eigenvalues located along the diagonal of the resulting matrix [S] N 3 . However, in the present analysis, the multicomplex eigenvalues along the diagonal are replaced with CR matrices. The number of CR matrices inside [S] N 3 is equal to the number of degrees-of-freedom in the reduced system (Eq. (19)). To illustrate this point, Eq. (20) presents a system with N degrees-of-freedom, and each CR matrices along the diagonal of [S] N 3 represents bicomplex numbers. In this system, the eigenvalues and the corresponding imaginary terms with the sensitivities information are in the first column of each CR matrix. For this example, Im 1 is reserved for the first-order sensitivity with respect to any input parameter and denoted as Im S , and Im 2 is reserved for the imaginary axis required by the periodic boundary conditions from Bloch's theorem and denoted as Im B .
Given the matrix [S] N 3 , one can obtain the dispersion diagrams, ω(κ), and the corresponding sensitivities by calculating the square root of the resulting CR matrices of the eigenvalues using the Denman and Beavers iteration method [71]. Other methods used for calculating the square root of matrices [72] are known to fail when CR is used to represent multicomplex numbers [56].
By applying Eqs. (10) and (11) to the color-highlighted terms shown in Eqs. (20) and (21), the eigenvalues and their corresponding first-order sensitivities are obtained as follows: For reference, the implementation to compute second-order sensitivities of this example is shown in the Supplementary Material I on the ASME Digital Collection.
The CR matrix notation of multicomplex numbers allows to expand the complex standard eigenvalue problem presented in Ref. [73] to include general damping in the formulation of the sensitivity analysis. However, the QR algorithm produced inaccurate results when compared with analytical expressions of PMs with damping. Hence, further research is necessary to better understand the impact of the eigenvalue problem solver used in the analysis of PMs with damping. For lucidity, we relegate this procedure for future work and present the implementation of ZTSE for linear undamped discrete systems.

Extension to Continuum Finite
Elements. The use of ZTSE can be extended to continuum finite elements, and this has been called the multicomplex finite element method (ZFEM) [74].
They key aspect of ZFEM is that the nodal coordinates are multicomplex variables, having real and imaginary components. The real components define the physical location of the nodes (similarly as in traditional finite elements method), and the imaginary components define the "perturbations" of the real (physical) nodal coordinates, with the specific perturbation defining the derivative to be computed. To illustrate this concept, Fig. 4 shows an example of a 16-noded quadrilateral element with four real nodes, four nodes corresponding to the first imaginary axis, four nodes corresponding to the second imaginary axis, and four more nodes corresponding to the mixed imaginary axis 1,2. The other input parameters, such as material properties (i.e., density, Young's modulus, Poisson's ratio) and boundary conditions are also multicomplex variables. In the process of computing the mass and stiffness matrices, each multicomplex variable is represented with the CRM notation of multicomplex numbers.
ZFEM can be implemented in the sensitivity analysis of structural systems with large number of degrees-of-freedom, while preserving high accuracy; however, it is important to reiterate that the convergence rate of the QR algorithm is N 3 , where N is the total number of degrees-of-freedom, and as ZTSE multiplies the number of real degrees-of-freedom by 2 η , where the total number of imaginary axis η is given by the highest order of sensitivity to compute plus the imaginary axis required by Bloch's theorem. This means that the computational runtime grows as function of the number of degrees-of-freedom used and the highest order of sensitivity to compute. An additional consideration is that the accuracy of the derivatives is dependent on the accuracy of the real part; hence, increasing the number of nodes will increase the accuracy of the real part and its sensitivities. In this study, we limit our analysis to discrete systems due to their simplicity allowing an easy understanding and demonstration of the proposed methodology. The implementation and verification using other continuum elements will be relegated to future works.

Case Studies
To illustrate this new methodology for the sensitivity analysis of PMs, we apply the method to the sensitivity analysis of the diatomic and a 2D squared lattice in the following sections. The numerical results from the ZTSE and Bloch's theorem methodology are compared with analytical equations, and the errors are measured using the normalized root-mean-square deviation (NRMSD) as described in Eq. (23). Also, a comparison of the runtime between the proposed methodology and FDM is presented. 4.1 Diatomic Lattice. The diatomic lattice is an infinite array of alternating masses m 1 and m 2 interconnected with springs with linear stiffness c. Also, the two consecutive masses of the same kind are separated a distance a from each other, and each mass m 1 and m 2 possesses one translational degree-of-freedom u s and v s , respectively. Figure 5 shows a schematic of the diatomic lattice. The analysis of wave motion for this system can be easily derived by obtaining the equations of motion for two neighboring masses, assuming the absence of external forcing, and assuming harmonic motion. Then, after applying the periodic conditions from Bloch's theorem, the dispersion relation for this system is given by Since the diatomic lattice has two active degrees-offreedom, the dispersion relation has two branches called optical and acoustic modes denoted in Eq. (24) as ω O and ω A , respectively.
Next, we present the detailed numerical implementation of our method as illustrated in the flowchart of Fig. 2(a). First, we perturb the variable of interest, and second, we use CR matrix representation for each input variable. Out of briefness, we only present the procedure for first-order sensitivities of m 2 , but the procedure can be extrapolated for the other variables and higherorder sensitivities as well. In this example, we use the first imaginary axis for the sensitivities and the second imaginary axis for imposing Bloch's periodic boundary conditions denoted with i B . With this in mind, the CR matrix representation for the input parameters is as follows:    The next step is to assemble the global mass and stiffness matrices of the unit cell. For this example, we use the spring finite element formulation. The resulting generalized eigenvalue problem is expressed as follows: where the vectors {u s }, {v s }, and {u s+1 } contain the displacement degrees-of-freedom of the corresponding masses presented in Fig. 5. For the diatomic lattice, the wave propagation is in one dimension only, and therefore, a single wavenumber is enough to describe the irreducible Brillouin zone. For this reason, to apply Bloch-Floquet periodic boundary conditions, we have: The matrix exponential function algorithm used is further detailed in Refs. [75,76] and yields correct results for CR matrix representation of multicomplex numbers. Finally, we solve the generalized eigenvalue problem for the reduced system: Equation (28) is solved following the procedure in Eq. (19). The information of the modes and its sensitivities are found in the system matrix [S] as previously stated in Eqs. (20)- (22).
Journal of Applied Mechanics FEBRUARY 2022, Vol. 89 / 021007-9 4.1.1 Sensitivity Analysis of Diatomic Lattice. We implemented the multicomplex-step perturbation method to compute the dispersion relation in Eq. (24) and its sensitivities using the following input values, m 1 = 1, m 2 = 2, c = 1, and a = 1. With perturbation magnitude, h = 10 −10 times the nominal value of the perturbed variable. The error quantification was performed with Eq. (23). Figure 6 shows the comparison of the dispersion relation and firstand second-order sensitivities obtained using the symbolic algebra package of MATLAB and the numerical results obtained with ZTSE. Excellent agreement was found with both methods, resulting in NRMSD values of 5 × 10 −14 for the dispersion relation, 6 × 10 −12 , and 10 −10 for the firstand second-order sensitivities with respect to m 2 , respectively, and 9 × 10 −11 for the second-order mixed sensitivity with respect to m 2 and a. The remaining comparisons for the first-, second-, and third-order sensitivities can be found in the Supplementary Material II on the ASME Digital Collection. The largest NRMSD value obtain for third-order derivatives was 9 × 10 −9 .

Runtime Comparison for Diatomic Lattice.
A comparison was conducted between ZTSE and FDM comparing the runtime for computing the first-, second, and third-order sensitivity of the dispersion relation with respect to the mass m 2 . Hence, the mean runtime of 1000 runs were used, and the results were normalized with respect to the runtime required by ZTSE to compute the first-order sensitivity. The results are listed in Table 1. We see a lower runtime necessary for ZTSE than FDM to compute the firstand second-order sensitivities. However, for third-order sensitivities, FDM took less computational runtime than ZTSE. These runtimes are equivalent for sensitivities on other input parameters.

Perturbation
Step Convergence for Diatomic Lattice. Figure 7 shows the behavior of the error of the sensitivities as a function of the perturbation step size for FDM and ZTSE. The truncation error decreases as the perturbation step decreases for both methods until a critical point, from which the subtractive error starts dominating the inaccuracy of FDM, increasing the error of FDM as h decreases. In contrast, for ZTSE, the truncation error falls below machine precision, and the sensitivity reaches its maximum precision, from which further decreases in h do not result in accuracy increases for ZTSE.

2D Square Lattice.
A rectangular cell was analyzed that consisted of four masses m connected with horizontal and vertical springs with stiffness c x and c y , respectively. Springs with effective stiffness c d connect the diagonals of the lattice. The rectangular cell has horizontal length a and vertical length b. The rectangular cell shown in Fig. 8 represents the unit cell of the infinitely periodic cell.
The following expression gives the analytical dispersion relation for the 2D square lattice with diagonal springs: Next, we present the detailed numerical implementation of our method, as illustrated in the flowchart of Fig. 2(a). First, we perturb the variable of interest, and second, we use CR matrix representation for each input variable. Out of briefness, we only present the procedure for first-order sensitivities of m, but the procedure can be extrapolated for the other variables and higher-order sensitivities. In this example, we use the first imaginary axis for the sensitivities and the second imaginary axis for imposing Bloch's periodic boundary conditions. With this in mind, the CR matrix representation for the input parameters is expressed as follows: 021007-10 / Vol. 89, FEBRUARY 2022 Transactions of the ASME The next step is to assemble the global mass and stiffness matrices of the unit cell. For this example, we use the spring finite element formulation. The resulting generalized eigenvalue problem is expressed as follows: For the 2D square lattice, it is necessary to use two wavenumbers to describe the first Brillouin zone. For this reason, to apply Bloch-Floquet periodic boundary conditions, we have: Finally, we solve the generalized eigenvalue problem for the reduced system: Equation (36) is solved following the procedure in Eq. (19). For this problem, each degree-of-freedom carries the information of the wave propagating along the x and y directions. Since our interest is on the magnitude of the frequency, it is necessary to add both degrees-of-freedom before using the matrix square root. The information of the modes ω 2 x and ω 2 y and their sensitivities are found in the system matrix [S] as previously stated in Eqs. (20)- (22).
4.2.1 Sensitivity Analysis of the 2D Square Lattice. We implemented the multicomplex-step perturbation method to compute the dispersion relation in Eq. (30) and its sensitivities using the following input values, m = 1, c x = 1, c y = 1, c d = 1, a = 1, and b = 1, with perturbation magnitude h = 1 × 10 −10 . The error quantification was performed with Eq. (23). Figure 9 shows excellent agreement between the numerical solution obtained with ZTSE and MATLAB's symbolic algebra for the dispersion relation and firstand second-order sensitivities. It is important to note here that when κ x = κ y , for the input parameter values selected, ω x = ω y , which means the system produces repeated eigenvalues that were successfully obtained with our method. Excellent agreement was found with both methods, resulting in NRMSD values of 10 −16 for the dispersion relation, 10 −15 for both, the first and second-order sensitivities with respect m, respectively, and 5 × 10 −15 for the second-order mixed sensitivity with respect to m and a. The remaining first-, second-, and third-order sensitivities can be found in the Supplementary Material III on the ASME Digital Collection. The largest NRMSD value obtain for third-order derivatives was 3 × 10 −8 .

Runtime
Comparison for the 2D Square Lattice. A runtime comparison between ZTSE and FDM for the first-, second-, and third-order sensitivities of the dispersion relation with respect to the mass m was conducted. For this purpose, the mean runtime of 1000 runs were used and normalized with respect to the runtime required by ZTSE to compute the first-order sensitivities. The results are listed in Table 2. Overall, a lower runtime is required by ZTSE than FDM to compute the sensitivities of the dispersion relation for first-, second-, and third-order sensitivities. In contrast to the diatomic lattice, in the 2D square lattice, ZTSE took less computational runtime to compute the third-order sensitivity than FDM.

Perturbation
Step Convergence for the 2D Square Lattice. Figure 10 shows the behavior of the sensitivities error as a function of the perturbation step size for FDM and ZTSE. The truncation error decreases as the perturbation step decreases for both methods until a critical point, from which the round-off error starts dominating the inaccuracy of FDM, increasing the error of FDM as h decreases. In contrast, for ZTSE, the truncation error falls below the machine precision number, and the sensitivity reaches its maximum precision, from which further decrease in h does not result in an accuracy increase for ZTSE.

Conclusions
In this study, we presented a new method that integrates the multicomplex Taylor series expansion with Bloch's theorem to enable the arbitrary-order sensitivity analysis of phononic metamaterials. The capability of obtaining high-accuracy higher-order sensitivities has potential applications in improving current optimization strategies in the design of PMs and enabling the design of parameter insensitive PMs. We verified this new methodology by analyzing the sensitivities of two discrete PMs: (i) the diatomic lattice and (ii) a 2D squared lattice with diagonal springs. Both systems' dispersion relations and sensitivities were calculated using our new methodology and compared with the finite differences method (FDM) and analytical functions.
We found that numerical sensitivities obtained with ZTSE and Bloch are significantly more accurate than those obtained with FDM. ZTSE proved to be perturbation step size independent. In contrast, FDM proved to be prone to truncation and round-off errors, leading to grossly inaccurate sensitivities. Moreover, first-, second-, and third-order numerical sensitivities obtained with ZTSE took less computational runtime than those obtained with FDM.
Additional advantages of the proposed method are that it does not require the computation of eigenvectors to obtain the sensitivities of the eigenvalues, nor to know the sensitivities of the stiffness and mass matrices a priori. Also, ZTSE allows for simultaneous computation of eigenvalues and their respective sensitivities in a single procedure, and it is general enough to compute arbitrary-order sensitivities. Furthermore, we showed that the coupled ZTSE-Bloch methodology can take advantage of the Cauchy-Riemann matrix notation for multicomplex numbers, eliminating the need to use eigenvalue problem solvers that handle complex variables. Discrete systems with local resonators can also be studied with ZTSE by increasing the number of degrees-of-freedom.
We posit that the method presented in this study will help close the breach between theoretical predictions and practical developments of PMs by providing a handy yet highly accurate method to perform sensitivity analysis of the dynamical behavior of PMs. Moreover, this methodology will help to understand better the impact of defects and other unintended variations in PMs and will aid the design of new PMs with unique dynamic properties.
Future research efforts will focus on extending this method into available real-value finite element software, including sensitivities to perform uncertainty quantification in PMs, and integrating this strategy to guide the design of new PMs with novel properties. The inclusion of viscous damping promises the ability to compute accurate sensitivities in more realistic analysis. However, further research is necessary to better understand the impact of the eigenvalue problem solver used for the analysis of PMs with damping. Similarly, the presented methodology is limited to the analysis of linear PMs. Future research efforts will commit on the implementation of ZTSE for the sensitivity analysis of nonlinear PMs. a = unit cell length c = Springs' stiffness h = perturbation step magnitude n = order of sensitivities ℝ = set of real numbers I = identity matrix i B = imaginary axis reserved for Bloch's theorem i s = imaginary axis reserved to the computation of sensitivities m 1 = mass corresponding to node 1 u R = reduced displacement vector x i1 = imaginary part of x corresponding to first imaginary axis x Re = real component of x z η = multicomplex number of order η ℂ 1 = set of complex numbers ℂ 2 = set of bicomplex numbers K R = reduced stiffness matrix M R = reduced mass matrix f (x) = arbitrary function evaluated at x i 1 , i 2 , i 3 = imaginary axes 1,