Non-reciprocal wave propagation in discretely modulated spatiotemporal plates

TL;DR

This paper introduces a general plane wave expansion method for analyzing non-reciprocal wave propagation in discretely modulated spatiotemporal plates, demonstrating how such materials can break symmetry and violate reciprocity in wave transmission.

Contribution

It presents a novel, general formulation of the Plane Wave Expansion Method applicable to 2D spatiotemporal unit cells and explores materials that break mirror symmetry in momentum space.

Findings

Validated dispersion profiles with numerical simulations.

Demonstrated non-reciprocal wave propagation in modulated plates.

Showed the ability to violate reciprocity principle in wave transmission.

Abstract

We investigate non-reciprocal wave propagation in spatiotemporal phononic plates. In particular, the first goal of this manuscript is to present a general formulation of the Plane Wave Expansion Method (PWEM) that, in contrast with previous works, is applicable to any class of 2D spatiotemporal unit cells whose properties can be expanded in traveling plane waves. The second goal is to exploit this analysis tool in order to study a new class of materials capable of violating mirror symmetry in momentum space, therefore breaking reciprocity principle along different wave propagation directions. This is obtained by considering the plate elastic properties to be discretely modulated in space and continuously in time. Theoretical dispersion profiles are validated and compared with numerical simulations.

Equations

$$\frac{\partial^{2}m_{x}}{\partial x^{2}}+2\frac{\partial^{2}m_{xy}}{\partial x\partial y}+\frac{\partial^{2}m_{y}}{\partial y^{2}}=-\frac{\partial}{\partial t}\bigg{[}s(x,y)\rho(x,y,t)\frac{\partial w(x,y,t)}{\partial t}\bigg{]}$$

$$\begin{split}m_{x}&=B(x,y,t)\bigg{[}\frac{\partial^{2}w(x,y,t)}{\partial x^{2}}+\nu\frac{\partial^{2}w(x,y,t)}{\partial y^{2}}\bigg{]}\\[4.0pt] m_{y}&=B(x,y,t)\bigg{[}\frac{\partial^{2}w(x,y,t)}{\partial y^{2}}+\nu\frac{\partial^{2}w(x,y,t)}{\partial x^{2}}\bigg{]}\\[4.0pt] m_{xy}&=(1-\nu)B(x,y,t)\frac{\partial^{2}w(x,y,t)}{\partial x\partial y}\end{split}$$

$$\begin{array}{r}B(x,y,t)=\sum _{m,n,v=-\infty }^{\infty }{\hat{B}}_{m,n,v}{\mathrm{e}}^{\mathrm{j}({\mathbf{\kappa }}_{\mathbf{m}}\cdot \mathbf{r}-v{\omega }_{m}t)}\\ G(x,y,t)=\sum _{m,n,v=-\infty }^{\infty }{\hat{G}}_{m,n,v}{\mathrm{e}}^{\mathrm{j}({\mathbf{\kappa }}_{\mathbf{m}}\cdot \mathbf{r}-v{\omega }_{m}t)}\end{array}$$

$$\begin{array}{r}{\hat{B}}_{m,n,v}=\frac{1}{T_m{\lambda }_{mx}{\lambda }_{my}}{\int }_{D}B(x,y,t){\mathrm{e}}^{-\mathrm{j}({\mathbf{\kappa }}_{\mathbf{m}}\cdot \mathbf{r}-v{\omega }_{m}t)}\mathrm{dD}\\ {\hat{G}}_{m,n,v}=\frac{1}{T_m{\lambda }_{mx}{\lambda }_{my}}{\int }_{D}G(x,y,t){\mathrm{e}}^{-\mathrm{j}({\mathbf{\kappa }}_{\mathbf{m}}\cdot \mathbf{r}-v{\omega }_{m}t)}\mathrm{dD}\end{array}$$

