Introduction

Metamaterials and quantum graphs

Since Veselago published his pioneering paper in 1968—The Electrodynamics of Substances with Simultaneously Negative Values of Permeability and Permittivity1—it has been understood that the manipulation of electromagnetic material properties can give rise to exotic wave effects. In this paper, Veselago showed that when a material has both permittivity and permeability less than zero, its refractive index would be negative. By balancing the wave vector and the Poynting vector, it was shown that the electromagnetic waves within would have anti-parallel phase and group velocities. As a result, Snell’s law and the Doppler and Vavilov–Cherenkov effects are reversed. These results were considered purely theoretical, since no such material was known to exist.

In 2000, Pendry demonstrated a remarkable application of negative refractivity. Showing, given an ideal lossless slab of material with a refractive index of \(n=-1\) in a medium with equal and opposite refractive index, a lens can be theoretically constructed that could focus light perfectly2. Crucially such a material could overcome the diffraction limit of a traditional lens. Pendry presents a practical way one could engineer such a material. Proposing the periodic arrangement of unit cells made from C-shaped metal elements or “split-ring-resonators” with wires, giving rise to an effective negative refractive index within some frequency domain. This was experimentally demonstrated by Smith3. Since then, these man-made materials or “metamaterials” have been investigated extensively. There have been countless proposals for different resonator designs and arrangements, giving rise to a large number of different wave effects.

Metamaterials function due to the interplay between the wavelength and the scale of the unit cell. For wavelengths of the order or less than that of the unit cell, the waves undergo Bragg scattering interacting directly with each resonant element. However, in the long wavelength regime, the material appears continuous with properties owing to the underlying structure4,5. By varying these structural or resonant elements, one can achieve the desired wave effects. The required constituents of a metamaterial continue to be a matter of debate and various modelling techniques are used such as, transmission line models6, boundary element models7 and Finite Element Analysis8 to gain insight into the design changes required to achieve the desired effects. These simulation techniques can be quite time consuming when setting up the models, let alone considering design modifications, so there is a need for simple and computationally cheap models for proof of principle studies. The proposed quantum graph approach can provide this simplified tool for designing metamaterial with a large range of different properties.

Since their introduction in 19979, quantum graphs have become a most valuable model for studying quantum and wave effects in the context of quantum chaos and beyond. The theory of quantum graphs describe a network constructed from a set of vertices that act as point scatterers connected by a set of one-dimensional edges (bonds), wherein waves are free to travel. The wave transport through the network depends on the chosen properties of the point scatterers characterised by unitary scattering matrices. By changing the elements of the scattering matrix, one can customise the wave transport through the network. The versatility and ease of construction of quantum graphs, as well as the finite dimensionality and the exactness of ‘semiclassical’ expressions, make them ideal toy models to test ideas and research hypothesis’—see10 and11 for an overview. Besides quantum chaos, applications encompass modelling the vibrations of coupled plates12, formulating quantum random walks13,14 as well as quantum search algorithms15. One advantage of a quantum graph formulation is, that eigenvalue conditions can be written in terms of a secular equation involving the determinant of a matrix of finite dimension. Similarly, the scattering matrix, describing the properties of open quantum graphs, can be given in terms of a closed form expression involving finite dimensional matrices.

Naturally the concise language of quantum graph theory lends itself well to describing metamaterials. The lattice nodes act as scatterers representing sub-wavelength resonant elements and the edges model phase modulation between each resonant element. For this work, the metamaterial is modeled as an infinite square periodic quantum graph embedded in a 2D real space. A generalisation to 3D is straightforward by adding additional edges and vertices in the extra dimension. With this simple quantum graph construction, a multitude of nontrivial wave effects can be modelled by adjusting the parameters of the graph metric and the node scattering matrix. By manipulating the resulting band diagrams, we will demonstrate some of these effects, namely positive and negative refraction as well as beam steering.

The paper is structured as follows: in section “Waves in finite open quantum graphs—modelling resonant element” we will describe the set up of open quantum graphs which will serve as the resonant elements in the metamaterials. We will then construct period quantum graphs in section “A quantum graph model for a metamaterial” and derive the equations for constructing plane wave solutions and the dispersion curves in section “Wave propagation in metamaterials—plane wave solutions”. In section “Wave propagation in metamaterials—Gaussian beams”, we consider wave energy transport based on a Gaussian beam representation of the wave solutions and focus then on the available parameter space used for band engineering in section “Band engineering”. In the last two sections, section “Wave refraction between metamaterials” and section “Waves in N layered metamaterials”, we focus on reflection/transmission at “interfaces” between different metamaterials both in the two and N-layer case.

Waves in finite open quantum graphs—modelling resonant element

Metamaterials are typically constructed from a periodic arrangement of sub-wavelength resonant elements. We consider here in particular two-dimensional metamaterials, a generalisation to three dimensions is straight-forward. Using quantum graph theory, we can model a single resonant element as an open quantum graph, which is then periodically arranged to construct the metamaterial.

Figure 1
figure 1

Three examples of open quantum graphs are shown. (a) Represents a cross resonator, (b) a ring resonator and (c) some zig zag structured resonator. The edges between vertices are depicted in black and open leads in grey.

Full size image

Consider a connected and open quantum graph \(\Gamma (V,E,L)\) as a model for a resonator. The graph is constructed from a set of bidirectional edges \(E = \{e_{0},e_{1},\ldots ,e_{|E|-1}\}\) with a corresponding metric \({\mathcal {L}} = \{l_{e_{0}},l_{e_{1}},\ldots ,l_{e_{|E|-1}}\}\), connected by a set of vertices \(V = \{V_{0},V_{1},\ldots ,V_{|V|-1}\}\) with corresponding boundary conditions. In addition, a set of four leads \(L = \{L_{l}, L_{r},L_{d},L_{u}\}\) (semi-infinite edges) are imposed in the left(l), right(r), down(d) and up(u) directions serving as incoming and outgoing channels in the two dimensional lattice, see the examples in Fig. 1. Both the edges in E and leads in L are endowed with a one-dimensional wave equation,

$$\begin{aligned} \left( \frac{\partial ^{2}}{\partial z^{2}} + k^{2} \right) \varvec{\psi }(z) = 0. \end{aligned}$$
(1)

