Parabolic Metamaterials and Dirac Bridges

A new class of multi-scale structures, referred to as `parabolic metamaterials' is introduced and studied in this paper. For an elastic two-dimensional triangular lattice, we identify dynamic regimes, which corresponds to so-called `Dirac Bridges' on the dispersion surfaces. Such regimes lead to a highly localised and focussed unidirectional beam when the lattice is excited. We also show that the flexural rigidities of elastic ligaments are essential in establishing the `parabolic metamaterial' regimes.


Introduction
Metamaterials are synthetic materials that are designed to exhibit surprising properties that are rarely found in nature. Following the first practical demonstration of negative refraction [1] and the perfect lens [2][3][4], there has been an explosion of interest in metamaterials. One of the most unusual classes of metamaterial structures is the class of hyperbolic metamaterials, which currently enjoys significant interest in optics and plasmonics [5], as well as acoustics and solid mechanics [6][7][8][9][10][11][12]. Optical hyperbolic metamaterials are characterised by their having effective permittivity or permeability tensors with principle components that differ in sign. It is this feature that allows lenses created with hyperbolic metamaterials to break the diffraction limit. More generally, in the setting of time-harmonic fields, hyperbolic metamaterials are governed by effective hyperbolic partial differential equations (PDEs) as opposed to elliptic PDEs for classical wave-like systems.
In the present paper we describe a new class of dynamic multi-scale structures, which we refer to as parabolic metamaterials. This novel group of materials is characterised by having effective partial differential equation of the parabolic type and can exhibit a spatially localised unidirectional wave propagation. Similar effects have been observed in a different physical configuration, for flexural waves in periodically pinned Kirchhoff plates [13], but were not associated with the notion of parabolic metamaterials.
The concept of parabolic metamaterials is associated with the presence of Dirac points on the dispersion surfaces. Dirac points have been extensively studied in both the physics and mathematics literature (see [14,15], among others), motivated by the novel electronic transport properties of graphene [16]. Dirac cones have also recently been observed in phononic crystals by exploiting accidental degeneracies caused by appropriately tuning material and geometric parameters in order to induce degeneracies at the centre of the Brillouin zone [15]. In the present paper, the Dirac cones we exploit result from symmetries of the lattice system (so-called essential degeneracies). These Dirac cones are connected by relatively flat narrow regions on the dispersion surfaces, which we term Dirac bridges. These Dirac bridges possess resonances where the dispersion surfaces are locally parabolic, which give rise to highly localised unidirectional wave propagation.
We emphasise that the effects described in the present paper are dynamic in nature and the parabolic metamaterial behaviour is associated with resonances of the system. With this in mind, the novel method of high frequency homogenisation developed in a recent series of papers [17][18][19][20][21][22], will prove useful in the description and analysis of these parabolic metamaterials. Classical multi-scale asymptotic homogenisation approaches for microstructured media (see [23,24] among others) usually involve the study of a class of static model problems on an elementary cell and can be considered valid in the long-wave quasi-static regime. In contrast the approach of high frequency homogenisation considers perturbations away from known resonances and is therefore ideally suited to the study of dynamic parabolic wave propagation and the associated problems examined in the present paper. The approach yields effective partial differential equations with constant coefficients that govern the dynamic macroscopic behaviour in the vicinity of resonances. The methodology has proved very useful for a range of physical systems including continuous [20] and discrete elasticity [22], electromagnetism [21], platonics [19], as well as scalar problems [17] where the approach has been particularly useful in the study of hyperbolic dispersion surfaces, dynamic anisotropy, and star-shaped waveforms [25].
The concept of dynamic anisotropy and localisation in periodic lattices has been studied in several earlier papers [6][7][8][9][10][11][12]. In the context of the theoretical studies [9][10][11][12], hyperbolic metamaterials can be viewed as multi-scale structures that support a class of Bloch-Floquet waves and exhibit locally hyperbolic behaviour in the neighbourhood of certain resonances. In such cases, strong localisation can occur with wave propagation only permitted along directions associated with the principle curvatures of the hyperbolic surface. The paper [11] provides an important connection between the dispersive properties of Bloch-Floquet waves in an infinite lattice and the solutions of forced problems, which exhibit spatially localised waveforms of different types.
The present paper is developed as follows. The formulation of the problem and discussion of dispersion surfaces are given in section 2. The notions of parabolic and hyperbolic metamaterials are discussed in section 3, where we also discuss the method of high frequency homogenisation. The theoretical discussion is complemented with numerical computations illustrating the highly localised wave propagation associated with parabolic metamaterials. The paper is finalised in section 4 with some discussion and concluding remarks.