$$w(x,y,t)=\hat{w}(x,y,t){\mathrm{e}}^{\mathrm{j}(\mathbf{k}\cdot \mathbf{r}-\omega t)},\hat{w}(x,y,t)=\sum _{p,q,r=-\infty }^{\infty }{\hat{W}}_{p,q,r}{\mathrm{e}}^{\mathrm{j}({\mathbf{\kappa }}_{\mathbf{m}}\cdot \mathbf{r}-r{\omega }_{m}t)}$$

$$\left[{\tilde{\mathbf{L}}}_0(k_x,k_y)+{\tilde{\mathbf{L}}}_1\omega +{\tilde{\mathbf{L}}}_2{\omega }^2\right]\tilde{\mathbf{w}}=0$$

$$\begin{split}&E(x,y,t)=E_{0}\bigg{\{}1+\frac{\alpha_{m}}{2}\left[\cos\bigg{(}(i-1)\frac{2\pi}{R_{s}}-\omega_{m}t\bigg{)}\right]+\frac{\alpha_{m}}{2}\left[\cos\bigg{(}(j-1)\frac{2\pi}{R_{s}}-\omega_{m}t\bigg{)}\right]\bigg{\}}\quad\\[4.0pt] &where:\;\;\;i=1,...,R_{s}\hskip 28.45274ptj=1,...,R_{s}\end{split}$$

$${\kappa }_m=\frac{{\kappa }_{mx}{\kappa }_{my}}{\sqrt{{\kappa }_{mx}^2+{\kappa }_{my}^2}}{\mu }_{\gamma }=\sqrt{{\mu }_x^2+{\mu }_y^2}={\mu }_x\sqrt{1+\tan {\gamma }^2}$$

$${\nu }_{\gamma }=\mathbf{\nu }\cdot \frac{\mathbf{\kappa }}{|\mathbf{\kappa }|}$$

$$\displaystyle B\bigg{[}\frac{\partial^{4}w}{\partial x^{4}}+2\frac{\partial^{4}w}{\partial x^{2}\partial y^{2}}+\frac{\partial^{4}w}{\partial y^{4}}\bigg{]}+2\frac{\partial B}{\partial x}\frac{\partial}{\partial x}\bigg{[}\frac{\partial^{2}w}{\partial x^{2}}+\frac{\partial^{2}w}{\partial y^{2}}\bigg{]}+2\frac{\partial B}{\partial y}\frac{\partial}{\partial y}\bigg{[}\frac{\partial^{2}w}{\partial x^{2}}+\frac{\partial^{2}w}{\partial y^{2}}\bigg{]}+$$

$$\displaystyle+\bigg{[}\frac{\partial^{2}B}{\partial x^{2}}+\frac{\partial^{2}B}{\partial y^{2}}\bigg{]}\bigg{[}\frac{\partial^{2}w}{\partial x^{2}}+\frac{\partial^{2}w}{\partial y^{2}}\bigg{]}-(1-\nu)\bigg{[}\frac{\partial^{2}B}{\partial x^{2}}\frac{\partial^{2}w}{\partial y^{2}}-2\frac{\partial^{2}B}{\partial x\partial y}\frac{\partial^{2}w}{\partial x\partial y}+\frac{\partial^{2}B}{\partial y^{2}}\frac{\partial^{2}w}{\partial x^{2}}\bigg{]}=$$

$$=-\frac{\partial G}{\partial t}\frac{\partial w}{\partial t}-G\frac{{\partial }^{2}w}{\partial t^2}$$

$$\displaystyle\sum\limits_{m,n,v,p,q,r}\hat{B}_{m,n,v}\hat{W}_{p,q,r}\bigg{\{}\bigg{[}(pk_{mx}+k_{x})^{2}+(qk_{my}+k_{y})^{2}\bigg{]}\bigg{[}([m+p]k_{mx}+k_{x})^{2}+([n+q]k_{my}+k_{y})^{2}\bigg{]}+$$

$$\displaystyle-(1-\nu)\bigg{[}(mk_{mx})(qk_{my}+k_{y})-(nk_{my})(pk_{mx}+k_{x})\bigg{]}^{2}\bigg{\}}\cdot{\rm{e}}^{{\rm{j}}([m+p]k_{mx}x+[n+q]k_{my}y-[v+r]\omega_{m}t)}=$$