Here \(\varvec{\psi }(z) = (\varvec{\psi }_{L}(z),\varvec{\psi }_{E}(z))^{T}\) is the vector of all lead and edge solutions and k is the wave number. The solution on a given edge \(\psi _{e_{j}}(z)\) exists for \(z \in [0,l_{e_{j}}]\), while the solution on a given lead \(\psi _{L_{j}}(z)\) exists for \(z \in [0,\infty )\). (For ease of notation, z represents the spacial coordinates on any edge or lead). The eigenfunction solution \(\varvec{\psi }(z)\) is expressed as a superposition of counter propagating plane waves,

$$\begin{aligned} {\varvec{\psi }}(z) = {\varvec{a}}^{+} e^{ikz} + {\varvec{a}}^{-} e^{-ikz}. \end{aligned}$$
(2)

Here \({\varvec{a}}^{\pm } = ({\varvec{a}}_{L}^{\pm },{\varvec{a}}_{E}^{\pm })^{T}\) represents the vector of all wave amplitudes heading in (−) or out (\(+\)) of a vertex between the leads and edges. The vertices in V are generally considered to be point scatterers, and the mapping from all incoming to all outgoing wave amplitudes at the vertices is described by a unitary matrix \({\hat{S}}\) , that is,

$$\begin{aligned} {\varvec{a}}^{+} = {\hat{S}}{\varvec{a}}^{-}. \end{aligned}$$
(3)

We assume Neumann boundary conditions at individual vertices enforcing wave continuity and flux conservation at each vertex10 although more general vertex scattering matrices can be considered16). Enforcing Neumann boundary conditions, we obtain for the \(pq\hbox {th}\) element of the scattering matrix associated with vertex \(V_{j}\)

$$\begin{aligned} \left\{ {\hat{S}}_{V_{j}}\right\} _{pq} = \frac{2}{v_{_{V_j}}} - \delta _{pq}. \end{aligned}$$
(4)

Here, \(\delta _{pq}\) is the Kronecker delta, and \(v_{_{V_j}}\) represents the valency (number of edges and leads attached to a given vertex) at vertex \(V_{j}\).

To treat the connected graph as a model for a resonator, we wish to understand how waves incoming from any of the leads are redistributed and modulated into outgoing waves leaving through the leads. To do this a four dimensional graph scattering matrix \({\hat{S}}_{\Gamma }\) is constructed,

$$\begin{aligned} {\varvec{a}}_{L}^{+} ={\hat{S}}_{\Gamma }(k) {\varvec{a}}_{L}^{-}, \end{aligned}$$
(5)

reducing the entire connected graph in \(\Gamma (V,E,L)\) to a single scatterer illustrated in Fig. 2.

Figure 2
figure 2

(a) Some arbitrary compact quantum graph \(\Gamma (V,E,L)\). (b) The graph having been reduced to a single vertex. Noted are the lead wave amplitudes \({\varvec{a}}_{L}^{\pm }\).

Full size image

To construct \({\hat{S}}_{\Gamma }\) from the underlying connected quantum graph, one decomposes the scattering matrix into discrete events. First, one accounts for the prompt reflections from the leads back to the leads in a matrix \({\hat{S}}_{LL}\). Second, the waves on the leads are coupled to the edges in E with a matrix \({\hat{S}}_{EL}\). Third, the edge dynamics are described by an infinite series of scattered paths expressed as a Neumann series, \(\sum _{n=0}^\infty ({\hat{P}}(k) {\hat{S}}_{EE})^n {\hat{P}}(k) := [ {\hat{\mathbb{I}}} - {\hat{P}}(k){\hat{S}}_{EE}]^{-1} {\hat{P}}(k)\). Here, \({\hat{P}}(k)\) maps the outgoing wave amplitudes on the edges \({\varvec{a}}^{+}_{E}\) to the incoming wave amplitudes \({\varvec{a}}^{-}_{E}\), taking account of the phase modulation owing to the metric \({\mathcal {L}}\), and \({\hat{S}}_{EE}\) represents the vertex scattering between edges. Fourth and finally, the waves on the edges are coupled back onto the leads using a matrix \({\hat{S}}_{LE}\). Explicitly, we can write

$$\begin{aligned} {\hat{S}}_{\Gamma }(k) = {\hat{S}}_{LL} + {\hat{S}}_{LE}\left[ {\hat {{\mathbb{I}}}} -{\hat{P}}(k){\hat{S}}_{EE}\right] ^{-1}{\hat{P}}(k){\hat{S}}_{EL}, \end{aligned}$$
(6)

see Kottos and Smilansky17 for details. The matrix terms in \({\hat{S}}_{\Gamma }\) can be deduced from the open graph scattering matrix \({\hat{S}}\) in (3) using the block-matrix representation,

$$\begin{aligned} \begin{pmatrix} {\varvec{a}}_{L}^{+} \\ \hline {\varvec{a}}_{E}^{+} \end{pmatrix} = \left( \begin{array}{c|c} {\hat{S}}_{LL} &{}\quad {\hat{S}}_{LE} \\ \hline {\hat{S}}_{EL} &{}\quad {\hat{S}}_{EE} \\ \end{array} \right) \begin{pmatrix} {\varvec{a}}_{L}^{-} \\ \hline {\varvec{a}}_{E}^{-} \end{pmatrix} \end{aligned}$$
(7)

In the next step, we will use this graph based resonator to construct a periodic metamaterial.

A quantum graph model for a metamaterial

By periodically arranging these graph based resonant elements we will now model metamaterial behaviour. Consider placing each resonant element in a mesh, spaced by edges of length l. To keep track of the wave location in the graph, indices nm are introduced, where n representing steps in the horizontal direction and m in the vertical, as illustrated in Fig. 3. We use z again for the coordinate on the edges with \(z=0\) at a given resonant element and \(z = l\) at the neighbouring element for all edge (where we leave out the indices nm for convenience).

Figure 3
figure 3

A 2D square periodic arrangement of quantum graph based resonant elements \(\Gamma (V,E,L)\) connected by edges of length l as a model for a 2D metamaterial. Location \(n = m = 0\) shown in the subplot with labeled wave amplitudes, \({\varvec{a}}^\pm\) and \({\varvec{c}}^-\), traveling on the edges in the neighbourhood of the central vertex.

Full size image