Formulation of the problem
In the present paper we study the in-plane motion of an elastic triangular lattice which, initially, might seem to be a conventional and well studied mechanical system. Conventionally such structures would be considered to be so-called "stretch dominated" [26], that is, the lattice remains statically determinate if the flexural rigidity of the links is neglected. Therefore, the bending effects in asymptotic models of thin elastic ligaments are often neglected for triangular lattices, both in statics and dynamics. However, here we show that this asymptotic approximation requires an adjustment in the dynamic regime and the effects of bending, generally neglected, play a fundamental role leading to a novel class of metamaterial behaviour. In the present paper we introduce the notion of a new type of elastic metamaterials referred to as parabolic metamaterials.
The problem is formulated using the same framework as established in [11,22]. The geometry is shown in figure 1, the elastic ligaments (denoted by grey links) are treated as massless Euler-Bernoulli beams, which allow for both flexural and extensional modes, and the junction points have both translational and rotational inertia. It is convenient, without loss of generality, to t 1 t 2 Figure 1: The uniform triangular lattice with the black balls denoting concentrated masses and the grey lines massless elastic beams. The primitive cell is shaded in red and the direct lattice vectors are shown, also in red.
normalise the system such that both the length of the ligaments and the mass of each junction is unity. Each node is labelled by the integer multi-index m ∈ Z 2 such that the spatial position of the m th node is x(m) = Tm, where T = [t 1 , t 2 ] is matrix with columns formed from the direct lattice vectors t 1 = [1, 0] T and t 2 = [1/2, √ 3/2] T . Each node is connected to its set of nearest neighbours, denoted by N = {p : x(m + p) − x(m) ≤ 1}. We emphasise that 0 ∈ N , i.e. each node is connected to itself. The inertial properties of each node are described by the matrix M = diag [1, 1, J] where J is the polar moment of inertia of each node. Finally, the elastic stiffness of the ligament connecting node m + p to node m is described by the matrix C(p). Using this notation, the equations of motion for time-harmonic waves propagating in the lattice can be conveniently written in the form where u(m) = [u 1 (m), u 2 (m), θ(m)] is the generalised displacement of each node: u 1 and u 2 are the usual Cartesian displacements and θ is the angle of in-plane rotation. The stiffness matrices C(p) can be deduced from those given in our earlier paper [11, Eq. 9] by taking the limit → 0. For the convenience of the reader, we also provide them in appendix A. We introduce the discrete Fourier transform and its inverse for the uniform triangular lattice under consideration. Applying the discrete Fourier transform to the equations of motion (1), we find whence the dispersion equation is immediately obtained The dispersion equation is a cubic polynomial in ω 2 and its roots can therefore be expressed in closed form; the expressions for the roots are somewhat tedious and are therefore omitted for brevity. The dispersion surfaces shown in figure 2 are a visual representation of the roots of the dispersion equation and show the frequency ω as a function of Bloch vector k. Figure 2 shows two cases: the first, in panel 2a, is for the case when the flexural modes of the elastic ligaments are neglected and only extensional modes are considered; the second, in panel 2b, shows the same configuration but with the bending modes of the links taken into account. The additional, predominantly flat band, shown in 2b is a consequence of the bending modes of the ligaments and represents a rotationally dominant mode. It should be emphasised that, in the former case where only extensional modes are considered, the shape of the dispersion surfaces is essentially determined by the geometry. Altering the stiffness of the elastic ligaments and inertia of the junctions would fundamentally only act to stretch the dispersion surfaces. This is in contrast to the latter case where the flexural deformation of the elastic ligaments is taken into account. In this non-dimensional form there are two free parameters which affect the shape and curvature of the dispersion surfaces: β, roughly speaking the ratio of bending stiffness to extensional stiffness of the ligaments, and J, which characterises the ratio of rotational inertia to translational inertia of the junctions. In terms of physical parameters, β = 2I/(S 2 ) and J =J/(m 2 ) where m andJ are the mass and polar moment of inertia of the lattice junctions respectively; I, S, and are the area moment of inertia, cross-sectional area, and length of the lattice ligaments respectively. These surfaces have several interesting features including the characteristic semi-infinite stop band associated with discrete structures, flat bands which can give rise to slow waves, and the presence of Dirac cones that have now become synonymous with hexagonal and triangular geometries. However our interest lies with a novel feature of the dispersion surfaces that, although related to Dirac cones, has received little attention. These features are only present in the case when the flexural effects of the elastic ligaments are included, as shown in figure 2b. In particular, it will be shown that the points indicated by B (±j) correspond to resonances and, at these points, the dispersion surfaces are locally parabolic leading to what we term parabolic metamaterial effects. These parabolic regions connect neighbouring Dirac points, hence we refer to them as Dirac bridges.

