Anomalous energy transport in laminates with exceptional points

TL;DR

This paper investigates unusual energy transport phenomena in elastic laminates, revealing modes that attenuate and propagate simultaneously, and demonstrating exceptional points leading to metamaterial effects like negative refraction and beam steering.

Contribution

It uncovers the existence of exceptional points in in-plane elastic waves in laminates, enabling novel metamaterial behaviors not seen in anti-plane shear wave models.

Findings

Modes that attenuate and propagate simultaneously

Existence of exceptional points where modes coalesce

Negative refraction and beam steering observed near these points

Abstract

Recent interest in metamaterials has led to a renewed study of wave mechanics in different branches of physics. Elastodynamics involves a special intricacy, owing to a coupling between the volumetric and shear parts of the elastic waves. Through a study of in-plane waves traversing periodic laminates, we here show that this coupling can result with unusual energy transport. We find that the corresponding frequency spectrum contains modes which simultaneously attenuate and propagate, and demonstrate that these modes coalesce to purely propagating modes at exceptional points, a property that was recently reported in parity-time symmetric systems. We show that the laminate exhibits metamaterial features near these points, such as negative refraction, and beam steering and splitting. While negative refraction in laminates has been demonstrated before by considering pure shear waves impinging an interface with multiple layers, here we realize it for coupled waves impinging a simple single-layer interface. This feature, together with the appearance of exceptional points, are absent from the model problem of anti-plane shear waves which have no volumetric part, and hence from the mathematically identical electromagnetic waves. Our work further paves the way for applications such as asymmetric mode switches, by encircling exceptional points in a tangible, purely elastic apparatus.

Equations

$$\begin{array}{r}\left(\lambda {}^{\left(\text{p}\right)}+2\mu {}^{\left(\text{p}\right)}\right)u_{1,11}{}^{\left(n\right)}+\left(\lambda {}^{\left(\text{p}\right)}+\mu {}^{\left(\text{p}\right)}\right)u_{2,21}{}^{\left(n\right)}+\mu {}^{\left(\text{p}\right)}{u}_{1,22}{}^{\left(n\right)}=\rho {}^{\left(\text{p}\right)}\ddot{u_1}{}^{\left(n\right)},\\ \left(\lambda {}^{\left(\text{p}\right)}+2\mu {}^{\left(\text{p}\right)}\right)u_{2,22}{}^{\left(n\right)}+\left(\lambda {}^{\left(\text{p}\right)}+\mu {}^{\left(\text{p}\right)}\right)u_{1,12}{}^{\left(n\right)}+\mu {}^{\left(\text{p}\right)}{u}_{2,11}{}^{\left(n\right)}=\rho {}^{\left(\text{p}\right)}\ddot{u_2}{}^{\left(n\right)},\end{array}$$

$$\begin{array}{rl}{u}_1{}^{\left(n\right)}={\varphi }_{,1}{}^{\left(n\right)}+{\psi }_{,2}{}^{\left(n\right)},& u_2{}^{\left(n\right)}={\varphi }_{,2}{}^{\left(n\right)}-{\psi }_{,1}{}^{\left(n\right)}\end{array}.$$

$${\nabla }^2\varphi {}^{\left(n\right)}=\frac{1}{c_L^{2\left(\text{p}\right)}}\ddot{\varphi }{}^{\left(n\right)},$$

$$\varphi {}^{\left(n\right)}=A_{L+}^{\left(\text{n}\right)}{e}^{i\left(\omega t+k_L^{\left(\text{p}\right)}{x}_1-k_2{}^{\left(\text{p}\right)}{x}_2\right)}+A_{L-}^{\left(\text{n}\right)}{e}^{i\left(\omega t-k_L^{\left(\text{p}\right)}{x}_1-k_2{}^{\left(\text{p}\right)}{x}_2\right)},$$

$$\psi {}^{\left(n\right)}=A_{S+}^{\left(\text{n}\right)}{e}^{i\left(\omega t+k_S^{\left(p\right)}{x}_1-k_2{}^{\left(\text{p}\right)}{x}_2\right)}+A_{S-}^{\left(\text{n}\right)}{e}^{i\left(\omega t-k_S^{\left(p\right)}{x}_1-k_2{}^{\left(\text{p}\right)}{x}_2\right)},$$

$$k_L^{2\left(\text{p}\right)}=\frac{{\omega }^2}{c_L^{2\left(\text{p}\right)}}-k_2^2{}^{\left(\text{p}\right)},$$

$$\mathsf{s}{}^{\left(n\right)}={\left(u_1{}^{\left(n\right)},{\sigma }_{11}{}^{\left(n\right)},u_2{}^{\left(n\right)},{\sigma }_{21}{}^{\left(n\right)}\right)}^{\mathsf{T}};$$

$$\mathsf{s}{}^{\left(n\right)}\left(x_1=x_0{}^{\left(n\right)}-h{}^{\left(n\right)}\right)=e^{i{k}_{\mathrm{B}}h}{\mathsf{s}}^{\left(n+1\right)}\left(x_1=x_0{}^{\left(n\right)}+h^{\left(n+1\right)}\right).$$

$$\left(\begin{array}[]{ll}-\mathsf{I}_{2}&\mathsf{H}_{11}\\ \mathsf{0}_{2}&\mathsf{H}_{21}\end{array}\right)\cdot\mathsf{s_{m}}\left(x_{1},x_{2}\right)=e^{ik_{B}h}\left(\begin{array}[]{ll}-\mathsf{H}_{12}&\mathsf{0}_{2}\\ -\mathsf{H}_{22}&\mathsf{I}_{2}\end{array}\right)\cdot\mathsf{s_{m}}\left(x_{1},x_{2}\right),$$

$$\cos k_{\mathrm{B}}h=\frac{-a_2\pm \sqrt{a_2^2-4a_1\left(a_3-2a_1\right)}}{4a_1},$$

$${\left\{\left(u_{i,j}+u_{j,i}\right)\mu \left(x_1\right)+\lambda \left(x_1\right){\delta }_{ij}{u}_{k,k}\right\}}_{,j}=\rho \left(x_1\right){\ddot{u}}_i,$$

$$u_i={\tilde{u}}_i\left(x_1\right)e^{i\left(\omega t-k_{\mathrm{B}}{x}_1-k_2{x}_2\right)},{\tilde{u}}_i\left(x_1\right)={\tilde{u}}_i\left(x_1+h\right),$$

$$\zeta \left(x_1\right)=\sum _{\mathrm{G}}\zeta \left(\mathrm{G}\right)e^{-i\mathrm{G}{x}_1},\zeta =\rho ,\mu ,\lambda ,{\tilde{u}}_1,{\tilde{u}}_2,$$

$$\zeta \left(\mathrm{G}\right)=\frac{1}{h}{\int }_0^h\zeta \left(x_1\right)e^{i\mathrm{G}{x}_1}\mathrm{d}{x}_1.$$