As before, all graph edges are endowed with the one dimensional wave equation (1). The vector of solutions is restricted to the edges (l, r, d, u) in each unit cell of the material at location nm,

$$\begin{aligned} \varvec{\psi }_{nm}(z) = \begin{pmatrix} \psi _{nm,l}(z) \\ \psi _{nm,r}(z) \\ \psi _{nm,d}(z) \\ \psi _{nm,u}(z) \end{pmatrix}. \end{aligned}$$
(8)

To solve for the wave solution across the metamaterial, we begin by evaluating the scattering at the central vertex \(n=m=0\). Let \({\varvec{a}}^\pm = (a_l^\pm , a_r^\pm , a_d^\pm , a_u^\pm )^T\) be the vector of wave amplitudes heading in \((-)\) or out \((+)\) of the central vertex, see Fig. 3. We have

$$\begin{aligned} {\varvec{a}}^+ = {\hat{S}}_{\Gamma }(k) {\varvec{a}}^-. \end{aligned}$$
(9)

The waves once scattered, undergo phase modulation as they move along each edge allowing one to express the outgoing waves at the central vertex \({\varvec{a}}^+\) in terms of the incoming wave amplitudes \({\varvec{c}}^-= (c_l^- , c_r^- , c_d^- , c_u^-)^T\) at the neighboring vertices,

$$\begin{aligned} {\varvec{c}}^- = e^{ikl} {\varvec{a}}^+. \end{aligned}$$
(10)

Due to the periodicity of the structure we can use Bloch’s theorem to describe the wave solution in any unit cell nm in terms of the solution at the central vertex \(nm=00\)18, that is,

$$\begin{aligned} \varvec{\psi }_{nm}(z) = e^{i(\kappa _{x}n + \kappa _{y}m)l} \varvec{\psi }_{00}(z). \end{aligned}$$
(11)

Here, \(\kappa _{x}\) and \(\kappa _{y}\) are the Bloch wave numbers in the horizontal and vertical direction, respectively. The incoming wave amplitudes at a given vertex and its neighbours are thus related by

$$\begin{aligned} {\varvec{a}}^- = {\hat{B}}(\kappa _x,\kappa _y) {\varvec{c}}^- \end{aligned}$$
(12)

where

$$\begin{aligned} {\hat{B}}(\kappa _x,\kappa _y) := \begin{pmatrix} 0 &{}\quad e^{-i\kappa _xl} &{}\quad 0 &{}\quad 0 \\ e^{ i\kappa _xl} &{}\quad 0 &{}\quad 0 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad 0 &{}\quad e^{-i\kappa _yl} \\ 0 &{}\quad 0 &{}\quad e^{ i\kappa _yl} &{}\quad 0 \\ \end{pmatrix}. \end{aligned}$$
(13)

By combining Eqs. (9), (10) and (12), we find the condition

$$\begin{aligned} \left[ {\hat{\mathbb{I}}} - e^{ikl}{\hat{B}}(\kappa _x,\kappa _y) {\hat{S}}_{\Gamma }(k)\right] {\varvec{a}}^- = {\varvec{0}}, \end{aligned}$$
(14)

which is satisfied if the wave vector k on each edge and the quasi-wave vectors \(\kappa _x\) and \(\kappa _y\) fulfill the secular equation

$$\begin{aligned} \det \left[ {\hat{\mathbb{I}}} - e^{ikl}{\hat{B}}(\kappa _x,\kappa _y){\hat{S}}_{\Gamma }(k)\right] = 0. \end{aligned}$$
(15)

In satisfying this condition, the incoming amplitudes at the central vertex \({\varvec{a}}^-\) can be found by solving for the eigenvector in Eq. (14). We then obtain the wave solutions in each unit cell nm throughout the metamaterial in the form

$$\begin{aligned} \varvec{\psi }_{nm}(z;k,\kappa _x,\kappa _y) = e^{i(\kappa _{x}n + \kappa _{y}m)l}\big ({\varvec{a}}^+(k,\kappa _x,\kappa _y)e^{ikz} + {\varvec{a}}^-(k,\kappa _x,\kappa _y)e^{-ikz}\big ), \end{aligned}$$
(16)

where \({\varvec{a}}^+\) is obtained via (9).

Wave solutions for the infinite lattice

In the following, we will show how to construct plane wave and then beam-like solutions for the infinite lattice. The starting point for the construction are the solutions of Eq. (15) giving rise to dispersion curves of the form \(k(\kappa _x, \kappa _y)\) or \(\kappa _x(k, \kappa _y)\). We encounter real or imaginary solutions of \(\kappa _x(k, \kappa _y)\) corresponding to propagating or evanescent waves, respectively, the latter occurring in band gaps. The full wave solution can then be constructed from the corresponding eigenvectors in (14) extended throughout the lattice using the Bloch condition (11), as formulated in (16).

Wave propagation in metamaterials—plane wave solutions

We start by constructing plane wave solutions and will discuss a particular simple case here, where the dispersion curves can be given analytically. That is, we consider a lattice without resonators where each edge is directly connected through the vertex. Here, the vertex acts as a point scatterer, where the scattering matrix \({\hat{S}}_{\Gamma }(k)\) in (6) is given by the Neumann scattering matrix (4) with \(v_V = 4\), that is,

$$\begin{aligned} {\hat{S}}_{\Gamma _{point}}= \frac{1}{2} \begin{pmatrix} -1 &{}\quad 1 &{}\quad 1 &{}\quad 1 \\ 1 &{}\quad -1 &{}\quad 1 &{}\quad 1 \\ 1 &{}\quad 1 &{}\quad -1 &{}\quad 1 \\ 1 &{}\quad 1 &{}\quad 1 &{}\quad -1 \end{pmatrix}. \end{aligned}$$
(17)

A generalisation to arbitrary \({\hat{S}}_{\Gamma }(k)\) is straightforward, examples will be shown in later sections. By solving Eq. (15), the dispersive properties of the lattice are given by the relation

$$\begin{aligned} 2\text{ cos }(kl) = \text{ cos }(\kappa _{x}l) + \text{ cos }(\kappa _{y}l). \end{aligned}$$
(18)
Figure 4
figure 4

The band diagram for a lattice without resonant elements. Shown are iso-frequency contours for various values of k in the \(\kappa _{x}\), \(\kappa _{y}\) plane.