Constructing Dirac bridges
The presence of Dirac cones is generally associated with symmetries of the system through its geometry (see, for example, [14]) and therefore we would always expect regular triangular geometries to possess Dirac cones. The presence of parabolic Dirac bridges, in contrast, requires both the existence of Dirac cones in addition to the careful tuning of material parameters. Although one can obtain explicit expressions for the dispersion surfaces and, hence, their curvature in closed form, these expressions are often tedious and it can be difficult to extract any physical insight from them. Here, we employ the recently developed theory of high frequency homogenisation [17,18]. The multiple-scales approach, originally developed in [17], involves perturbations away from resonances and provides the effective continuous partial differential equation that governs the long scale behaviour of the medium. The rapid oscillations associated with the resonance and multiple scattering are encoded in the material parameters of the effective medium and reflected in the coefficients of the homogenised PDE. The approach has recently been extended to elastic lattices [22], which makes it an ideal tool to analyse Dirac bridges for the structures considered here.

High frequency homogenisation
We will now briefly summarise the essential features of the high frequency homogenisation methodology. For further details and a complete description of the approach, the reader is referred to [22] and reference therein.
We begin by introducing two scales: the short-scale discrete variable m, and the long-scale continuous variable η = εTm. It is assumed that the small parameter 0 < ε 1 characterises the short-scale of the lattice. For example, if the lattice is formed by an N × N grid of particles, where N 1, then ε = 1/N . The displacement is then considered as a function of two independent vector-valued variables: u = u(m, η). Expanding the Fourier-transformed equations of motion (2) in the continuous variable η for small ε we find and ∇ acts on the continuous long-scale variable η. Next, we take the following Ansätze for the Fourier-transformed displacement and frequency-squared which yields a hierarchy of algebraic problems in ascending orders of ε. For most of what follows, it is necessary to consider only the first three problems, where S i = ω 2 i M − σ i . Considering the leading order problem (7), we observe that S 0 is independent of the long-scale variable η and its derivatives. Therefore fixing k, which is equivalent to setting the phase-shift across the elementary cell, reduces (7) to an eigenvalue problem which, for simple eigenvalues, admits solutions of the form where we normalise such that |U F 0 | = 1. Indeed, for the present configuration, it will be sufficient to only consider non-degenerate eigenvalues. Examining the first order problem (8), we find that it is solvable if and only if (see [22]) where the dagger denotes the Hermitian Transpose. For the cases of non-degenerate eigenvalues considered in the present paper we find that U F 0 (k) † σ 1 (ω, k, ∇)U F 0 (k) = 0 which implies ω 2 1 = 0 and the first order problem admits the solution where the superscript plus sign denotes the Moore-Penrose pseudoinverse, ψ is an arbitrary vector, and I is the identity matrix of appropriate dimension.
Moving to the second order problem, applying the solvability condition, and using the fact that I − S + 0 S 0 is the orthogonal projector onto the kernel of S 0 , we obtain the second order partial differential equation for the leading order envelope function φ 0 (η) and the correction to the frequency where we suppress the arguments of the operators. Equation (11) represents the homogenised partial differential equation that governs the long-scale response of the medium. We emphasise that the envelope function in (11) is a function of the long-scale continuous variable η only. The short-scale behaviour is encapsulated by the bracketed operator in (11); this operator corresponds to a second order partial differential equation with constant coefficients and will be used in the following section to construct the parabolic metamaterial.