$$\displaystyle=\sum\limits_{m,n,v,p,q,r}\hat{G}_{m,n,v}\hat{W}_{p,q,r}\bigg{\{}(r\omega_{m}+\omega)([v+r]\omega_{m}+\omega)\bigg{\}}\cdot{\rm{e}}^{{\rm{j}}([m+p]k_{mx}x+[n+q]k_{my}y-[v+r]\omega_{m}t)}$$

$$\displaystyle\sum\limits_{p=-P}^{P}\sum\limits_{q=-Q}^{Q}\sum\limits_{r=-R}^{R}\hat{B}_{a-p,b-q,c-r}\hat{W}_{p,q,r}\bigg{\{}\bigg{[}(pk_{mx}+k_{x})^{2}+(qk_{my}+k_{y})^{2}\bigg{]}\bigg{[}(ak_{mx}+k_{x})^{2}+(bk_{my}+k_{y})^{2}\bigg{]}+$$

$$\displaystyle-(1-\nu)\bigg{[}([a-p]k_{mx})(qk_{my}+k_{y})-([b-q]k_{my})(pk_{mx}+k_{x})\bigg{]}^{2}\bigg{\}}=$$

$$\displaystyle=\sum\limits_{p=-P}^{P}\sum\limits_{q=-Q}^{Q}\sum\limits_{r=-R}^{R}\hat{G}_{a-p,b-q,c-r}\hat{W}_{p,q,r}\bigg{\{}(r\omega_{m}+\omega)(c\omega_{m}+\omega)\bigg{\}}$$

$$\displaystyle\sum\limits_{p=-P}^{P}\sum\limits_{q=-Q}^{Q}\bm{\tilde{K}}_{a-p,b-q}\bigg{\{}\bigg{[}(pk_{mx}+k_{x})^{2}+(qk_{my}+k_{y})^{2}\bigg{]}\bigg{[}(ak_{mx}+k_{x})^{2}+(bk_{my}+k_{y})^{2}\bigg{]}+$$

$$\displaystyle-(1-\nu)\bigg{[}([a-p]k_{mx})(qk_{my}+k_{y})-([b-q]k_{my})(pk_{mx}+k_{x})\bigg{]}^{2}\bigg{\}}\bm{\tilde{w}}_{p,q}=$$

$$\displaystyle=\sum\limits_{p=-P}^{P}\sum\limits_{q=-Q}^{Q}\bigg{\{}\bm{\tilde{M}}^{0}_{a-p,b-q}+\bm{\tilde{M}}^{1}_{a-p,b-q}\omega+\bm{\tilde{M}}^{2}_{a-p,b-q}\omega^{2}\bigg{\}}\bm{\tilde{w}}_{p,q}$$

$${\tilde{\mathbf{K}}}_{a-p,b-q}=$$

$${\tilde{\mathbf{M}}}_{a-p,b-q}^0=$$

$${\tilde{\mathbf{M}}}_{a-p,b-q}^1=$$

$${\tilde{\mathbf{M}}}_{a-p,b-q}^2=$$

$${\tilde{\mathbf{w}}}_{p,q}=\{w_{p,q,-R},\dots ,w_{p,q,+R}{\}}^T$$

$$\sum\limits_{p=-P}^{P}\bigg{(}\bm{\tilde{Z}}^{0}_{a-p}+\bm{\tilde{Z}}^{1}_{a-p}\omega+\bm{\tilde{Z}}^{2}_{a-p}\omega^{2}\bigg{)}\bm{\tilde{w}}_{p}=0$$

$${\tilde{\mathbf{Z}}}_{a-p}^0=$$

$${\tilde{\mathbf{Z}}}_{a-p}^1=$$

$${\tilde{\mathbf{Z}}}_{a-p}^2=$$