Full size image

The resulting surface, or band, is shown in Fig. 4; it is periodic in \(\varvec{\kappa }\)-space over the Brillouin Zone (BZ), where \(BZ \in [-\frac{\pi }{l}, \frac{\pi }{l}]^{2}\). The shape of the band is a function of the square periodic graph topology and the symmetry of the scattering matrix. To generate a solution consider taking a single value of k, resulting in a discrete ring of \(\varvec{\kappa }\) solutions on an iso-frequency contour. By picking a point on the contour, one constructs the full wave solution from the corresponding eigenfunction in Eq. (16) as plotted in Fig. 5.

Figure 5
figure 5

Iso-frequency contour of Fig. 4 with possible wave vectors \(\varvec{\kappa } = (\kappa _{x}^{\rightarrow },\kappa _{y}')^{T}\) in (a) and \(\varvec{\kappa } =(\kappa _{x}^{\leftarrow },\kappa _{y}')^{T}\) in (b), both shown in green, with corresponding Poynting vector \({\varvec{J}} = (J_{x},J_{y})^{T}\) and \({\varvec{J}} = (-J_{x},J_{y})^{T}\) normal to the contour shown in blue. (See next section for the choice of notation). The resulting real components of the eigenfunction solutions \(\varvec{\psi }_{nm}^{\rightleftharpoons }(z;k,\kappa _{x},\kappa _{y})\) for the values of the wave vectors used in (a) and (b) are also shown.

Full size image

The phase of the resulting eigenfunction is given by the chosen Bloch wave vector \(\varvec{\kappa } = (\kappa _{x},\kappa _y)^{T}\), while the energy flow is given by the Poynting vector \({\varvec{J}} = (J_x,J_y)^T\) normal to the iso-frequency contour. To construct the components of \({\varvec{J}}\) explicitly, one evaluates the flux on the graph edges. The 1D flux J of the wave function \(\varvec{\psi }_{nm}\) on edge p is

$$\begin{aligned} J(\psi _{nm,p}(z)) := \Re \left( \bar{\psi }_{nm,p}(z)\frac{1}{i}\frac{\partial \psi _{nm,p}(z)}{\partial z}\right) = k(|a_{p}^{+}|^{2} - |a_{p}^{-}|^{2}) \end{aligned}$$
(19)

with \(\bar{\psi }\) denoting the complex conjugate of \(\psi\). The horizontal and vertical components of the Poynting vector can then be evaluated in terms of the waves on edges right(r)(or left(l)) and up(u)(or down(d)), respectively, that is

$$\begin{aligned} {\varvec{J}} = \begin{pmatrix} J_{x} \\ J_{y} \end{pmatrix} := k \begin{pmatrix} |a_{r}^{+}|^{2} - |a_{r}^{-}|^{2} \\ |a_{u}^{+}|^{2} - |a_{u}^{-}|^{2} \end{pmatrix} = k \begin{pmatrix} |a_{l}^{-}|^{2} - |a_{l}^{+}|^{2} \\ |a_{d}^{-}|^{2} - |a_{d}^{+}|^{2} \end{pmatrix}. \end{aligned}$$
(20)

A point on notation

In the example shown in Fig. 5, for every value of \(\kappa _{y}\), there are two corresponding values of \(\kappa _{x}\). Naturally, the choice of \(\kappa _{x}\) informs the direction of energy flow. To delineate between waves traveling in opposite horizontal directions, the following notation is used: eigenfunction solutions with Poynting vector heading to the right are given index \(\rightarrow\) and the corresponding eigenvector in Eq. (14) is expressed as \({\varvec{a}}^{-}\) as before where \(\varvec{a^{\pm }} := \varvec{a^{\pm }}(k,\kappa _{x}^{\rightarrow },\kappa _{y})\); eigenfunction solutions with Poynting vector heading to the left are given index \(\leftarrow\) with the corresponding eigenvector in (14) relabeled as \({\varvec{b}}^{-}\) where \({\varvec{b}}^{\pm } := {\varvec{a}}^{\pm }(k,\kappa _{x}^{\leftarrow },\kappa _{y})\). Explicitly,