$$\zeta \left(\mathrm{G}\right)=\{\begin{array}{ll}\frac{h^{\left(\text{a}\right)}}{h}{\zeta }^{\left(\text{a}\right)}+\frac{h^{\left(\text{b}\right)}}{h}{\zeta }^{\left(\text{b}\right)},& \mathrm{G}=0,\\ {\left(\pi m\right)}^{-1}\left({\zeta }^{\left(\text{a}\right)}-{\zeta }^{\left(\text{b}\right)}\right)\left(1-e^{i\mathrm{G}{h}^{\left(\text{a}\right)}}\right),& \mathrm{G}\ne 0.\end{array}$$

$$\begin{array}{rl}& \sum _{\mathrm{G},{\mathrm{G}}^{\prime }}\{-k_2\left[k_2{\tilde{u}}_1\left({\mathrm{G}}^{\prime }\right)+\left({\mathrm{G}}^{\prime }+k_{\mathrm{B}}\right){\tilde{u}}_2\left({\mathrm{G}}^{\prime }\right)\right]\mu \left(\mathrm{G}\right)-\\ & -\left(\mathrm{G}+{\mathrm{G}}^{\prime }+k_{\mathrm{B}}\right)\left[k_2\lambda \left(\mathrm{G}\right){\tilde{u}}_2\left({\mathrm{G}}^{\prime }\right)+\left({\mathrm{G}}^{\prime }+k_{\mathrm{B}}\right)\left(\lambda \left(\mathrm{G}\right)+2\mu \left(\mathrm{G}\right)\right){\tilde{u}}_1\left({\mathrm{G}}^{\prime }\right)\right]\}{e}^{-i\left(\mathrm{G}+{\mathrm{G}}^{\prime }\right)x_1}=\\ & =-\sum _{\mathrm{G},{\mathrm{G}}^{\prime }}{\omega }^2\rho \left(\mathrm{G}\right){\tilde{u}}_1\left({\mathrm{G}}^{\prime }\right)e^{-i\left(\mathrm{G}+{\mathrm{G}}^{\prime }\right)x_1},\\ & \sum _{\mathrm{G},{\mathrm{G}}^{\prime }}\{-k_2\left[\left({\mathrm{G}}^{\prime }+k_{\mathrm{B}}\right)\lambda \left(G\right)+\left(\mathrm{G}+{\mathrm{G}}^{\prime }+k_{\mathrm{B}}\right)\mu \left(\mathrm{G}\right)\right]{\tilde{u}}_1\left({\mathrm{G}}^{\prime }\right)-k_2^2\lambda \left(\mathrm{G}\right){\tilde{u}}_2\left({\mathrm{G}}^{\prime }\right)-\\ & -\mu \left(\mathrm{G}\right)\left[2k_2^2+\left({\mathrm{G}}^{\prime }+k_{\mathrm{B}}\right)\left(G+{\mathrm{G}}^{\prime }+k_{\mathrm{B}}\right)\right]{\tilde{u}}_2\left({\mathrm{G}}^{\prime }\right)\}{e}^{-i\left(\mathrm{G}+{\mathrm{G}}^{\prime }\right)x_1}=\\ & =-\sum _{\mathrm{G},{\mathrm{G}}^{\prime }}{\omega }^2\rho \left(\mathrm{G}\right){\tilde{u}}_2\left({\mathrm{G}}^{\prime }\right)e^{-i\left(\mathrm{G}+{\mathrm{G}}^{\prime }\right)x_1}.\end{array}$$

$$\left[{\mathsf{A}}_0+k_2{\mathsf{A}}_1+k_2^2{\mathsf{A}}_2\right]\tilde{\mathsf{u}}\left({\mathrm{G}}^{\prime }\right)={\omega }^2\mathsf{B}\tilde{\mathsf{u}}\left({\mathrm{G}}^{\prime }\right),$$

$$\left[\begin{array}[]{ll}\mathsf{0}&\mathsf{I}\\ \omega^{2}\mathsf{B}-\mathsf{A}_{0}&-\mathsf{A}_{1}\end{array}\right]\left(\begin{array}[]{l}\mathsf{\tilde{u}}\left(\mathrm{G^{\prime}}\right)\\ k_{2}\mathsf{\tilde{u}}\left(\mathrm{G^{\prime}}\right)\end{array}\right)=k_{2}\left[\begin{array}[]{ll}\mathsf{I}&\mathsf{0}\\ \mathsf{0}&\mathsf{A}_{2}\end{array}\right]\left(\begin{array}[]{l}\mathsf{\tilde{u}}\left(\mathrm{G^{\prime}}\right)\\ k_{2}\mathsf{\tilde{u}}\left(\mathrm{G^{\prime}}\right)\end{array}\right)$$

$$⟨\mathcal{P}⟩=-\frac{1}{2}⟨\mathrm{Re}\left[\mathbf{\sigma }{\dot{\mathbf{u}}}^{\ast }\right]⟩,$$

$$⟨\mathcal{P}⟩=\frac{\partial \omega }{\partial \mathbf{k}}E,E=\frac{1}{4}⟨{\sigma }_{ij}{u}_{i,j}^{\ast }⟩+\frac{1}{4}{\omega }^2⟨\rho u_k^{\ast }{u}_k⟩\equiv \frac{1}{2}{\omega }^2⟨\rho u_k^{\ast }{u}_k⟩,$$

$$\theta =\arctan \frac{⟨{\mathcal{P}}_1⟩}{⟨{\mathcal{P}}_2⟩}=\arctan \frac{\partial \omega /\partial k_{\mathrm{B}}}{\partial \omega /\partial k_2}.$$

$${\mathcal{P}}_2=\frac{1}{2}\omega \mathrm{Re}\left\{\mu \left(k_2{\tilde{u}}_1+k_{\mathrm{B}}{\tilde{u}}_2+i{\tilde{u}}_2^{\prime }\right){\tilde{u}}_1^{\ast }+\left(k_{\mathrm{B}}\lambda {\tilde{u}}_1+k_2\left(\lambda +2\mu \right){\tilde{u}}_2+i\lambda {\tilde{u}}_1^{\prime }\right){\tilde{u}}_2^{\ast }\right\},{\left(\cdot \right)}^{\prime }\equiv \frac{\partial \left(\cdot \right)}{\partial x_1},$$

$${\mathcal{P}}_1=\frac{1}{2}\omega \mathrm{Re}\left\{\mu \left(k_2{\tilde{u}}_1+k_{\mathrm{B}}{\tilde{u}}_2+i{\tilde{u}}_2^{\prime }\right){\tilde{u}}_2^{\ast }+\left(\lambda +2\mu \right)\left(k_{\mathrm{B}}{\tilde{u}}_1+i{\tilde{u}}_1^{\prime }\right){\tilde{u}}_1^{\ast }+\lambda k_2{\tilde{u}}_2{\tilde{u}}_1^{\ast }\right\}.$$

$$\begin{array}[]{ll}\mu^{\left(1\right)}=78.85\,\mathrm{GPa},\lambda^{\left(1\right)}=118.27\,\mathrm{GPa},&\rho^{\left(1\right)}=7800\,\mathrm{kg\,m^{-3}},h^{\left(1\right)}=1.3\,\mathrm{mm},\\ \mu^{\left(2\right)}=1.21\,\mathrm{GPa},\ \lambda^{\left(2\right)}=4.86\,\mathrm{GPa},&\rho^{\left(2\right)}=1200\,\mathrm{kg\,m^{-3}},h^{\left(2\right)}=3\,\mathrm{mm}.\end{array}$$

$$V=I{e}^{i\left(\omega t-k\sin {\theta }_i{x}_1-k\cos {\theta }_i{x}_2\right)},$$

$$\mathbf{T}=\sum _{m=1}^2{T}_m{\tilde{\mathbf{u}}}^{\left(m\right)}\left(x_1\right)e^{i\left(\omega t-k_{\mathrm{B}}{x}_1-k\cos {\theta }_i{x}_2\right)},$$

$${\varphi }_R^{\left(0\right)}=R^{\varphi }{e}^{i\left(\omega t+k_L^{\left(0\right)}{x}_1-k\cos {\theta }_i{x}_2\right)},{\psi }_R^{\left(0\right)}=R^{\psi }{e}^{i\left(\omega t+k_S^{\left(0\right)}{x}_1-k\cos {\theta }_i{x}_2\right)},$$

$${\mu }^{\left(0\right)}=0.178\mathrm{GPa},{\lambda }^{\left(0\right)}=0.714\mathrm{GPa},{\rho }^{\left(0\right)}=3000kg/{\mathrm{m}}^3.$$

$$\frac{1}{{\left|I\right|}^2{⟨{\mathcal{P}}_1⟩}_I}\left[\sum _{m=1}^2{\left|T_m\right|}^2{⟨{\mathcal{P}}_1⟩}_{T_m}^{\left(1\right)}-{\left|R^{\varphi }\right|}^2{⟨{\mathcal{P}}_1⟩}_{\varphi }-{\left|R^{\psi }\right|}^2{⟨{\mathcal{P}}_1⟩}_{\psi }\right]=1,$$

$$\mathbf{T}=\sum _{m=1}^{\infty }{T}_m{\tilde{\mathbf{u}}}^{\left(m\right)}\left(x_1\right)e^{i\left(\omega t-k_{\mathrm{B}}{x}_1-k_m{x}_2\right)},$$

$${\varphi }_R=\sum _{l=-\infty }^{\infty }{R}_l^{\varphi }{U}_l\left(x_1\right)e^{i\left(\omega t-k_{\mathrm{B}}{x}_1+{\kappa }_l^{\varphi }{x}_2\right)},{\psi }_R=\sum _{l=-\infty }^{\infty }{R}_l^{\psi }{U}_l\left(x_1\right)e^{i\left(\omega t-k_{\mathrm{B}}{x}_1+{\kappa }_l^{\psi }{x}_2\right)},$$

Full text

Keywords: metamaterial, negative refraction, wave propagation, Bloch-Floquet waves, composite, laminate, phononic crystal, exceptional points

1 Introduction

Metamaterials possess properties not found in nature, stemming from their architectured microstructure (Wegener, 2013, Kadic et al., 2019). Perhaps the most prominent thrust in metamaterials research aims at controlling waves for potential applications such as lensing, cloaking and noise reduction (Milton et al., 2006, Chen et al., 2010, Parnell and Shearer, 2013, Bigoni et al., 2013, Colquitt et al., 2014, Cummer et al., 2016). This interest led to a renewed study of wave mechanics in optics (Markos and Soukoulis, 2008, Banerjee, 2011), acoustics (Craster and Guenneau, 2012, Deymier, 2013), mechanical lattices (Phani, 2011, Raney et al., 2016, Phani and Hussein, 2017, Zelhofer and Kochmann, 2017, Ma et al., 2018) and elastodynamics (Brun et al., 2010, Shmuel and Band, 2016, Chen and Elbanna, 2017, Aghighi et al., 2019, Li and Reina, 2019).

While electromagnetic-, sound- and anti-plane shear waves are mathematically identical (Adams et al., 2008, Torrent and Sánchez-Dehesa, 2011), in-plane elastic waves are physically richer (Sigalas and Economou, 1992), since they comprise both volumetric and distortional parts, coupled through interfaces in the transmission medium. In this work, we show that in the simplest elastic composite—a laminate—this coupling gives rise to anomalous energy transport, which in other systems is achieved by significantly more complicated means. The anomalies reported here and their connection with recent studies in the field are summarized next.

Firstly, we show that the spectrum of in-plane waves in laminates exhibits exceptional points, accessible in a purely elastic setting. Exceptional points are states of a system at which two (or more) of its normal modes coalesce, together with their natural frequencies (Moiseyev and Friedland, 1980, Ding et al., 2015), and are the source of counterintuitive phenomena such as enhanced sensitivity, wave stopping and asymmetric transmission (Hodaei et al., 2017, Achilleos et al., 2017, Goldzak et al., 2018, Merkel et al., 2018). Exceptional points occur only in non-Hermitian systems (Moiseyev, 2011), a property which usually describes systems that interact with the environment. The current paradigm111To break reciprocity, we note that another emerging paradigm is to employ spatiotemporal composites (Trainiti and Ruzzene, 2016, Nassar et al., 2017, Milton and Mattei, 2017) to access exceptional points is by balancing external gain and loss through the system, to create parity-time ($\mathcal{PT}$) symmetry (Rüter et al., 2010, Shi et al., 2016, El-Ganainy et al., 2018). This procedure brings with additional complexity to the elastic medium, as such realizations require incorporating optomechanical, acoustoelectric or piezoelectric elements into the system (Xu et al., 2015, Christensen et al., 2016, Hou and Assouar, 2018). Our findings thus suggest a simpler, purely elastic setting to access exceptional points, thereby evading these complexities.

Recently, it was demonstrated that encircling these points leads to fascinating asymmetric mode switching of microwaves in a metallic waveguide (Doppler et al., 2016). We argue that our setting constitutes a tangible platform to realize analogous encirclement for elastic wave switching, having the wavenumber, which is related to the excitation angle, as the trajectory parameter. As we show in the sequel, this is made possible owing to the intrinsic Riemann surface structure of our spectrum near these points in the complex wave vector space. In contrast with cases where this structure is artificially obtained by an analytical continuation of an arbitrary parameter (e.g., the frequency in Shanin et al., 2018), here, complex wavenumbers are an intrinsic and accessible part of the spectrum.

Secondly, we demonstrate that the exceptional points foreshadow anomalous energy transport, including negative refraction. In this regard, perhaps our most striking result is the excitation of negative refraction in laminates by an incoming wave from a homogeneous medium whose interface with the laminate is parallel to the layers (Fig.  1). To put this result into context, we recall that the first report of negative refraction was for electromagnetic waves, and required a two-dimensional composite made of complicated split ring resonators (Smith et al., 2000, Shelby et al., 2001, Pendry, 2004). Interest in this phenomenon has been disseminated to elastodynamics, accompanied with ongoing studies (Craster and Guenneau, 2012, Chen et al., 2017, Bordiga et al., 2019, Nemat-Nasser, 2019). Notably, Willis (2013) was the first to show that a simple laminate is capable of negatively refracting anti-plane shear waves. While in the conventional arrangement (interface parallel to the layers), refraction is always positive, Willis realized that negative refraction is possible when the interface is normal to the layers (Fig.  1c). Based on this interface configuration, further studies of such waves in laminated media were carried out (Willis, 2015, Nemat-Nasser, 2015, Srivastava, 2016, Morini et al., 2019). Hence, we demonstrate that in-plane waves may refract negatively in the conventional arrangement, without the need for this complex apparatus, nor for gain and loss (cf.  Hou et al., 2018).

For completeness, we also analyze the transmission problem through an interface normal to the layers. As highlighted by Srivastava and Willis (2017) and the references therein, a wave incident to such interface induces an infinite number of transmitted waves, which are required to satisfy corresponding continuity conditions across the interface. We suggest a method to calculate the resultant normal mode decomposition owing to incident in-plane waves, based on suitable orthogonality conditions (Mokhtari et al., 2019). Subsequently, we demonstrate analogous phenomena to those reported by Srivastava (2016), who studied the anti-plane problem. These include beam steering—small changes in the incident angle leading to large changes in the transmission angle; beam splitting—an incident wave transmitted as simultaneous negative and positive refracted beams; and pure negative refraction.

Our study is detailed in the forthcomings Secs.  as follows. Sec.  2 firstly revisits the equations governing in-plane waves in infinite periodic laminates, and subsequently provides two formulations—the extended plane wave expansion method and the hybrid matrix method—for their solution (Laude et al., 2009, Tan, 2010). By analyzing the pertinent equations, we explain why in contrast with the anti-plane problem, we here obtain complex wavenumbers, exceptional points, and negative refraction in the simple configuration. We numerically solve the eigenvalue problems for an exemplary infinite laminate, and present its spectrum in Sec.  3. As predicted, the complex spectrum contains exceptional points with hallmarks of negative refraction. Sec.  4 links the studied spectrum and eigenmodes to their excitation by incoming waves from a homogeneous half-space that shares an interface with the exemplary laminate. Finally, our main results and conclusions are summarized in Sec. 5.

2 Equations of in-plane waves in laminated media

The equations governing in-plane waves in elastic solids can be found in Graff (1975), and their extension to laminated media appears, e.g., in Lowe (1995) and the references therein (see also Adams et al., 2009 for in-plane Bloch waves in infinitely periodic strips). In this Sec., we firstly revisit the equations for an infinite laminate comprising two alternating layers, and provide two formulations to determine the resultant normal modes via different eigenvalue problems. Based on these formulations, we explain next why in contrast with anti-plane waves, in-plane waves in laminates admit complex eigenvalues and exceptional points, and connect these features of the spectrum to the transport of energy.

2.1 Governing Equations

We study an infinite elastic laminate made of a periodic repetition of phases a and b in the $x_{1}$ direction (Fig. 1a). The thickness, mass density and Lamé coefficients of each phase are $h^{\left(\emph{p}\right)},\rho^{\left(\emph{p}\right)},\lambda^{\left(\emph{p}\right)}$ and $\mu^{\left(\emph{p}\right)}$, respectively, where $\emph{p}=\emph{a}$ or b, denoting the corresponding phase.

\floatsetup

[figure]style=plain,subcapbesideposition=top

The objective is to determine the propagation of free time-harmonic in-plane waves in the laminate. To this end, we seek solutions for the displacements $u_{1}{}^{\left(n\right)}$ and $u_{2}{}^{\left(n\right)}$ in each $n^{\mathrm{th}}$ layer using the Naiver-Lamé equations

[TABLE]

subjected to boundary conditions that will be specified later. We employ the Helmholtz decomposition to write the in-plane components of $\bm{u}{}^{\left(n\right)}$ in terms of scalar potential $\phi{}^{\left(n\right)}$ and $\psi{}^{\left(n\right)}$, namely,

[TABLE]

This decomposition simplifies Eq.  (1) to the form

[TABLE]

where $c_{L}^{\left(\emph{p}\right)}=\sqrt{\left(\lambda^{\left(p\right)}+2\mu^{\left(p\right)}\right)/\rho^{\left(p\right)}}$ and $c_{S}^{\left(\emph{p}\right)}=\sqrt{\mu^{\left(p\right)}/\rho^{\left(p\right)}}$ are the pressure and shear wave velocities of phase p, respectively. At a fixed frequency $\omega$, Eq. (3) is solved by

[TABLE]

where

[TABLE]

and $\left\{A_{L/S,\pm}{}^{\left(n\right)}\right\}$ are integration constants to be determined from the conditions on the boundaries of layer $n$. These correspond to the continuity of the traction and displacements at interfaces between adjacent layers, which immediately requires $k_{2}^{\left(\emph{a}\right)}=k_{2}^{\left(\emph{b}\right)}\eqqcolon k_{2}$. The corresponding equations are compactly written in terms of the state vector $\mathsf{s}{}^{\left(n\right)}\left(x_{1},x_{2},t\right)$

[TABLE]

say the boundary between layers $n$ and $n+1$ is at $x_{1}=x_{0}^{\left(n\right)}$, then the continuity conditions are simply $\mathsf{s}{}^{\left(n\right)}\left(x_{1}=x_{0}{}^{\left(n\right)}\right)=\mathsf{s}^{\left(n+1\right)}\left(x_{1}=x_{0}{}^{\left(n\right)}\right)$. The coupling between shear and pressure modes enters through the latter condition, since the components of $\mathsf{s}{}^{\left(n\right)}$ depend both on $\phi{}^{\left(n\right)}$ and $\psi{}^{\left(n\right)}$. The remaining equations stem from the Bloch-Floquet theorem, which states that over the course of one period the governing fields are related via the Bloch wavenumber $k_{\mathrm{B}}$, namely,

[TABLE]

2.2 The hybrid matrix method to determine $k_{\mathrm{B}}\left(\omega\right)$

Eqs. (6) and (7) can be combined in different ways to deliver an eigenproblem for $k_{\mathrm{B}}$ and $\mathsf{s}{}^{\left(n\right)}$ as functions of real $\omega$ at prescribed $k_{2}$, given the laminate composition. The transfer matrix formulation is the most intuitive and common approach, however it suffers from numerical instabilities (Pérez-Álvarez et al., 2015). Here, we employ the stable hybrid matrix method (Tan, 2010), which reads

[TABLE]

where $\mathsf{0}_{2}$ (resp. $\mathsf{I}_{2}$) is the $2\times 2$ zeros (resp. unit) matrix, and $\mathsf{s_{m}}$ and $\mathsf{H}_{ij}$ are given in Appendix A, together with a detailed derivation of Eq. (8). The resultant quartic equation for $k_{\mathrm{B}}$ is solved by

[TABLE]

where $a_{i}$ are given in Appendix A. Eq.  (9) provides the dispersion relation which relates the microstructure and mechanical properties of the laminate to the waveform at each frequency. The structure of Eq.  (9) implies that if $k_{\mathrm{B}}$ is a solution, then so are $\pm k_{\mathrm{B}}+2\pi m$ for any $m\in\mathbb{Z}$. Hence, the real part of all dispersion curves is representable over the irreducible $1^{\mathrm{st}}$ Brillouin zone222See, *e.g., *Zhang (2019) for a discussion on the reciprocal space symmetries and degeneracies in the general case. $0<k_{\mathrm{B}}<\pi/h$. The wavenumber $k_{\mathrm{B}}$ can be real, pure imaginary (henceforth referred as imaginary) or complex, as demonstrated in the sequel, and contrary to the case of anti-plane shear (Willis, 2015, Srivastava, 2016). This property is essential for the spectrum to exhibit a Riemann structure, and, in turn, exceptional points. Real $k_{\mathrm{B}}$ corresponds to a propagating Bloch mode along $x_{1}$, where imaginary $k_{\mathrm{B}}$ corresponds to attenuating modes, with the exponential decay $e^{-\left|k_{\mathrm{B}}\right|x_{1}}$; since $e^{in\pi}=\pm 1$ for integer $n$, we will also refer to $k_{\mathrm{B}}$ with Re$\,k_{\mathrm{B}}h=n\pi$ as imaginary, and interpret corresponding modes as non-propagating. We note that modes with $k_{\mathrm{B}}h=n\pi$ are at the boundary of the Brillouin zones, and represent standing waves. Complex $k_{\mathrm{B}}$ describes a progressive mode that exponentially decay according to Im$\,k_{\mathrm{B}}$. We emphasize that wave attenuation is not an indication of energy loss, as our system is non-dissipative; it is the result of a gradual scattering of energy to incoherent waves with zero mean. Bands of frequencies without real $k_{\mathrm{B}}$ roots are termed directional band gaps (or simply gaps) since there is a gap in the spectrum in the $x_{1}$ direction. When there are no real $k_{2}$ solutions across these bands, they are termed complete gaps, since all directions of propagation are prohibited.

2.3 An eigenvalue problem for $k_{2}$ using a plane wave approach

The standard plane wave expansion method to obtain the dispersion relation in elastodynamics dates back to Sigalas and Economou (1992) and Kushwaha et al. (1993). In photonics, the method has been extended by Hsue et al. (2005) to a formulation in which the Bloch wavenumber is the eigenvalue, and later on by Laude et al. (2009) for elastodynamics. A general analysis of the method, associated linear operators, and the properties of the eigenvalues was carried out recently by Mokhtari et al. (2019). Building upon the approach of Laude et al. (2009), we formulate for our settings an eigenvalue problem in which $k_{2}$ is the eigenvalue333We recall that $k_{2}$ is not a Bloch wavenumber and thus is not identified with any irreducible Brillouin zone (or, conversely, identified with an infinite one), as the laminate is homogeneous in the $x_{2}$ direction. . We begin by deriving the standard plane wave method by substituting into the Cauchy equations of motion

[TABLE]

the Bloch form

[TABLE]

where here and henceforth, the periodic part of $\left(\circ\right)$ is denoted by $\tilde{\left(\circ\right)}$. Since $\rho\left(x_{1}\right)$, $\mu\left(x_{1}\right)$, $\lambda\left(x_{1}\right)$ and $\tilde{u}_{i}\left(x_{1}\right)$ are periodic in $x_{1}$ with a period $h$, they can be written as

[TABLE]

where $\mathrm{G}=\frac{2\pi m}{h}$, $m\in\mathbb{Z}$, with the Fourier coefficients $\left\{\zeta\left(\mathrm{G}\right)\right\}$

[TABLE]

We have that

[TABLE]

Substituting Eqs.  (11) and (12) into Eq. (10) and factoring out $e^{i\left(\omega t-k_{\mathrm{B}}x_{1}-k_{2}x_{2}\right)}$ yield

[TABLE]

We multiply Eqs. (15) by $e^{i\mathrm{G^{\prime\prime}}x_{1}}$ and integrate over one period. Due to Fourier orthogonality, only terms satisfying $\mathrm{G^{\prime\prime}}=\mathrm{G}+\mathrm{G^{\prime}}$ remain, and we end up with an infinite set of equations that can be cast in matrix form as

[TABLE]

where $\mathsf{\tilde{u}}\left(\mathrm{G^{\prime}}\right)$ is a column vector comprising the Fourier coefficients of $\tilde{u}_{1}$ and $\tilde{u}_{2}$, and the matrices $\mathsf{A}_{0},\mathsf{A}_{1},\mathsf{A}_{2}$** **and $\mathsf{B}$ are given in Appendix B. Eq.  (16) constitutes a generalized eigenproblem for $\omega^{2}$ and $\mathsf{\tilde{u}}\left(\mathrm{G^{\prime}}\right)$ at prescribed $k_{2}$ and $k_{\mathrm{B}}$. A mixed generalized eigenproblem for $k_{2}$ at prescribed $\omega$ and $k_{\mathrm{B}}$ follows from Eq.  (16) using a state space-like formulation (Chapt.  10 in Deymier, 2013, see also Hussein et al., 2014), namely,

[TABLE]

For completeness, note that an equivalent formulation can be obtained using the displacement together with the stress and the eigenvector (Mokhtari et al., 2019). For computational purposes, the number of terms in the Fourier series is truncated, say by $-N\leq m\leq N$. The matrices $\mathsf{A}_{i},\mathsf{B},\mathsf{0}$ and $\mathsf{I}$ are accordingly of dimension $\left(2N+1\right)\times\left(2N+1\right)$.

2.4 Some properties of $k_{\mathrm{B}}$

and $k_{2}$ and their effect on energy transport

To highlight how the coupling between shear and pressure waves affect $k_{\mathrm{B}}$ and $k_{2}$, we record next their properties in the case of anti-plane shear in laminates, where such coupling is absent, and subsequently point out the differences. These differences have significant implications on frequency spectrum, and in turn the energy flow, since the slopes of propagating branches are a measure of its mean. More formally, the average energy flow is

[TABLE]

and satisfies (Willis, 2015)

[TABLE]

where $\left\langle\circ\right\rangle$ denotes averaging over one spatial and temporal period, and $E$ is the total mean energy density; $\bm{\mathcal{P}}$ is called the acoustic Poynting vector, and the $\bm{k}$-gradient of $\omega$ is identified as the group velocity.

Exceptional points.—In anti-plane shear, the size of the transfer matrix is $2\times 2$, whose two roots are $\lambda_{1,2}=e^{ik_{B1,B2}h}$. When the transfer matrix is real, one can show that if $\lambda$ is a solution, then so are $\lambda^{*}$ and $\lambda^{-1}$ . Since there are only two roots, either $\lambda=\lambda^{*}$ and then $k_{\mathrm{B}}$ is imaginary, or $\lambda^{*}=\lambda^{-1}$ and then $k_{\mathrm{B}}$ is real. Accordingly, eigenvalue degeneracies occur only at the edge of the irreducible Brillouin zones, where the modes are standing and not propagating. In this usual degeneracy, the eigenmodes remain linearly independent, and the splitting of the eigenvalues from the so-called diabolic points scales linearly.

The situation is significantly different when there are two coupled displacements as considered here. The size of the corresponding transfer matrix is $4\times 4$, thus has four roots. Hence, the fact that $\lambda^{*}$ and $\lambda^{-1}$ are also solutions does not enforce that $\lambda=\lambda^{*}$ or $\lambda^{*}=\lambda^{-1}$, and, in turn the exclusion of complex $k_{\mathrm{B}}$ roots. As we show in Sec. 3, not only such complex roots exist, they coalesce together with their eigenmodes at exceptional points inside the irreducible $1^{\mathrm{st}}$ Brillouin zone, according to a square-root scaling, constituting a Riemann surface structure.

Negative refraction.—We quantify the flow direction of each mode that is propagating in the laminate plane using the angle (Fig.  1)

[TABLE]

Srivastava (2016) showed analytically for anti-plane shear waves that $k_{2}$ and $\mathcal{P}_{2}$ share the same sign, and hence so does $\left\langle\mathcal{P}_{2}\right\rangle$. Therefore, the macroscopic transport of energy of anti-plane shear waves in the $x_{2}$ direction is aligned with the local transport. This implies that in the canonical configuration of excitation (Joseph and Craster, 2015),* i.e.*, when the laminate is impinged by a wave at a boundary normal to the lamination direction (Fig. 1b), anti-plane shear will always refract positively. The reason is that continuity requires the excited Bloch waves to share the same vertical wavenumber as the incident wave, and hence the incident angle and the angle of the Bloch waves have the same sign. In the present problem, where $\mathcal{P}_{2}$ is

[TABLE]

this is no longer the case; as we will demonstrate in the sequel, there are positive $k_{2}$ for which simultaneously $\left\langle\mathcal{P}_{2}\right\rangle<0$ and $\left\langle\mathcal{P}_{1}\right\rangle>0$, where

[TABLE]

This implies that negative refraction is realizable in the simple interface configuration.

Willis (2015) devised a complex configuration to achieve negative refraction of anti-plane shear waves, by considering waves impinging at an interface parallel to the lamination direction (Fig. 1c). We will additionally show that in-plane waves can refract negatively in this excitation setup as well.

3 Mode spectrum of an exemplary laminate

In this Sec., we study the frequency spectrum of in-plane waves propagating through an infinite laminate comprising steel and acrylic layers, whose properties are

[TABLE]

The same laminate was considered by Nemat-Nasser (2015), Willis (2015) and Srivastava (2016) to study anti-plane shear waves. The forthcoming study of the spectrum is the basis for the transmission analysis in Sec.  4. To simplify the presentation of the spectrum, we fix one of the three parameters $\left\{\omega,\ k_{\mathrm{B}},\ k_{2}\right\}$, and evaluate the relation between the remaining two.

\floatsetup

[figure]style=plain,subcapbesideposition=top

Our analysis begins by plotting in Fig. 2 the ordinary frequency $f=\omega/2\pi$ versus $0<k_{\mathrm{B}}h<\pi$ at prescribed $k_{2}$ values. Panels 2(a)-(c) correspond to $k_{2}h=0.5,\ 0.5i,$ and $0.5+0.5i$, respectively. Re$\,k_{\mathrm{B}}h$ and |Im$\,k_{\mathrm{B}}h$| (Im$\,k_{\mathrm{B}}h$ in panel c) are shown in solid and dashed curves, respectively, where the real and imaginary parts of a certain branch are plotted using the same color. Gaps are indicated by the red shading.

Panel 2(a) exhibits the following notable features. Firstly, there are no branches with real $k_{\mathrm{B}}$ from $f=0$ to $14.1\,\mathrm{kHz}$. In a study of the dependency of this gap on Re$\,k_{2}h$ (not shown here), we found that it widens as Re$\,k_{2}h$ increases. Across $136.2<f<181.2\ \mathrm{kHz}$, the branches are complex conjugates of each other; from $f=159.3\ \mathrm{kHz}$ the imaginary parts decrease until they vanish at $f\approx 181.2\ \mathrm{kHz}$. At this exceptional point the eigenvalues coalesce. Beyond this point there is a special narrow range where both branches are real, while their slopes have an opposite sign, i.e., the modes are propagating in opposite directions; this is the fingerprint of exceptional points.

Panel 2(b) exhibits different characteristics. Firstly, across the studied frequencies, one of the branches always has an imaginary part, hence the maximal number of modes that propagate in $x_{1}$ at any frequency is one. Contrary to the case in panel 2(a), here the first gap emerges above $f=0$, namely, at $f=180.6\ \mathrm{kHz}$. In a study of the dependency of this gap on Im$\,k_{2}h$ (not shown here), we found that gaps starting at $f=0$ emerge at higher values of Im$\,k_{2}h$. Thus, the diagram starts with a propagating band, where the horizontal group velocity changes sign at $f=148\ \mathrm{kHz}$. This flip of sign *inside *the irreducible first Brillouin zone is unique to in-plane waves.

Panel 2(c) shows an anomalous scenario without passbands, except at a discrete frequency ($f\approx 54.2\,$kHz) for which one branch is propagating. To facilitate the visual identification of this frequency, here we plot the signed imaginary part instead of its absolute value. Accordingly, this frequency is spotted by the zero crossing of Im$\,k_{\mathrm{B}}h$. There is no counterpart to this phenomenon in anti-plane waves, where complex $k_{2}$ are not accessible.

\floatsetup

[figure]style=plain,subcapbesideposition=top

Fig. 3 shows the (a) real and (b) imaginary parts of $k_{\mathrm{B}}h$ as functions of $f$ at $k_{2}h=0.5$, when we carry out an analytic continuation to the dispersion relation, such that the domain of the spectrum is formally extended to complex-valued frequencies. We emphasize that here, this is an artificial continuation (see, e.g., Shanin et al., 2018), as complex frequencies are not accessible (we are concerned with time-harmonic waves in a system with no dissipation). Lu and Srivastava (2018), for example, considered an analytic continuation of the shear modulus, which has the physical interpretation of a system with gain or loss.

The spectrum exhibits a structure of a Riemann surface in the vicinity of the exceptional point: Fig. 2(a) is thus a section of that surface, at the plane Im$\,f=0$. In view of the duality of between the frequency and the wavenumbers (Torrent et al., 2018) in the different forms of the eigenvalue problem (Mokhtari et al., 2019), a Riemann surface is expected when $f$ is fixed and the roots of $k_{\mathrm{B}}$ are evaluated against $k_{2}$, as we will demonstrate in the sequel. Importantly, the states of the system in the complex wave vector space are accessible, and are not merely a formal extension.

\floatsetup

[figure]style=plain,subcapbesideposition=top

Towards this end, we fix $f$ and examine the spectrum in the complex $k_{\mathrm{B}}-k_{2}$ space. Figs. 4(a) and 4(b) show branches of purely real $k_{\mathrm{B}}$ versus Re$\,k_{2}h$ and Im$\,k_{2}h$, at $f=100\ \mathrm{kHz}$ and $160\ \mathrm{kHz}$, respectively444The diagrams were evaluated using the extended plane wave method with 51 plane waves in the expansion. A comparison with the exact hybrid matrix method can be found in Appendix C.. Segments of real, imaginary and complex $k_{2}$ are denoted by blue, red, and grey, respectively. Notably, the number of modes with real or imaginary $k_{2}$ is finite. In our study (not shown here) the number of modes with real or imaginary $k_{2}$ increases with the frequency. We further note that in our study on branches with complex $k_{2}$, we found that their number is infinite, although it cannot be observed from the truncated diagram we show. By contrast, Srivastava (2016) showed that in the anti-plane motion there are no complex $k_{2}$, and the number of branches with imaginary $k_{2}$ is infinite. In both cases (anti-plane and in-plane waves), the existence of an infinite number of decaying modes conforms with the need of such a set in satisfying the continuity of the displacement and traction across certain interfaces (Srivastava and Willis, 2017). \floatsetup[figure]style=plain,subcapbesideposition=top

Modes that propagate along both $x_{1}$ and $x_{2}$ are studied in Fig. 5, by restricting attention to the plane Im$\,k_{2}h=0$, and evaluating $k_{\mathrm{B}}h$ versus Re$\,k_{2}h$ at prescribed frequencies. The real and imaginary parts of $k_{\mathrm{B}}$ are shown in solid and dashed curves, respectively, where each color corresponds to a different branch.

At $f=60\ \mathrm{kHz}$ (panel 5a), two modes of different branches are propagating when $0<k_{\mathrm{B}}h<0.919$; only one mode propagates when $0.919<k_{\mathrm{B}}h<3.101$, and at larger Bloch wavenumbers there are no propagating modes. Here and at higher frequencies, we find that beyond the illustrated $k_{2}$-interval there are no propagating branches—solutions at higher $k_{2}$ values are always given by $k_{\mathrm{B}}$ with an imaginary part, such that $k_{\mathrm{B}}$ is either imaginary or complex.

At a slightly higher frequency ($f=60.1\ \mathrm{kHz}$, panel 5b), the whole $1^{\mathrm{st}}$ Brillouin zone admits real solutions, contrary to the case at $f=60\ \mathrm{kHz}$. This change in the spectrum between the frequencies has a significant implication on beam steering, as we will show in Sec. 4.2. In a narrow range near $k_{2}h=0$, the blue mode has a non-zero imaginary part and its real part equals $\pi$. Thus, for the frequencies in panels 5(a)-(b) and displayed range of $k_{2}h$ values, the Bloch wavenumber of the modes is either real, or has an imaginary part with Re$\,k_{\mathrm{B}}h=0$ mod $\pi$.

Panel 5(c), which evaluates the modes at $f=100\ \mathrm{kHz}$, also shares the latter characteristic. Here, however, the real part of the blue is a non-monotonic function of $k_{2}h$, as it changes trend twice. This is the fingerprint of an exceptional point, and an indicator of negative refraction as will be demonstrated in Sec. 4.2. Similarly to the observation made regarding the modes near the exceptional point in Fig. 2, the slopes of the two modes have different sign. There are accordingly three propagating modes across the range $1.87\apprle k_{\mathrm{B}}h\lesssim 2.03$ associated with the same branch, which exhibit different vertical wavelengths. These features are unique to the in-plane motion, as in the anti-plane motion the real part of $k_{\mathrm{B}}h$ is either monotonically increasing to $\pi$ or monotonically decreasing to [math] (Srivastava, 2016).

At $f=160\ \mathrm{kHz}$ (panel 5d) across the interval $0<\mathrm{Re}\,k_{2}h\lesssim 0.79$, there are two complex conjugate pairs of $k_{\mathrm{B}}$ solutions with and identical real part, up to $k_{2}h\approx 0.79$. At this exceptional point, the imaginary part of all the branches vanishes, and the modes coalesce. Beyond this point, the real part of the branches diverges in a different directions, similarly to the feature observed in Fig. 2(a). At $k_{2}h\approx 0.8$, the blue branch becomes attenuating, and at $k_{2}h\approx 1$ the green branch also becomes imaginary. The blue mode is propagating again for $3.08<k_{2}h<3.22$; beyond that range there no more propagating modes.

When $f=181.3\,\mathrm{kHz}$ (panel 5e), there are two exceptional points at $k_{2}h=0.503$ and $0.93$. Again, in the vicinity of these points there are modes with slopes of different sign. As in panel 2(d), there is a range of Bloch wavenumbers with three propagating modes when $k_{\mathrm{B}}h>1.26$, where for lower $k_{\mathrm{B}}h$ there are two propagating modes.

We recall that in Fig.  2(a), we exhibited an exceptional point for $f=181.2\,\mathrm{kHz}$ and $k_{2}h=0.5$, i.e., when perturbing one of the parameters of the system about this point (e.g., $k_{2}h$ is perturbed from 0.5 to 0.503), another exceptional point emerges at a variation of another parameter of the system (e.g., is perturbed $f$ from 181.2 to $181.3\,\mathrm{kHz}$). This is not accidental; in fact, the exceptional points we demonstrate for each set of parameters comprise exceptional curves in a higher dimensional space.

At $f=340\,\mathrm{kHz}$ (panel 5f), the propagating green mode is a non-monotonic function of $k_{2}$, where initially it has a positive slope up to a maximal point, beyond which the slope changes sign, until the mode becomes imaginary at $k_{2}h=1.72$. The blue mode is either real or imaginary with $\mathrm{Re}\,k_{\mathrm{B}}h$ mod $\pi=0$. The number of propagating modes for prescribed $k_{\mathrm{B}}$ is higher than at lower frequencies, namely, a minimal number of three and a maximal number of five across $1.624<k_{\mathrm{B}}h<1.88$.\floatsetup[figure]style=plain,subcapbesideposition=top

As mentioned earlier, we expect the spectrum to exhibit the structure of a Riemann surface near exceptional points in the complex $k_{\mathrm{B}}-k_{2}$ space. This is shown in Fig. 6, where the (a) real and (b) imaginary parts of $k_{\mathrm{B}}h$ are plotted against Re$\,k_{2}h$ and Im$\,k_{2}h$, at $f=160\,$kHz: Fig. 5(d) is thus the section Im$\,k_{2}=0$ of this surface. We emphasize again the the states of the system over this manifold are accessible since they correspond to realizable wavenumbers, and are not merely a formal extension.

\floatsetup

[figure]style=plain,subcapbesideposition=top

Finally, in Fig. 7 we fix $k_{\mathrm{B}}h$ and evaluate the spectrum in the $f-k_{2}$ space. Specifically, panels  7(a)-(c) correspond to $k_{\mathrm{B}}h=0.5,\,0.5i$ and $0.5+0.5i$, respectively; segments of real, imaginary and complex $k_{2}$ are denoted by blue, red, and grey, respectively. We observe that the number of modes in general—and propagating in particular—increases with frequency.

In panel 7(a) we observe an exceptional point when $k_{2}h=0.9i$ at 13 kHz, where a single complex mode branches to two imaginary modes; as frequency increases to 33 kHz, one of these branches becomes a propagating mode555Laude et al. (2009) discussed how the shifting of modes between the complex and real domains serves as a mechanism to preserve the total number of modes at a prescribed frequency. . A similar exceptional point is identified in panel 7(b), when $k_{2}h=2.3i$ at $72\,$kHz. There are no propagating modes in panel 7(c).

4 Transmission across an interface of a semi-infinite laminate

As mentioned, modes with either prescribed $k_{\mathrm{B}}$ or prescribed $k_{2}$ are excitable through two different configurations. We analyze the transmission and excited modes in these configurations, linking them to the exceptional points, negative refraction, and beam steering found in the spectrum, as reported in Sec.  3.

4.1 Metamaterial refraction via an interface parallel to the layers

To excite modes with designated vertical lengths, we connect the laminate at $x_{1}=0$ to a homogeneous half-space whose properties are denoted by the script 0, from which an incident wave propagates in an angle $\theta_{i}$ towards the interface (Fig. 1b). The wave is described by the potential $V$

[TABLE]

where the wavenumber $k=\omega/c^{\left(0\right)}$ is either $\omega/c^{\left(0\right)}_{S}$ for a shear wave, or $\omega/c_{L}^{\left(0\right)}$ for a pressure wave. The incident displacements are derived according to Eq. (2), by setting $\psi^{\left(0\right)}=V$ for incident shear, or $\phi^{\left(0\right)}=V$ for incident pressure wave. The incident wave is partially reflected back by the interface, and partially transmitted to the laminate as Bloch modes. The continuity conditions at $x_{1}=0$ enforces all these waves to have an $x_{2}$-dependency in the form $e^{-ik\cos\theta_{i}x_{2}}$. Corresponding Bloch modes are extracted by the procedure in Sec. 2.1 with $k_{2}=k\cos\theta_{i}$, in the form

[TABLE]

where the reflected waves are derived from

[TABLE]

such that $k_{L}^{\left(0\right)}$ and $k_{S}^{\left(0\right)}$ are related to $k_{2}=k\cos\theta_{i}$ via Eq. (5). Finally, the reflection coefficients $R_{l}^{\phi},R_{l}^{\psi}$ and transmission coefficients $T_{m}$ are determined from the continuity of the state vector at $x_{1}=0$.

We numerically demonstrate how negative refraction is obtained, using a half-space with the properties

[TABLE]

We set an incident pressure wave at $\theta_{i}=77.22^{\circ}$ and $f=100\,\mathrm{kHz}$, for which $k_{2}h=1$. The resultant Bloch wavenumbers are extractable from Fig. 5, where at the pertinent $k_{2}$ there is one attenuating mode and one propagating mode. We recall that the slope in Fig. 5 of the propagating mode at this $k_{2}$ is positive, while at higher $k_{2}$ the slope is negative, which is an indicator of negative refraction. This is verified by evaluating Eqs. (21)-(22), to find that the propagating angle is $\theta=107.03{}^{\circ}$. As a consistency check, we examine the balance of energy at $x_{1}=0$, namely,

[TABLE]

where the scripts $T_{m},\ \phi,\ \psi$ and $I$ denote quantities associated with the corresponding fields. Indeed, we find that the terms in left hand side are $0.617,\ 0.316$ and $0.067$, respectively, which sum to $1$. Note that only the propagating mode contributes to the first term, since the energy flux of the mode with the imaginary $k_{\mathrm{B}}$ vanishes.

As a second example, we set an incident shear incident wave at $\theta_{i}=89.62{}^{\circ}$ and $f=340\,\mathrm{kHz}$, for which $k_{2}h=0.25$. The resultant Bloch wavenumbers are extractable from Fig.  5(f). Here again, at the prescribed $k_{2}$ there is one attenuating mode, and one propagating mode with a positive slope, that will change sign at greater $k_{2}$. Calculation of Eqs. (21)-(22) confirms a negative refraction at $\theta=106.2^{\circ}$. Energy balance is also verified, where the calculation of the terms in the left hand side of Eq. (28) provides $0.148,\ 0.007$ and $0.845$, respectively, which again sum to $1$.

Fig.  8 shows frequency domain finite element simulations using COMSOL Multiphysics*®* the energy flux along $x_{2}$ ($\mathcal{P}_{2}$) resulting from an (a) incident pressure wave at $\theta_{i}=77.22^{\circ}$ and $f=100\,\mathrm{kHz}$, and (b) an incident shear incident wave at $\theta_{i}=89.62{}^{\circ}$ and $f=340\,\mathrm{kHz}$. We used low reflecting boundary conditions to avoid reflections due to the finite truncation of the spatial domain. The incident wave was excited using load control over a defined line. We used triangular elements, whose maximal size was third of the unit-cell thickness. Indeed, the resultant negative angles of refraction agree with the analytical prediction.

\floatsetup

[figure]style=plain,subcapbesideposition=top

4.2 Negative refraction via an interface normal to the layers

It is possible to excite designated Bloch wavenumbers using an incident wave from a homogeneous half-space that is bonded to the laminate by an interface along the lamination direction, say at $x_{2}=0$ (Fig. 1c). This configuration was proposed by Willis (2013) to achieve negative refraction of anti-plane shear waves, and was studied in greater detail by Willis (2015), Nemat-Nasser (2015), Srivastava (2016), and Morini et al. (2019). We recall that in this setup, an incident wave will excite an infinite number of decaying waves in the laminate, required to enforce continuity conditions (Srivastava and Willis, 2017). Here, we firstly develop a method to resolve the normal mode decomposition and energy partition of the excited in-plane waves in this configuration.

Consider again an incident wave according to the potential (24), now impinging on the horizontal interface $x_{2}=0$. In this setting, the continuity conditions at the interface enforce the reflected and transmitted waves to share the same horizontal length with the incident wave. Accordingly, the transmitted waves have the form

[TABLE]

where $k_{\mathrm{B}}=k\sin\theta_{i}$, and $\left\{k_{m}\right\}$ are extractable from Eq. (17). In view of the foregoing observations, the reflected waves is constructible from potentials that must have the form

[TABLE]

where

[TABLE]

and $U_{l}\left(x_{1}\right)=e^{-2l\pi ix_{1}/h}$.

The transmission and reflection coefficients $\left\{T_{m},R_{l}^{\phi},R_{l}^{\psi}\right\}$ are determined by the continuity conditions at $x_{2}=0$, which are compactly written as ****

[TABLE]

using

[TABLE]

For computational purposes, we truncate the sums over $l$ and $m$ such that $-N_{R}\leq l\leq N_{R}$ and $1\leq m\leq 2\left(2N_{R}+1\right)$, and apply the elegant orthogonality relation666Alternatively, we can use the Fourier orthogonality $\int_{0}^{h}U_{\alpha}\left(x_{1}\right)U_{\beta}^{*}\left(x_{1}\right)dx_{1}=h\delta_{\alpha\beta}$. (Mokhtari et al., 2019)

[TABLE]

to obtain an algebraic system of equations for $\left\{T_{m},R_{l}^{\phi},R_{l}^{\psi}\right\}$ in the form

[TABLE]

with diagonal submatrices whose elements are given in Appendix D; thus, the size of $\mathsf{M}_{TT}^{\mathsf{s}}$ is $2\left(2N_{R}+1\right)\times 2\left(2N_{R}+1\right)$, of $\mathsf{M}_{T\phi}^{\mathsf{s}}$ and $\mathsf{M}_{T\psi}^{\mathsf{s}}$ is $2\left(2N_{R}+1\right)\times\left(2N_{R}+1\right)$, of $\mathsf{M}_{\phi T}^{\mathsf{s}}$ and $\mathsf{M}_{\psi T}^{\mathsf{s}}$ is $\left(2N_{R}+1\right)\times 2\left(2N_{R}+1\right)$, and the size of the remaining submatrices is $\left(2N_{R}+1\right)\times\left(2N_{R}+1\right)$. In summary, Eq. (35) constitutes $4\left(2N_{R}+1\right)$ algebraic equations for $4\left(2N_{R}+1\right)$ scattering coefficients. The resultant coefficients should satisfy the following energy balance

[TABLE]

An example of the fields at the interface as obtained from the method is provided in Appendix E.

We proceed to analyze the transmission of Bloch modes in laminate (23) that are induced by impinging pressure wave from the half-space

[TABLE]

To this end, we plot in Fig. 9 the propagation angle $\theta$ of the modes analyzed in Fig. 5, as functions of the incoming angle $\theta_{i}$ and its corresponding $k_{\mathrm{B}}h$, indicated by the upper and lower axes, respectively. The calculation of $\theta$ was carried out through direct evaluation of Eqs. (21)-(22) for calculate $\arctan\frac{\left\langle\mathcal{P}_{1}\right\rangle}{\left\langle\mathcal{P}_{2}\right\rangle}$.

Panel 9(a) is for $f=60\,$kHz, which corresponds to Fig. 5(a), and uses the same color legend. Our first observation concerns the change of sign in $\theta$ with respect to $k_{\mathrm{B}}h=\pi$, owing to the $2\pi$-periodicity of the spectrum in Re$\,k_{\mathrm{B}}h$, and its and reflection symmetry within that Brillouin zone period. Accordingly, $\partial\omega/\partial k_{\mathrm{B}}$ changes sign, and hence so does $\left\langle\mathcal{P}_{1}\right\rangle$, as noted by Willis (2015) in his analysis of anti-plane shear. The diagram exhibits 2-fold rotational symmetry about $\left(\pi,0\right)$, which corresponds to $\theta_{i}\approx 30^{\circ}$. Notably, the transmission angle of the blue mode increases very fast near $k_{\mathrm{B}}h=3.101$, and then discontinuously flips from $90^{\circ}$ to $-90^{\circ}$ at $k_{\mathrm{B}}h=3.182$. This discontinuity occurs since the blue branch does not exist (a gap) when $3.101<k_{\mathrm{B}}h<\pi$, as highlighted in Fig. 5(a). Therefore, in the vicinity $\theta_{i}\approx 30^{\circ}$, it is sufficient to slightly change the incident angle in order to significantly steer the transmitted mechanical beam. Wide steering by small changes of the incident angle, albeit less extreme, occurs also about $\theta_{i}\approx 10^{\circ}$ and $55^{\circ}$. We further note that since there is a range of $k_{\mathrm{B}}h$—and hence of incident angles—in which there is only one propagating mode (see Sec. 3), it is possible to achieve purely negative transmission. Analogous observations were made by Srivastava (2016) for anti-plane shear waves.

Panel 9(b) is for $f=60.1\,$kHz, which corresponds to Fig. 5(b). We observe that the discontinuous sign flip of $\theta$ and significant beam steering are lost, as the $\theta\left(k_{\mathrm{B}}h\right)$-curve smoothly passes through the point $\left(0,\pi\right)$. This result agrees with the extension of the blue branch in Fig. 5(b) to $\pi$ near $k_{2}h=0$, where it has an imaginary part, hence becomes attenuating. Thus, in the vicinity of $f=60\,\mathrm{kHz}$, energy transport is also very sensitive to the excitation frequency.

Fig. 10 is for $f=181.3\,$kHz, corresponding to Fig. 5(e), which we recall exhibits two exceptional points; a zoom in about these points is shown in the right side of the figure. The color legend is different here: purple denotes solutions with $k_{2}h<0.5$, i.e., left to the first exceptional point in Fig. 5(e), where grey denotes solutions right to the second exceptional point, with $0.93<k_{2}h<1.21$. (Note that we do not include the third propagating branch with $3.4<k_{2}h$.) We mark the segments above and below the exceptional points Fig. 5(e) using circle marks and diamond marks, respectively.

Up to the first exceptional point $\left(k_{\mathrm{B}}h\approx 2.564,\theta_{i}\approx 7.753^{\circ}\right)$, the green branch of Fig. 5(e) admits two propagation angles which—remarkably—one of them is negative inside the first Brillouin zone. Beyond the exceptional point, mode switching occurs between the green branch (purple diamonds) and the blue branch (purple circles) of Fig. 5(e), as the former becomes attenuating and the latter becomes propagating with positive refraction. Similar switching occurs between the green and blue branches beyond the second exceptional point $\left(k_{\mathrm{B}}h\approx 3.021,\theta_{i}\approx 9.146^{\circ}\right)$. These phenomena are unique to the in-plane motion, owing to exceptional points, and are not accessible in the anti-plane setting. Finally, we note that here again, the diagram exhibits 2-fold rotational symmetry about $\left(\pi,0\right)$, which corresponds to $\theta_{i}\approx 9.514^{\circ}$, followed by the inversion of the aforementioned phenomena with respect to the incident angle.s

5 Conclusions

We have revisited the problem of in-plane waves propagation in periodic laminates, were the motivation was twofold. Firstly, this study is a necessary complement to the reports of Willis (2015), Nemat-Nasser (2015) and Srivastava (2016) on metamaterial phenomena in the model problem of anti-plane waves traversing elastic laminates. Secondly, this problem allows us to show that the coupling between shear and pressure parts can be harnessed for anomalous energy transport, within a relatively simple analytical study.

We have shown that the corresponding spectrum contains exceptional points at which two Bloch modes coalesce, and further showed that these points are the source of anomalous energy transport. To this end, we have determined how energy is scattered when an incident in-plane wave impinges the laminate in two model interface problems. We found metamaterial transmission through the laminate, such as pure negative refraction, and beam splitting and steering, at states of the system near the exceptional points. Notably, we achieved negative refraction in the canonical transmission problem, where the laminate is impinged by an incoming wave from a homogeneous medium whose interface with the laminate is parallel to the layers.

We emphasize again that these phenomena emerge from the unique coupling in elastodynamics between the volumetric and distortional modes of deformation, which cannot be observed in sound and light waves. Our work further paves the way for encircling exceptional points in a tangible, purely elastic apparatus, for applications such as asymmetric mode switches.

Acknowledgments

We are grateful to Profs. Nimrod Moiseyev, Alexei Mailybaev and Ankit Srivastava for fruitful discussions. We thank Dr.  Pernas-Salomón for his help regarding the hybrid matrix method. We acknowledge the support of the Israel Science Foundation, funded by the Israel Academy of Sciences and Humanities (Grant no. 1912/15), and the United States-Israel Binational Science Foundation (Grant no. 2014358), and Ministry of Science and Technology (Grant no. 880011).

Appendix A

The forthcoming formulation is an adaptation of the developments of Shmuel and Pernas-Salomón (2016) and Pernas-Salomón and Shmuel (2018) to the current problem. Here, the modified state vector—consisting of quantities which are continuous across the interface between two adjacent layers—is

[TABLE]

where the normalizing coefficients required to stabilize subsequent calculations are

[TABLE]

and

[TABLE]

The state vector can be written as follows

[TABLE]

where

[TABLE]

$\mathsf{Q}\left(x_{1},x_{2}\right)$ can also be defined as

[TABLE]

Using this definition the hybrid matrix of the $\emph{n}^{\mathrm{th}}$ layer is given by

[TABLE]

Using the definition in Eq. (A.9), we relate the field variables which appear in the modified state vector at the ends $x_{0}$ and $x_{0}+h^{\left(\emph{p}\right)}$ of the $\emph{n}^{\mathrm{th}}$ layer via

[TABLE]

where $\mathsf{H}^{\left(\emph{n}\right)}\left(h^{\left(\emph{p}\right)}\right)$ is the $4\times 4$ hybrid matrix of the $\emph{n}^{\mathrm{th}}$ layer. In order to calculate the total hybrid matrix of n layers, denoted by $\mathsf{H}^{\left(1,n\right)}$, a recursive algorithm is used in terms of the hybrid matrix of the first $\emph{n}-1$ layers ($\mathsf{H}^{\left(1,n-1\right)}$) and the hybrid matrix of the $n^{\mathrm{th}}$ layer $\mathsf{H}^{\left(\emph{n}\right)}$, namely,

[TABLE]

Where $\mathsf{H}_{11}^{\left(1,k\right)},\mathsf{H}_{12}^{\left(1,k\right)},\mathsf{H}_{21}^{\left(1,k\right)},\mathsf{H}_{22}^{\left(1,k\right)}$ denote the $2\times 2$ blocks of the hybrid matrix $\mathsf{H}^{\left(1,k\right)}$ of the first $k$ layers. The generalized eigenproblem presented in Eq. (8) yields the characteristic equation

[TABLE]

for $\lambda=e^{ik_{\mathrm{B}}h}$, where

[TABLE]

and $h_{ij}$ are the components of the total hybrid matrix. The latter two equalities imply that $\lambda^{-1}=e^{-ik_{\mathrm{B}}h}$ is also a solution, as expected. Eq. (A.12) yields a closed-form expression for $\cos k_{\mathrm{B}}h$, given in Eq. (9).

Appendix B

The matrices $\mathsf{A}_{0},\mathsf{A}_{1},\mathsf{A}_{2}$** **and $\mathsf{B}$ appearing in Eq. (16) are given by

[TABLE]

where the $\left(\mathrm{G},\mathrm{G^{\prime}}\right)$ component of each block is

[TABLE]

Appendix C

This appendix demonstrate the applicability of the extended plane wave expansion method, by comparing its results with the exact results of the hybrid matrix method. Fig. 11 shows real $k_{\mathrm{B}}h$ as function of the real and imaginary parts of $k_{2}h$ for laminate (23) and $f=100\ \mathrm{kHz}$. Points calculated with the extended plane wave expansion method and hybrid matrix method are depicted in purple and blue, respectively. We used 51 plane waves for the extended plane wave expansion method and points calculated using the hybrid matrix method differ from the results of the hybrid matrix by $5\times 10^{-2}$.

Appendix D

The components of the submatrices in Eq. 35 are

[TABLE]

where subscript numbers denote the component of $\mathsf{s}^{\perp}$, and subscript letters denote the potential from which this component is derived.

Appendix E

\floatsetup

[figure]style=plain,subcapbesideposition=top

Fig. 12 shows $\mathsf{s}_{1}^{\perp}$ (panel a), $\mathsf{s}_{2}^{\perp}$ (panel b), $\mathsf{s}_{3}^{\perp}$ (panel c, normalized by $k\bar{\mu}$), and $\mathsf{s}_{4}^{\perp}$ (panel d, normalized by $k\bar{\lambda}$) at the interface $x_{2}=0$ as calculated via Eq. (35). Specifically, solid and dashed curves correspond to the laminate and homogeneous half-space, respectively, where the real and imaginary parts of each field are denoted by blue and purple, respectively. The laminate and homogeneous half-space that was used have properties (23) and (27), respectively, and the calculation was carried out for $f=100\,\mathrm{kHz}$, $\theta_{i}=15.72^{\circ}$, and an incident shear wave. For this example, we chose $N_{R}=12$, which ensures all the propagating transmitted and reflected modes in the $x_{2}$ direction are included; the rest of the transmitted modes are those to have the lowest absolute value of $\mathrm{Im}k_{2}h$. The solution obtained does not match the fields perfectly, however it is sufficient on average provides a good approximation. The solution yields the following terms for the energy balance (36),

[TABLE]

which sum to 0.988: a difference of only 0.12%. In this case, the transmitted mode refracts positively with transmitted angle of $\theta=1.12^{\circ}$, however carries only a small fraction of the energy of the incident wave.