Hyperbolic, elliptic, and parabolic metamaterials
We now return to the points B (±j) , which lie on the boundary of the Brillouin zone. By considering the symmetries of the lattice, we can restrict our attention to the irreducible Brillouin zone formed by a triangle with vertices at k = [0, 0] T , [0, 2π/ √ 3] T , and [2π/3, 2π/ √ 3] T . Only B (1) lies within the closure of the first Brillouin zone; the remaining points can be obtained by symmetry operations. Therefore, the following will be restricted to analysis of the point B (1) . We are primarily interested in the behaviour of the upper two surfaces at this point as these exhibit the desired locally parabolic behaviour. In particular, on these surfaces and for certain values of β and J, the point B (1) lies at the centre of the Dirac bridge which connects adjacent Dirac points.
Using the approach outlined in §3.1, we now proceed to derive the effective partial differential equations that govern the response of the lattice in the neighbourhood of B (1) . The leading order problem (7) has three eigenvalues . Therefore, the first eigenvalue corresponds to the lower dispersion surface with the latter two associated with the upper two surfaces. Proceeding to higher order, we find that the second and third modes yield effective partial differential equations of the form where p = 1, 2, 3 and enumerates the modes. The bracketed differential operator in (12) is the explicit form of the same bracketed operator appearing in (11). The matrices A (p) are all diagonal and have the following forms where α 1 = (18J − 5)β + J and α 2 = (6J − 5)β + 3J. Upon examination of (13)- (15), it can be shown that the eigenvalues of A (p) can have different signs depending on the values of β and J, leading to effective partial differential equations. As remarked above, the first mode is typically associated with the lowest frequency band and although the desired locally parabolic behaviour is present, it will occur at frequencies where several other modes exist making the parabolic behaviour difficult to isolate. We therefore focus on the higher two modes for p = 2, 3. In order the classify the different behaviours expected from the homogenised partial differential equations, we introduce the following definition, that is consistent with the classical classification of linear second order partial differential equations.
Definition 3.1. Defining the discriminant of the homogenised partial differential equation (12) as ∆ (p) = det A (p) , the partial differential equation is then said to be Figure 3 shows a plot of sgn(∆ (p) ) for a range of (J, β) values, illustrating the different regimes. The grey (white) regions indicate that the effective partial differential equation is hyperbolic (elliptic). These regions are demarcated by solid black curves where the effective partial differential equation is parabolic. By tuning the material parameters of the triangular lattice, we can control the effective dynamic behaviour of the medium. For monochromatic waves, elliptic equations represent 'normal' wave propagation (e.g. the Helmholtz equation for acoustic waves), hyperbolic equations result in highly localised waves along a discrete set of principle directions (see, for example, [11]), and parabolic equations are associated with strongly localised beams where waves are only permitted to propagate along a single principle direction.

Parabolic profiles
Parabolic partial differential equations are characterised by vanishing discriminants, which implies that one of the eigenvalues of A (p) is zero. These eigenvalues are quadratic polynomials in β and therefore we can obtain their roots explicitly. Since β and J represent stiffness and mass respectively, we only consider positive roots. In particular, for the second eigenmode the effective PDE is parabolic when β (2) 11 = 7 18 and β where β (2) jj indicates that the coefficient of ∂ 2 /∂η 2 j vanishes for this value of β. Both β and J correspond to physical quantities (roughly speaking the ratio of flexural to compressional stiffness of the elastic ligaments and the ratio of translational to rotational inertia respectively) and therefore should take positive finite values. The branch of the square root is chosen accordingly.
For the third eigenmode we have similar results In order to satisfy the condition that β > 0, the second root β (11,2) is restricted to 0 < J < 29/42. These values are indicated by the solid black curves on figure 3 and demark the boundary between the regions where the effective PDE is either hyperbolic or elliptic. It is interesting to note that β (2) 11 is independent of J which suggests that, for this class of lattice, one can obtain parabolic behaviour by simply tuning the stiffness of the elastic ligaments without altering the inertial properties thus simplifying design and implementation of parabolic metamaterials.
Special attention is required at the poles of the tensors A (i) , that is, when α 1 = 0 and/or α 2 = 0. When α 1 = 0, the first and third eigenmodes coincide whilst the first and second eigenmodes coincide when α 2 = 0. All three modes coincide when α 1 = α 2 = 0 (J = 24/5). In such cases of degenerate eigenvalues, the leading order solution is a linear combination of the eigensolutions for each eigenvalue where U F (0,i) (k) is the eigenvector corresponding to ω 2 (0,i) , φ (0,i) (η) is the associated longscale envelope function, and D is the set of indices enumerating the degenerate modes. The order of degeneracy is then D . The solvability conditions for the first order problem are then  which together yield a coupled system of D first order partial differential equations for the longscale envelope functions and the first order correction to the frequency. In some cases the system may be explicitly decoupled. For example, when α 1 = 0 and the first and third modes coincide, we find that the coupled system reduces to a single uncoupled parabolic partial differential equation governing the envelope functions for the first and third eigenmodes This case (α 1 = 0) is associated with the lower solid curve in figure 3a. Similarly, when α 2 = 0 the first and second modes coincide and we find that the associated envelope functions are governed by another parabolic equation This case corresponds to the the middle solid curve in figure 3b.