$$\begin{aligned} \varvec{\psi }_{nm}(z;k,\kappa _{x},\kappa _{y}) := {\left\{ \begin{array}{ll} \varvec{\psi }_{nm}^{\rightarrow }(z;k,\kappa _{y})=e^{i(\kappa _{x}^{\rightarrow }n + \kappa _{y}m)l}\left( {\varvec{a}}^{+}e^{ikz} + {\varvec{a}}^{-}e^{-ikz}\right) ,&{}\quad J_{x} > 0\\ \varvec{\psi }_{nm}^{\leftarrow }(z;k,\kappa _{y})=e^{i(\kappa _{x}^{\leftarrow }n + \kappa _{y}m)l}\left( {\varvec{b}}^{+}e^{ikz} + {\varvec{b}}^{-}e^{-ikz}\right) ,&{}\quad J_{x} < 0\\ \end{array}\right. }. \end{aligned}$$
(21)

In choosing the value of \(\kappa _{x}\) (\(\kappa _{x}^{\rightarrow }\) or \(\kappa _{x}^{\leftarrow }\)), one implicitly chooses the wave direction. For this reason the \(\kappa _{x}\) dependence is dropped from the eigenfunction.

Wave propagation in metamaterials—Gaussian beams

Plane wave solutions travel in the direction of the wave vector \(\varvec{\kappa }\), which does not necessarily agree with the direction of energy flow given by the Poynting vector \({\varvec{J}}\) in (20). This can be observed when considering Gaussian beam solutions constructed via Fourier series with the eigenfunctions, shown above, as a basis. The solution of the Gaussian beam with focal point at \(n=n'\), expressed in components \(\varvec{\Phi }_{nm}(z) = (\Phi _{nm,l}(z),\Phi _{nm,r}(z),\Phi _{nm,d}(z),\Phi _{nm,u}(z))^{T}\), is given as,

$$\begin{aligned} \varvec{\Phi }_{nm}(z; k) = \frac{1}{\sqrt{2\pi }}\int _{\Omega }\alpha _{n'}(\kappa _y;k) \varvec{\psi }_{nm}^{\rightarrow }(z,\kappa _{y};k)d\kappa _{y}, \end{aligned}$$
(22)

where the integral is performed over the domain \(\Omega = \Omega (k)\) of the iso-frequency contour. Solutions outside \(\Omega\) are evanescent and so contribute nothing in the far field. The expansion coefficients \(\alpha _{n'}\) describing the beam profile at the focal point in terms of the eigenfunctions, given by the inverse transform,

$$\begin{aligned} \begin{aligned} \alpha _{n'}(\kappa _{y}; k) = \frac{1}{\sqrt{2\pi }}\sum _{m=-\infty }^{\infty }\Bigg \{&\int _{0}^{l}\bar{\psi }_{n'm,u}^{\rightarrow }(z,\kappa _y;k)\Phi _{n'm,u}(z;\kappa _{y})dz\Bigg \} \\ \end{aligned}. \end{aligned}$$
(23)

Here, \(\bar{\psi }_{n'm,u}^{\rightarrow }\) is the complex conjugate of the right moving eigenfunction expressed for all m at horizontal location \(n = n'\) along the upward edges u and \(\Phi _{n'm,u}(z)\) is the beam profile for all m with focal point \(n = n'\) expressed on edge u. We chose the beam profile to be Gaussian, such that

$$\begin{aligned} \Phi _{n'm,u}(z) = \frac{1}{\sqrt{\sigma \sqrt{\pi }}} e^{-\left( \frac{z + ml}{\sqrt{2}\sigma }\right) ^{2}} e^{i\kappa _{y}'l}, \end{aligned}$$
(24)

where \(\sigma\) is the width of the beam and the phase \(e^{i\kappa _{y}'l}\) determines the tilt angle of the beam with respect to the horizontal axis of the lattice. Since the space on the graph is discretised, so too are the spacial integrals. Figure 6 illustrates how varying the parameters \(\sigma\) and \(\kappa _y'\) effect the shape and direction of the resulting Gaussian beam for a given iso-frequency contour \(\Omega (k)\).

Figure 6
figure 6

Iso-frequency contour for a metamaterial with resonant elements \({\hat{S}}_{\Gamma _{point}}\) as defined in (17) at \(k = 1/l\) together with the expansion coefficients \(\alpha _{n'=0}\) in Eq. (23) and the resulting real components of the Gaussian beam profile \(\varvec{\Phi }_{nm}\) with focal point set to \(n' = 0\). The beam profiles, Eq. (24), shown here are characterised by the parameters \(\kappa _{y}' = 0\), 0 and 1/l and beam widths \(\sigma = 2.2l\), 6.6l and 6.6l for (a), (b) and (c), respectively.

Full size image

Band engineering

In a next step, we will use the quantum graph formulation for designing resonant elements leading to materials with “exotic” wave properties. As a simple example, we consider using the cross resonator illustrated in Fig. 1a as a resonant element for the metamaterial illustrated in Fig. 3. By setting the metric of the cross resonator to \({\mathcal {L}} = \{l_{x}/2,l_{x}/2,l_{y}/2,l_{y}/2\}\) as shown in Fig. 7a, one can modulate the phase in the horizontal and vertical direction across the resonant element thus breaking the scattering symmetry between the four leads. By enforcing Neumann boundary conditions at the vertices as defined by Eq. (4), the scattering matrix between the four leads L has solution,

$$\begin{aligned} {\hat{S}}_{\Gamma _{cross}}(k;l_{x},l_{y}) = \frac{1}{2} \begin{pmatrix} -e^{ikl_{x}} &{}\quad e^{ikl_{x}} &{}\quad e^{ik\left( \frac{l_{x} + l_{y}}{2}\right) } &{} \quad e^{ik\left( \frac{l_{x} + l_{y}}{2}\right) } \\ e^{ikl_{x}} &{} \quad -e^{ikl_{x}} &{}\quad e^{ik\left( \frac{l_{x} + l_{y}}{2}\right) } &{}\quad e^{ik\left( \frac{l_{x} + l_{y}}{2}\right) } \\ e^{ik\left( \frac{l_{x} + l_{y}}{2}\right) } &{}\quad e^{ik\left( \frac{l_{x} + l_{y}}{2}\right) } &{}\quad -e^{ikl_{y}} &{}\quad e^{ikl_{y}} \\ e^{ik\left( \frac{l_{x} + l_{y}}{2}\right) } &{}\quad e^{ik\left( \frac{l_{x} + l_{y}}{2}\right) } &{}\quad e^{ikl_{y}} &{}\quad -e^{ikl_{y}} \end{pmatrix}, \end{aligned}$$
(25)

found by Eq. (6). By solving Eq. (15) the dispersion relation is now given as

$$\begin{aligned} \hbox {cos}(\kappa _{x}l)e^{ik(l_{x}+l)}\left( e^{2ik(l_{y}+l)} - 1\right) + \hbox {cos}(\kappa _{y}l)e^{ik(l_{y}+l)}\left( e^{2ik(l_{x}+l)} - 1\right) -e^{2ik(l_{x} + l_{y} + 2l)} + 1 = 0. \end{aligned}$$
(26)

By comparing Fig. 7b with 7c, we see that varying the phase modulation across each resonant element only in one dimension breaks the symmetry of the scattering and gives rise to a saddle shaped band, which can be exploited to yield negative refraction as shown in the next section.

Figure 7
figure 7

(a) Shows the cross resonator plugged into the metamaterials unit cell with scattering matrix defined as \({\hat{S}}_{\Gamma _{cross}}(k;l_x,l_y)\) in (25). Both (b) and (c) show the resulting band for different values of \(l_x\) and \(l_y\). Also plotted in the \(\kappa _{x}\), \(\kappa _{y}\) plane are the iso-frequency contours of the second band for various values of k. Plot (b) is for \(l_{x} = 1\), \(l_{y} = 1\). Plot (c) is for \(l_{x} = 0\), \(l_{y} = 1\).

Full size image

Wave refraction between metamaterials

To exemplify the engineered refractive properties, consider two metamaterials represented by semi infinite square periodic quantum graphs, one with resonant elements given by \({\hat{S}}_{\Gamma _{1}}\) and the other with resonant elements given by \({\hat{S}}_{\Gamma _{2}}\). The two materials are connected along the y direction and material 1 exists for the set of unit cells \(n_{1} = (-\infty ,\ldots ,-2,-1\}\) and material 2 exists for the set of unit cells \(n_{2} = \{0,1,\ldots ,\infty )\), as illustrated in Fig. 8.

Figure 8
figure 8

Boundary region between metamaterials 1 and 2, understood as all right(r) edges for \(n = -1\) and all left (l) edges for \(n = 0\). Here, wave scattering from the boundary is divided into event I and II, where \(r_{p}\) and \(t_{p}\) represent reflection and transmission amplitudes for event \(p = I\) or II.

Full size image

The full wave solution \(\varvec{\Psi }_{nm}= (\Psi _{nm,l},\Psi _{nm,r},\Psi _{nm,d},\Psi _{nm,u})^{T}\) across the two materials can be constructed by a linear superposition of counter propagating eigenfunction solutions \(\varvec{\psi }_{j,nm}^{\rightleftharpoons }\) in material \(j = 1\) and 2, as expressed in Eq. (21), that is,

$$\begin{aligned} \begin{aligned} \varvec{\Psi }_{nm}(z,\kappa _{y};k)&= H({n}_{1})\left[ A_{1}(\kappa _{y};k)\varvec{\psi }_{1,nm}^{\rightarrow }(z,\kappa _{y};k) + B_{1}(\kappa _{y};k)\varvec{\psi }_{1,nm}^{\leftarrow }(z,\kappa _{y};k) \right] \\&\quad + H({n}_{2})\left[ A_{2}(\kappa _{y};k)\varvec{\psi }_{2,nm}^{\rightarrow }(z,\kappa _{y};k) + B_{2}(\kappa _{y};k)\varvec{\psi }_{2,nm}^{\leftarrow }(z,\kappa _{y};k) \right] , \end{aligned} \end{aligned}$$
(27)

where \(H({n}_{j})\) is the discretised Heavyside step function, that is,

$$\begin{aligned} H({n}_{j}) = {\left\{ \begin{array}{ll} 1,&{}\quad \forall n \in {n_{j}} \\ 0,&{}\quad \forall n \notin {n_{j}} \\ \end{array}\right. } \end{aligned}$$
(28)

for \(j=1,2\) and the coefficients \(A_j\) and \(B_j\) are associated with left and right moving waves, respectively. Thus, solutions with coefficients \(A_{1}\) and \(B_{2}\) represent waves incident on the interface, while solutions with coefficients \(B_{1}\) and \(A_{2}\) represent waves scattered from the interface. To determine the coefficients, we must satisfy the boundary conditions at the material interface. Naturally the wave solutions are given for a single value of k between the two materials enforcing \(k_{1} = k_{2} :=k\). As the system stays periodic in the y direction, the Bloch phase tangential to the interface also remain constant across the interface leading to the boundary condition \(\kappa _{1,y} = \kappa _{2,y} :=\kappa _y\). This condition is illustrated in Fig. 9 for a chosen value of \(\kappa _{y} = \kappa _{y}'\) by a horizontal dashed grey line connecting the iso-frequency contours at \(k = 1/l\) in material 1 and 2. The remaining unknowns are then the wave vectors normal to the interface \(\kappa _{j,x}^{\rightleftharpoons }(k,\kappa _{y}')\) which can be obtained from the dispersion curve determined by solving Eq. (15) for each material, see the vertical dashed grey lines illustrated below.

Figure 9
figure 9

The iso-frequency contours for \(k = 1/l\) for a scattering matrix with \({\hat{S}}_{\Gamma _{1}} := {\hat{S}}_{\Gamma _{point}}\) and \({\hat{S}}_{\Gamma _{2}} := {\hat{S}}_{\Gamma _{cross}}(k;0,1)\) representing the two metamaterials. The horizontal dashed grey line represents a single value of \(\kappa _{y} = \kappa _{y}'\). Intersections with contour lines give the corresponding values of \(\kappa _{j,x}^{\rightleftharpoons }\) for \(j=1\) and 2.

Full size image

To determine the scattering coefficients \(A_j\) and \(B_j\), we rescale the eigenfunction solutions of each material such that the magnitude of the horizontal flux of each component is equal, that is,

$$\begin{aligned} \begin{aligned}{}&J(\psi _{1,nm,l}^{\rightarrow }) = J(\psi _{2,nm,l}^{\rightarrow }) = -J(\psi _{1,nm,l}^{\leftarrow }) = -J(\psi _{2,nm,l}^{\leftarrow }) \end{aligned}, \end{aligned}$$
(29)

assuming \(\kappa _{j,x}^{\rightleftharpoons }\in \mathbb{R}\). Equally, the same condition can be enforced on the right edges, r. With the scaling choice (29), flux conservation across the interface

$$\begin{aligned} J\left( A_1{\psi }_{1,0m,l}^{\rightarrow } + B_{1}{\psi }_{1,0m,l}^{\leftarrow }\right) = J\left( A_2{\psi }_{2,0m,l}^{\rightarrow } + B_{2}{\psi }_{2,0m,l}^{\leftarrow }\right) \end{aligned}$$
(30)

reduces to

$$\begin{aligned} |A_{1}|^{2} + |B_{2}|^{2} = |B_{1}|^{2} + |A_{2}|^{2} \end{aligned}$$
(31)

and wave scattering at the interface can then be described in terms of a unitary scattering process. The corresponding interface scattering matrix \({\hat{S}}_{1,2}\) performing the mapping,

$$\begin{aligned} \begin{pmatrix} B_{1}\\ A_{2}\\ \end{pmatrix} ={\hat{S}}_{1,2} \begin{pmatrix} A_{1}\\ B_{2}\\ \end{pmatrix} \end{aligned}$$
(32)

can then be constructed by decomposing the interface scattering into two events. Event I describes a wave incident from material 1 onto material 2 with amplitude \(A_{1} = 1\) and \(B_{2} = 0\), producing a reflected and transmitted wave with respective amplitudes \(r_I\) and \(t_{I}\). Event II describes a wave incident from material 2 onto material 1 with amplitude \(B_{2} = 1\) and \(A_{1} = 0\), producing a reflected and transmitted wave with respective amplitudes \(r_{II}\) and \(t_{II}\), see Fig. 8. The interface scattering matrix takes on the form

$$\begin{aligned} {\hat{S}}_{1,2} = \begin{pmatrix} r_{I} &{}\quad t_{II} \\ t_{I} &{}\quad r_{II} \end{pmatrix}. \end{aligned}$$
(33)

To evaluate the matrix elements, consider first event I. For simplicity we choose to evaluate the waves at location \(n = 0\), at coordinate \(z = 0\) on edge l. Since the phase \(\kappa _{y}\) is the same in both materials, it is sufficient to evaluate the solutions at \(m = 0\). Here,

$$\begin{aligned} \begin{aligned}{}&{\Psi }_{(\forall n \le -1)0,l}(0)|_{n=0} = \left( a_{1,l}^{+} + a_{1,l}^{-}\right) + r_{I}\left( b_{1,l}^{+} + b_{1,l}^{-}\right) \\&{\Psi }_{(\forall n \ge 0)0,l}(0)|_{n=0} = t_{I}\left( a_{2,l}^{+} + a_{2,l}^{-}\right) . \end{aligned} \end{aligned}$$
(34)

\({\Psi }_{(\forall n \le -1)0,l}(0)|_{n=0}\) and \({\Psi }_{(\forall n \ge 0)0,l}(0)|_{n=0}\) both exist on the same edge in the boundary region, so to stay consistent, there must be an equivalence between the incoming(−) and outgoing(\(+\)) wave amplitudes of these solutions, that is,

$$\begin{aligned} \begin{aligned}{}&t_{I}a_{2,l}^{+} = a_{1,l}^{+} + r_{I}b_{1,l}^{+} \\&t_{I}a_{2,l}^{-} = a_{1,l}^{-} + r_{I}b_{1,l}^{-}. \end{aligned} \end{aligned}$$
(35)

This can be solved to give

$$\begin{aligned} \begin{aligned}{}&r_{I} = \frac{a_{2,l}^{+}a_{1,l}^{-} - a_{2,l}^{-}a_{1,l}^{+}}{a_{2,l}^{-}b_{1,l}^{+} - b_{1,l}^{-}a_{2,l}^{+}} \\&t_{1} =\frac{a_{1,l}^{-}b_{1,l}^{+} - b_{1,l}^{-}a_{1,l}^{+}}{a_{2,l}^{-}b_{1,l}^{+} - b_{1,l}^{-}a_{2,l}^{+}}. \end{aligned} \end{aligned}$$
(36)

Exactly the same procedure can be done for event II where

$$\begin{aligned} \begin{aligned}{}&{\Psi }_{(\forall n \le -1)0,l}(0)|_{n=0} = t_{II}\left( b_{1,l}^{+} + b_{1,l}^{-}\right) \\&{\Psi }_{(\forall n \ge 0)0,l}(0)|_{n=0} = \left( b_{2,l}^{+} + b_{2,l}^{-}\right) + r_{II}\left( a_{2,l}^{+} + a_{2,l}^{-}\right) , \end{aligned} \end{aligned}$$
(37)

which yields the equivalence condition

$$\begin{aligned} \begin{aligned}{}&t_{II}b_{1,l}^{+} = b_{2,l}^{+} + r_{II}a_{2,l}^{+} \\&t_{II}b_{1,l}^{-} = b_{2,l}^{-} + r_{II}a_{2,l}^{-} \end{aligned} \end{aligned}$$
(38)

with solutions

$$\begin{aligned} \begin{aligned}{}&r_{II} =\frac{b_{1,l}^{+}b_{2,l}^{-} - b_{1,l}^{-}b_{2,l}^{+}}{a_{2,l}^{+}b_{1,l}^{-} - a_{2,l}^{-}b_{1,l}^{+}} \\&t_{II} = \frac{a_{2,l}^{+}b_{2,l}^{-} - a_{2,l}^{-}b_{2,l}^{+}}{a_{2,l}^{+}b_{1,l}^{-} - a_{2,l}^{-}b_{1,l}^{+}}. \end{aligned} \end{aligned}$$
(39)

The results of this unitary interface scattering is shown in Fig. 10 for a Gaussian beam incident from material 1 with amplitude \(A_{1} = 1\) and no incident beam from material 2, that is, \(B_{2} = 0\).

Figure 10
figure 10

A refracted beam incident on a material interface with material properties defined by \({\hat{S}}_{\Gamma _{1}} := {\hat{S}}_{\Gamma _{point}}\) and \({\hat{S}}_{\Gamma _{2}} := {\hat{S}}_{\Gamma _{cross}}(k;0,1)\).

Full size image

The Gaussian Beam is constructed based on the method described in section “Wave propagation in metamaterials—Gaussian beams”, where the basis is the full wave field \(\varvec{\Psi }_{nm}(x,\kappa _{y};k)\) as defined in Eq. (27),

$$\begin{aligned} \varvec{\Phi }_{nm}(z) = \frac{1}{\sqrt{2\pi }}\int _{\Omega }\alpha _{_{-60}}(\kappa _y)\varvec{\Psi }_{nm}(z,\kappa _{y};k)d\kappa _{y}. \end{aligned}$$
(40)

The expansion coefficients \(\alpha _{_{-60}}\) are as in Eq. (23) for a beam focal point \(n' = -60\) and tilt given by \(\kappa _{y}' = 1/l\).

As was shown in section “Band engineering”, breaking the symmetry of the scattering matrix gives rise to a saddle shaped band resulting in the property of negative refraction at appropriate \(\kappa _y\) values, see Fig. 7c. When choosing such a material on the right hand side, one obtains the reflection/transmission behaviour shown in Fig. 11.

Figure 11
figure 11

Top: Iso-frequency contours for \(k = 1/l\) for a scattering matrix with \({\hat{S}}_{\Gamma _{1}} := {\hat{S}}_{\Gamma _{point}}\) and \({\hat{S}}_{\Gamma _{2}} := {\hat{S}}_{\Gamma _{cross}}(k;0,4.45)\). Bottom: The resulting real component of the Gaussian beam \(\varvec{\Phi }_{nm}\) constructed from the full wave field \(\varvec{\Psi }_{nm}\) across the two materials with incident wave amplitudes \(A_{1} = 1\) and \(B_{2} = 0\). The focal point of the Gaussian Beam is set to \(n' = -60\) with tilt given by \(\kappa _{y}' = 1/l\).

Full size image

Waves in N layered metamaterials

Consider now a system of N layered metamaterials as shown in Fig. 12, each with their own properties as defined previously.

Figure 12
figure 12

A system of N layered metamaterials with wave amplitudes, \(A_{j}\) and \(B_{j}\), noted at each interface.

Full size image

Each material spans a domain defined by the set of unit cells \(n_{j}\). Here the wave function across all N materials is expressed as a linear superposition of counter propagating eigenfunction solutions of each material, that is,

$$\begin{aligned} \varvec{\Psi }_{nm}(x,\kappa _{y};k) = \sum _{j = 1}^{N}H(n_{j}) \left[ A_{j}(\kappa _{y};k)\varvec{\psi }_{j,nm}^{\rightarrow }(x,\kappa _y;k) + B_{j}(\kappa _{y};k)\varvec{\psi }_{j,nm}^{\leftarrow }(x,\kappa _y;k)\right] . \end{aligned}$$
(41)

To evaluate the full wave function across all N materials, one must determine the coefficients \(A_{j}\) and \(B_{j}\) that satisfy the boundary conditions at all material interfaces. This is done by the Transfer Matrix Method19. In the previous section, the wave amplitudes between two materials were determined by constructing a scattering matrix mapping incoming to outgoing wave amplitudes between the metamaterials 1 and 2 as defined in Eq. (32). This matrix can be rearrange to give a transfer matrix that maps wave amplitudes from material 1 to material 2, \(\hat{\tilde{S}}_{1,2}\).

$$\begin{aligned} \begin{pmatrix} A_{2}\\ B_{2}\\ \end{pmatrix} =\hat{\tilde{S}}_{1,2} \begin{pmatrix} A_{1}\\ B_{1}\\ \end{pmatrix} = \begin{pmatrix} t_{I} - \frac{r_{I}r_{II}}{t_{II}} &{}\quad \frac{r_{II}}{t_{II}}\\ -\frac{r_{I}}{t_{II}} &{}\quad \frac{1}{t_{II}}\\ \end{pmatrix} \begin{pmatrix} A_{1}\\ B_{1}\\ \end{pmatrix}. \end{aligned}$$
(42)

This procedure can be generalised to any arbitrary materials j and \(j+1\) giving

$$\begin{aligned} \begin{pmatrix} A_{j+1}\\ B_{j+1}\\ \end{pmatrix} =\hat{\tilde{S}}_{j,j+1} \begin{pmatrix} A_{j}\\ B_{j}\\ \end{pmatrix}. \end{aligned}$$
(43)

Waves propagating across a given material j accumulate a Bloch phase which can be expressed in terms of the matrix

$$\begin{aligned} {\hat{P}}_{j}(\kappa _{y};k) = \begin{pmatrix} e^{i\kappa _{j,x}^{\rightarrow }(\kappa _{y};k)W_{j}} &{}\quad 0 \\ 0 &{}\quad e^{i\kappa _{j,x}^{\leftarrow }(\kappa _{y};k)W_{j}}\\ \end{pmatrix}, \end{aligned}$$
(44)

where \(W_{j}=(\text{ max }(n_{j}) - \text{ min }(n_{j}))l\) is the width of material j. Having now formulated both scattering and propagation, one can express the wave amplitudes in material N in terms of the wave amplitudes in material 1.

$$\begin{aligned} \begin{pmatrix} A_{N}\\ B_{N}\\ \end{pmatrix} = \hat{\tilde{S}}_{N-1,N}\left( {\hat{P}}_{N-1}\hat{\tilde{S}}_{N-2,N-1}\right) \ldots \left( {\hat{P}}_{3}\hat{\tilde{S}}_{2,3}\right) \left( {\hat{P}}_{2}\hat{\tilde{S}}_{1,2}\right) \begin{pmatrix} A_{1}\\ B_{1}\\ \end{pmatrix} := \hat{\tilde{S}}_{1,N} \begin{pmatrix} A_{1}\\ B_{1}\\ \end{pmatrix}. \end{aligned}$$
(45)

We can now rearrange this transfer operator \(\hat{\tilde{S}}_{1,N}\) such that the entire system of layered materials acts as a single point scatterer by introducing the scattering matrix \({\hat{S}}_{1,N}\) defined as

$$\begin{aligned} \begin{pmatrix} B_{1}\\ A_{N}\\ \end{pmatrix} = {\hat{S}}_{1,N} \begin{pmatrix} A_{1}\\ B_{N}\\ \end{pmatrix}. \end{aligned}$$
(46)

Now by setting the incident wave amplitudes \(A_{1}\) and \(B_{N}\), we determine the amplitudes of the scattered field by direct substitution into Eq. (46). Knowing the wave amplitudes in material 1, \(A_{1}\) and \(B_{1}\), it is trivial to determine all other amplitudes using \({\hat{P}}_{j-1}\hat{\tilde{S}}_{j,j+1}\). The results of this procedure is plotted in Fig. 13 for a three layered material.

Figure 13
figure 13

Top: The iso-frequency contours of three metamaterials with properties defined by scattering matrices \({\hat{S}}_{\Gamma _{point}}\), \({\hat{S}}_{\Gamma _{cross}}(k;0,4.45)\) and \({\hat{S}}_{\Gamma _{point}}\) for \(k = 1/l\). Bottom: The real component of a Gaussian beam \(\varvec{\Phi }_{nm}\) incident from material 1 constructed from the full wave field \(\Psi _{nm}\) for an incident wave from metamaterial 1, \(A_{1} = 1\), and no incident wave from material 3, \(B_{3} = 0\).

Full size image

Conclusion

We propose here a graph based technique for designing metamaterials allowing for a fast and flexible way to test resonant element proposals embedded periodically within a 2D square lattice. By modelling resonant elements in terms of open, graph-based scattering systems, we retain a connection with the underlying geometrical structure of the element which will inspire the engineering of physical metamaterials. This will make it possible to search for exotic band-diagrams, dispersion curves and wave effects linking these to an underlying array of resonant scattering systems. In this paper, we introduce the principles of graph based metamaterial construction including the set-up and we show how to obtain plane wave and beam-like solutions in the infinite medium. We then demonstrate the handling of boundary conditions at interfaces as well as the generalisation to N layered media. In each case, the computation can be reduced to low-dimensional matrix problems with the help of Bloch’s theorem. Further results modelling various special wave effects will be presented in forthcoming publications.