$$\left[\bm{\tilde{L}}_{0}(k_{x},k_{y})+\bm{\tilde{L}}_{1}\omega+\bm{\tilde{L}}_{2}\omega^{2}\right]\bm{\tilde{w}}=0\$$

$${\tilde{\mathbf{L}}}_0(k_x,k_y)=$$

$${\tilde{\mathbf{L}}}_1=$$

$${\tilde{\mathbf{L}}}_2=$$

Full text

1 Introduction

Non-reciprocal wave propagation in phononic crystals has drawn growing attention within the research community in the past years [1, 2]. Indeed, the opportunity to design elastic structures which support one-way wave propagation can be particularly relevant for next-generation applications involving elastic energy manipulation, waveguiding and conversion [3, 4].

In this context, 1D phononic waveguides have recently been designed to break the reciprocity principle. One possibility is to locally alter material properties mimicking the propagation of a wave [5, 6, 7, 8, 9], which generally requires active elements. A second possibility is to leverage 2D lattice structures with spinning gyroscopes, establishing elastic analogues of the Quantum Hall Effect (QHE) [10, 11, 12]. Specifically, this strategy relies on the exploitation of one-way topologically protected edge waves, thus not involving the bulk of the material. Systems with broken time reversal symmetry can be also designed using three port devices [13] or exploiting nonlinearities [14, 15].

Non-reciprocal elastic wave manipulation in 2D structures has been accomplished by Attarzadeh et al. [16] leveraging spatiotemporally modulated membranes to break time reversal symmetry, continuously biasing the applied traction field and material density in space and time. Compared with 1D structures, the addition of one spatial dimension makes wave propagation analysis more complex, but also opens new promising ways to manipulate elastic waves along different directions.

Although the aforementioned examples offer a theoretical framework to study nonreciprocity in a specific class of 2D spatiotemporal mechanical structures, additional efforts are required for the analysis of arbitrarily shaped modulations. In pasticular, in this manuscript our aim is to bridge the gap between fully continuous modulations (for which material properties are pointwise varying in time) and those unit cell profiles that can be more easily realized in engineering applications - which must be studied with appropriate tools.

Motivated by this consideration, we investigate discretely modulated plates. This new class of materials is characterized by spatially discrete and temporally continuous elasticity profiles, mimicking the propagation of a plane wave in two dimensions.

In contrast with already studied spatiotemporal materials, the proposed structure allows for a more feasible implementation, which is supposed to be provided using established techniques, such as negative capacitance shunts applied to piezoelectric materials bonded on a 2D substrate [17, 18]. The corresponding non-symmetric band diagram is computed using a generalization of the Plane Wave Expansion Method (PWEM) provided in [16] extending our previous work [7] to 2D structures. The generalized PWEM can be thus applied to any kind of 2D spatiotemporal unit cells whose material properties can be projected in a Fourier basis - and not limited to any specific class of unit cell profiles - therefore assuming that the modulation can be mapped as a series of plane waves characterized by different propagation speed.

We demonstrate that discretely modulated plates support non-reciprocal wave propagation. Specifically, this is accomplished highlighting the leading Bloch-wave component [7, 6, 15] of the dispersion diagram, which is characterized by a wide set of supported modes. Once the diagram is filtered to the most relevant terms, associated directivity plots and group velocities are computed to predict directional wave propagation phenomena. Theoretical dispersion plots are verified through numerical simulation of wave propagation and actual dispersion, which is reconstructed from the displacement response of a plate under tone burst excitation.

The article is organized as follows: in section 2 we describe the analytical procedure for the band diagram computation. In section 3, theoretical and numerical results are compared and discussed. Concluding remarks are presented in section 4.

2 Analytical procedure for non-symmetric band diagram computation

Consider the general equation describing the out of plane dynamics of a Kirchoff plate, therefore neglecting in-plane polarized motion:

[TABLE]