Localisation and wave beaming in parabolic metamaterials
For the purposes of illustration, we now consider the forced problem where the infinite lattice is excited at the origin by a time-harmonic force. In this case, the equations of motion are where F is the forcing amplitude and δ is the discrete Dirac delta. The response of the lattice can be characterised by the Green's function, which we express in the form of a double Fourier integral Here we are concerned with frequencies that lie within the stop band for which det σ(ω, k) = 0 and the integrand is well behaved. In general, equations of the form (18), cannot be evaluated in closed form, although for lower dimensions or simple configurations some progress can be made (see, for example, [27,28]). However, the Green's function (18) is amenable to numerical evaluation, particularly outside the pass band. This is the approach that we take here where we directly evaluate the double integral (18) using iterative quadrature. Figure 4 shows two configurations where the material parameters have been tuned to achieve an effective parabolic metamaterial. A point force generates "beam-like" fields, where the energy is highly localised within a narrow beam; this effect is due to the parabolic nature of the lattice. The lattice structure is forced by a point load in the bottom left hand corner and a single unidirectional beam is emitted. The beam consists of plane waves propagating along the direction defined by the effective parabolic PDE and modulated by a Gaussian-like beam. In each case we show only the dominant component of the displacement; the omitted components are small at the forcing point and decay rapidly away from the loading. Figure 4a corresponds to the second mode where we choose β = β (2) 11 = 7/18. As remarked earlier, this configuration has the interesting feature that it is parabolic regardless of the chosen value of J, meaning that only the stiffness of the lattice links needs to be tuned and the mass of the lattice nodes can be chosen arbitrarily. This feature makes the practical implementation of parabolic metamaterials more straightforward. The second panel, figure 4b, corresponds to β = β It is interesting to note, however, that the resulting wave pattern remains essentially unchanged regardless of the chosen orientation of forcing; the highly localised unidirectional beam is still present, but the dominant component of displacement may change.
We remark that the numerical evaluation of (18), although straightforward, is computationally intensive. Indeed, the computations shown in 4 consists of an 80 × 80 grid of nodes at each of which three double integrals corresponding to each component of (18) is evaluated. The computation was parallelised and implemented in MATLAB. Although we have not carried out an explicit analysis of the computational cost, it is illuminating to consider that each of the two computations shown in figure 4 required approximately 24 hours to complete using 8 dedicated physical cores running at 2.5GHz. This emphasises the usefulness of the elegant closed form asymptotic solutions presented in this paper.

Concluding remarks
In conclusion we emphasise that, for an elastic triangular lattice, the flexural rigidities of the elastic ligaments can be neglected and the longitudinal stiffness plays the dominant role in the determination of the effective material properties. Such an approximation is fully applicable to static problems. However, we have demonstrated that the flexural rigidities of the elastic beams within a triangular lattice play an important role in higher-frequency regime linked to description of a new class of metamaterials, which we refer to as 'parabolic metamaterials'.
The notion of 'hyperbolic metamaterials' is well-established (see, for example, [5]) and can be associated with important features of dynamic anisotropy of periodic multi-scale structured media (see [6][7][8][9][10][11][12], among others). However the new class of parabolic metamaterials, identified for the first time in the present paper, is characterised by parabolic effective equations on the long scale, and hence a forced input would produce a highly localised unidirectional Gaussianlike beam in the appropriate dynamic regime. The notion of a 'Dirac bridge' is introduced to describe such dynamic regimes. These bridges are sections of the dispersion surfaces that connect adjacent Dirac cones. At the centre of each Dirac bridge is a stationary point associated with a resonance of the system; the dispersion surfaces are locally parabolic in the neighbourhood of these resonances which are associated with the described parabolic behaviour.
The high frequency homogenisation method employed in the present paper leads to explicit formulae for the coefficients of the envelope equations derived here. It has also been demonstrated how mechanical parameters of the elastic ligaments, such as the moment of inertia of the cross section as well as the ratio of elastic stiffnesses, influence the macroscopic of the metamaterial. It is interesting to note that, for the triangular lattice, these Dirac bridges only appear when the flexural interaction of the thin elastic ligaments is considered; this is despite the common convention [26] to ignore bending moments for so-called 'stretch-dominated structures'. Whilst this may be valid in the classical long wave regime, we have demonstrated here that such interaction play a vital role in the dynamic regime and can lead to novel and exciting effects.
Finally, this analytical study opens perspectives for dynamic analysis of defects and exponential boundary layers at the edges of large structured clusters.