where $w(x,y,t)$ is the out of plane displacement field, $s(x,y)$ and $\rho(x,y,t)$ are the plate thickness and material density functions respectively. $m_{i,j}$ are the bending moment stress resultants, which read:

[TABLE]

where $E(x,y,t)$ is the Young’s Modulus and $B(x,y,t)=\frac{E(x,y,t)}{1-\nu^{2}}\frac{s^{3}(x,y)}{12}$ is the bending stiffness, which is a general function of space and time. Merging Eqs. 1 and 2 gives the PDE governing elastic wave propagation in modulated plates, reported in Appendix A for the sake of brevity. Now, under the assumption of space-time periodic elasticity and material density, $B(x,y,t)$ and $G\left(x,y,t\right)=s(x,y)\rho(x,y,t)$ can be written in terms of exponential functions:

[TABLE]

where the term $\sum\limits_{m,n,v=-\infty}^{\infty}$ denotes a nested summation over the indexes $m,n,v$. $\bm{\kappa_{m}}=\left(mk_{mx},nk_{my}\right)$ and $\bm{r}=\left(x,y\right)$ are the modulation wavevector and spatial coordinates mapped within the unit cell domain. In the case at hand, we define $\kappa_{mx}=\frac{2\pi}{\lambda_{mx}}$, $\kappa_{my}=\frac{2\pi}{\lambda_{my}}$ and $\omega_{m}=\frac{2\pi}{T_{m}}$ as spatial and temporal modulation wavenumbers and frequency, with associated wavelengths $\lambda_{mx}$ and $\lambda_{my}$ and temporal period $T_{m}$. Corresponding Fourier coefficients $\hat{B}_{m,n,v}$ and $\hat{G}_{m,n,v}$ are obtained by numerical integration within the domain $D=\left[-\displaystyle\frac{\lambda_{mx}}{2},\displaystyle\frac{\lambda_{mx}}{2}\right]\times\left[-\displaystyle\frac{\lambda_{my}}{2},\frac{\lambda_{my}}{2}\right]\times\left[-\displaystyle\frac{T_{m}}{2},\frac{T_{m}}{2}\right]$:

[TABLE]

A general formulation of the PWEM is thus presented here to compute the wave propagation properties associated with the spatiotemporal unitary cell, therefore assuming a solution $w(x,y,t)$ which is a combination of plane waves owning the same periodicity of the modulation:

[TABLE]

where $\bm{\kappa}=\left(\kappa_{x},\kappa_{y}\right)$ is the imposed wavevector field and $\bm{\kappa_{m}}=\left(pk_{mx},qk_{my}\right)$, differently from already existing - and less general - formulations in which the wave solution is reconstructed using plane waves having the same propagation speed along $x$ and $y$ directions. Indeed, thanks to the 3D Fourier transform of the unit cell profile, $B(x,y,t)$ and $G(x,y,t)$ are written as a series of traveling plane waves regardless the space-time function designed within the domain $D$ [7]. Upon combination of Eqs. 1-5 and enforcing orthogonality of exponential functions, the dispersion relation $\omega=\omega\left(\kappa_{x},\kappa_{y}\right)$ yields:

[TABLE]

which is a Quadratic Eigenvalue Problem (QEP) solved for $\omega$ imposing $\kappa_{x}$ and $\kappa_{y}$. $\bm{\tilde{L}}_{0}(k_{x},k_{y})$, $\bm{\tilde{L}}_{1}$ and $\bm{\tilde{L}}_{2}$ are full matrices of dimension $\Gamma_{o}=(2P+1)(2Q+1)(2R+1)$, being $P,Q,R$ are the Fourier expansion truncation orders. A complete description of the analytical procedure and QEP matrices is provided in Appendix A.

3 Wave propagation in discretely modulated plates

Consider a periodic plate whose elasticity is modulated according with the following law:

[TABLE]

where $R_{s}$ denotes the number of sub-cells along $x$ and $y$, $\alpha_{m}=E_{m}/E_{0}$ is the dimensionless modulation amplitude, thus $E_{m}$ and $E_{0}$ are the maximum and mean Young’s Modulus values. $i,j$ define the temporal phase shift between consecutive sub-cells along $x$ and $y$ directions respectively. In the case at hand, the modulation is designed to mimic a plane wave that propagates along $\pi/4$ for $R_{s}=3$, as shown in Fig. 1 for three consecutive time instants.

Plugging the Fourier expansion of $E\left(x,y,t\right)$ for $\nu=0.02$, $\alpha_{m}=0.8$ and $\lambda_{mx}=\lambda_{my}=\lambda_{m}=0.06\;m$ into Eq. 6 gives the associated dispersion relation, which is represented in Figs. 2-2 mapped along different wave propagation directions $\gamma$. Here $\Omega=\omega/\left(c_{0}\kappa_{m}\right)$, being $c_{0}=\sqrt{E_{0}/\rho}$ and:

[TABLE]

where $\gamma=\arctan{\left(\kappa_{y}/\kappa_{x}\right)}$ and $\mu_{x}=\kappa_{x}\lambda_{mx}$. Given the time-space periodicity of $w\left(x,y,t\right)$, the number of dispersion branches is consistent with the harmonics considered in the analysis, thus $P=3,Q=3,R=1$ and the series expansion is limited to $\left(2P+1\right)\left(2Q+1\right)$ spatial and $\left(2R+1\right)$ temporal harmonics. The leading terms are thus identified by the magnitude of the associated eigenvector components [15] which is shown by the colored plot in Figs. 2-2, therefore revealing the presence of non-reciprocal bandgaps. In the remainder of this section we restrict the analysis to these components, thus neglecting higher order harmonics, as waves propagate mainly as their leading terms. This assumption is confirmed from numerical simulations performed using the commercial software COMSOL Multiphysics. Specifically, the system is made of $90\times 90$ unit cells and forced in its central region using a wide spectrum tone burst excitation. The resulting displacement field is integrated within the spatial and temporal simulation domains to obtain the numerical Bloch diagram, proving good agreement with the analytical results, as shown in the comparison in Figs. 2-2 for $\gamma=0,\pi/4,3/4\pi$ respectively.

A more complete description of the wave propagation properties is presented in Figs. 3-3, in terms of directivity plot: $\Omega\left(\kappa_{x},\kappa_{y}\right)$ is mapped in a polar plane, thus $\Omega$ is the distance from the origin and $\gamma=\arctan{\left(\kappa_{y}/\kappa_{x}\right)}$ represents the angle. Specifically, for unmodulated, spatially and spatiotemporally modulated plates. When the spatial modulation is turned on, a Bragg-bandgap opens, whose position in the frequency domain varies as a function of the direction of interest $\gamma$. The addition of a temporal periodicity results into a modulation-induced tilting of the frequency-momentum space, which is evident from the asymmetric plot in Fig. 3. The amount of frequency bias is dictated by the scalar product between modulation velocity vector and wave propagation direction under investigation, thus:

[TABLE]

In case of $\gamma=\pi/4$, the velocity component is $\nu_{\pi/4}=\nu_{max}$, while $\nu_{3/4\pi}=0$. As a result, band diagram is tilted along every direction but $3/4\pi+i\pi$, as already pointed out in Ref. [16] for spatiotemporal membranes. Consider now the non-reciprocal wave propagation problem between two points A (emitter) and B (receiver). Point A is located in the central region of the plate, whereas B must be set sufficiently far from the excitation. In contrast with 1D waveguides, when the system is forced waves propagate along the direction identified by the group velocity field $C_{g}=\nabla\Omega\left(\mu_{x},\mu_{y}\right)$, which is shown in Fig. 3 for unmodulated, spatially and spatiotemporally modulated medium.

Thanks to the spatiotemporal modulation, this latter is characterized by angular regions in which wave propagation does not occur. A non-reciprocal device can be thus achieved when the receiver B is located within $\beta_{1}\in\left[160^{\circ},200^{\circ}\right]$ and $\beta_{2}\in\left[250^{\circ},290^{\circ}\right]$ for $\Omega\approx 0.065$, which is confirmed by the RMS of the displacement field obtained by numerical simulation under narrowband spectrum excitation as shown in Fig. 4. The resulting input/output transmitted signal is clearly asymmetric, proving that the phononic plate is able to break time reversal symmetry along the predicted directions in Fig. 3. Mirrored stop-regions within $\left[-20^{\circ},20^{\circ}\right]$ and $\left[70^{\circ},110^{\circ}\right]$ can be achieved inverting the modulation direction.

A similar behavior is also achieved along an angular region centered at $45^{\circ}$ for $\Omega\approx 0.14$, i.e. along the direction of propagation of the modulation, which is tailorable, as shown by 3 and 4 where corresponding group velocity and RMS vibrational content are shown.

4 Conclusions

In this work we presented a general formulation of the PWEM in order to study non-reciprocal wave propagation in spatiotemporally modulated plates. This analysis tool is applicable to any class of modulations that can be written as a series of traveling plane waves, thus extending the allowable modulation classes to any type unit cell profile. Specifically, not limiting the analysis to continuous plane wave modulations, which are difficult to be realized for practical applications.

In the second part of the manuscript, the extended PWEM has been applied to spatially discrete and temporally continuous elastic medium. Numerical solutions have been compared to theoretical results, proving that the PWEM is able to correctly describe wave propagation properties of arbitrary shaped unit cell profiles. Specifically, we demonstrated that discretely modulated plates are able to break the reciprocity principle within different angular regions. The proposed configuration has been tested in time domain as a non-reciprocal device, which is able to support one-way elastic signals between input-output ports.

Appendix A QEP matrices description

Consider the wave propagation problem described in Section 2. The general equation governing the out of plane dynamic for a Kirchoff plate reads:

[TABLE]

Plugging Eqs. 3-5 into Eq. 10 gives:

[TABLE]

which can be simplified exploiting the orthogonality of the Fourier basis, thus all the terms are multiplied by ${\rm{e}}^{{-\rm{j}}(ak_{mx}x+bk_{my}y-c\omega_{m}t)}$ and integrated over $D=\left[-\frac{\lambda_{mx}}{2},\frac{\lambda_{mx}}{2}\right]\times\left[-\frac{\lambda_{my}}{2},\frac{\lambda_{my}}{2}\right]\times\left[-\frac{T_{m}}{2},\frac{T_{m}}{2}\right]$. Eq. (LABEL:eq18), can now be conveniently rewritten as:

[TABLE]

Eq. 12 is a QEP, which can be written using a compact matrix notation, expanding the inner summation terms:

[TABLE]

where $\bm{\tilde{K}}_{a-p,b-q}$, $\bm{\tilde{M}}^{0}_{a-p,b-q}$, $\bm{\tilde{M}}^{1}_{a-p,b-q}$, $\bm{\tilde{M}}^{2}_{a-p,b-q}$ are full square matrices of size $(2R+1)$:

[TABLE]

and $\bm{\tilde{w}}_{p,q}$ accommodates the $(2R+1)$ time-harmonic components:

[TABLE]

In the same way, the summation term $q\in[-Q,Q]$ can be expanded, thus:

[TABLE]

where $\bm{\tilde{Z}}^{0}_{a-p}$, $\bm{\tilde{Z}}^{1}_{a-p}$, $\bm{\tilde{Z}}^{2}_{a-p}$ are full square matrices of order $(2Q+1)(2R+1)$ which take the following form:

[TABLE]

which is further expanded for $p\in[-P,P]$ wave components:

[TABLE]

where $\bm{\tilde{L}}_{0}(k_{x},k_{y})$, $\bm{\tilde{L}}_{1}$, $\bm{\tilde{L}}_{2}$ are full square matrices of order

$(2P+1)(2Q+1)(2R+1)$:

[TABLE]