Dynamic Green's functions in discrete flexural systems

The paper presents an analysis of the dynamic behaviour of discrete flexural systems composed of Euler--Bernoulli beams. The canonical object of study is the discrete Green's function, from which information regarding the dynamic response of the lattice under point loading by forces and moments can be obtained. Special attention is devoted to the interaction between flexural and torsional waves in a square lattice of Euler--Bernoulli beams, which is shown to yield a range of novel effects, including extreme dynamic anisotropy, non-reciprocity, wave-guiding, filtering, and the ability to create localised defect modes, all without the need for additional resonant elements or interfaces. The analytical study is complimented by numerical computations and finite element simulations, both of which are used to illustrate the effects predicted. A general algorithm is provided for constructing Green's functions as well as defect modes. This algorithm allows the tuning of the lattice to produce pass bands, band gaps, resonant modes, wave-guides, and defect modes, over any desired frequency range.


Introduction
In recent years, the study of lattice dynamics has undergone something of a renaissance due, in part, to the counter-intuitive effects associated with metamaterials, such as perfect lenses [1 ], negative refractive index materials [2 ], and invisibility cloaks [3 ]. The majority of scholarly work in this area has focused on photonic, plasmonic, and acoustic metamaterials, whilst the science of mechanical metamaterials attracts less attention. Nevertheless, elastic metamaterials have found applications a wide array of settings, from cloaking to energy dissipation in engineering structures and seismic protection, among many others [4 -11 ].
Much of the work on metamaterials is underpinned by the study of wave propagation in periodic media, which has been studied in various forms for centuries and remains an active area of research today [12 -14 ]. Dynamic vector systems corresponding to multiscale mechanical materials offer unique opportunities to study the flexural and rotational displacements that occur in elastic lattice systems [15 -22 ]. Devices that combine, for example, beams and plates with additional resonating structures have been used to produce novel effects associated with bending and rotation, leading to phononic properties not typically found in optical or plasmonic materials, such as uni-directional wave modes, allowing the transmission and localisation of energy [18 , 23 -27 ]. The Green's function is the canonical object of study for many problems associated with wave propagation in structured solids as it contains all the fundamental information (such as the dispersive properties) corresponding to the dynamic response of the system. There is a substantial amount of literature on Green's functions, particularly for scalar systems, such as the wide variety of topologies and wave propagation modes (free and forced) studied on mass-spring lattices, [19 , 28 -30 ]. Green's functions afford the means to control wave propagation and analyse defects that occur in lattice systems [16 , 31 ]. A method for the localisation of energy in scalar lattices was provided in [32 ], where the Green's function was used to customise evanescent defect modes that coincide with localised band gap modes.
Considering lattices of beams, various studies have explored systems of Rayleigh beams, including different topologies, the effects of rotational inertia, and the propagation of both in-plane and out-of-plane vibrations [33 -36 ]. Forcing in the form of applied moments is uncommon in flexural systems, but has been shown to produce interesting wave patterns with varying degrees of anisotropy [35 ].
Lattices of Euler-Bernoulli beams have been the subject of studies focusing on many different facets, including their topology and propagation modes, and have been combined with additional resonators to produce interesting phenomena that can be used as wave guides [19 , 20 , 36 -38 ]. One such example [34 ] demonstrated how the interface between a half-plane of Euler-Bernoulli beams, and a half-plane of Rayleigh beams can produce negative refraction for incoming lattice waves. Another example [18 ] combined a periodic array of Euler-Bernoulli beams on a plate with gyroscopic resonators, leading to 'chiral flexural waves' and one way edge waves at the chiral interface. Unlike lattices of Rayleigh beams, rotational inertia is usually neglected on lattices of Euler-Bernoulli beams as standard. In this study, we explore the effect of rotational inertia and torsional stiffness in lattices of Euler-Bernoulli beams with remarkable results. We emphasise that torsional interactions are often neglected when studying flexural systems insofaras the two are treated independently. However, we show that properly accounting for the interaction between these flexural and torsional effects leads to unprecedented levels of control over the dispersive properties of the lattice. Alongside these observations, we also compare the different waveforms generated in the lattices from applied forcing, versus applied moments and show that with applied moments, the lattice can act as a wave guide without the need for additional resonators or interfaces.
The structure of the paper is as follows: firstly, in § 2 we construct a class of Green's functions for a 2D square lattice of Euler-Bernoulli beams, where the nodes have mass and rotational inertia and the beams possess torsional stiffness. In § 2.4, we use applied forces and moments to induce dynamic anisotropy, including unusual asymmetric propagative (nonreciprocal) modes. In § 2.5 we show that altering the magnitude of the rotational inertia and torsional stiffness provides a range of possible (and significantly different) dispersion diagrams. In § 3 we study the related problem of a 1D chain of Euler-Bernoulli beams with nodes that possess rotational inertia, and show that the absence of torsional interactions allows us to evaluate the full class of Green's functions for the chain in closed form. Lastly, in § 4 we generalise the approach used in the previous sections to provide a method for constructing Green's functions on discrete lattices of any dimension, with lattice connections that take any form, and then go on to derive custom localised defect modes.
In this study, the control we demonstrate over the direction of wave propagation lays the foundation for designing metamaterials of Euler-Bernoulli beams that have tailor-engineered properties.

A square lattice of Euler-Bernoulli beams
In this section we consider an infinite square lattice of thin Euler-Bernoulli beams with junctions possessing both mass and rotational inertia. The beams have both flexural rigidity and torsional stiffness, meaning that the torsional rotation and flexural displacement at the junctions are coupled. In particular, whilst the flexural displacement and rotation within a beam are coupled, the torsional rotation is decoupled within the beam itself. At the junctions, however, the flexural rotation in a beam will induce a torsional rotation in the perpendicular beams, coupling the two classes of waves. This coupling of flexural and torsional interactions is often overlooked in flexural systems but, as will be shown later ( § 2.4), properly accounting for these interactions leads to novel and interesting effects. In Figure 1, the eigenmodes of the lattice are illustrated using a finite element model of the lattice unit cell. We construct the lattice in the xy-plane and therefore define the out-of-plane translational motion along the z-axis. For convenience, we choose the mass of the nodes, the length of the beams and the flexural rigidity as natural units. Each lattice node is labelled by the double-index (m, n) ∈ Z 2 , such that 0 denotes the node located at (x, y) = (0, 0). We introduce the generalised displacement vector u (m,n) = [w(m, n), θ x (m, n), θ y (m, n)] T to describe the displacement and rotation of each node. The first component of u (m,n) describes the translation w(m, n) of the (m, n) th node, while the second and third components, θ x (m, n) and θ y (m, n), characterise the rotational displacement of the nodes around the respective coordinate axes, as indicated in Figure 2b. Figure 2: (a) The Euler-Bernoulli beam lattice geometry, with (b) the corresponding coordinate axes and displacements.
The four beams connected at the node (m, n) are denoted by X + , X − , Y + and Y − corresponding to the axes on which the beams lie. Hence, the beam spanning m ≤ x ≤ m+1 is denoted X + and so forth, as shown in Figure 2a. The out-of-plane flexural deformation w(x, y) of the lattice links is governed by the fourth order differential equation for Euler-Bernoulli beams, ∂ 4 w ∂x 4 = 0, on X ± and ∂ 4 w ∂y 4 = 0, on Y ± . (2.1) The first-order spatial derivatives w ,x (x, y), w ,y (x, y), give the angles of flexural rotation, about the x-and y-axes respectively, associated with the flexural deformation; here subscript letters proceeded by commas indicate partial derivatives with respect to the relevant spatial variables. We define τ x (x, y) as the torsion angle that the X ± beams experiences about the x-axis, and likewise for τ y (x, y) about the y-axis. The torsion angles each satisfy To construct the equation of motion, we first consider the forces and bending moments that the X ± beams apply to the node (m, n). From equations (2.1) and (2.2), we see that the flexural deformation and torsion angles are cubic and linear polynomials respectively. The coefficients of these polynomials are found by imposing boundary conditions at the ends of the X ± beams; we use the following constants, x .
Careful attention is required when considering the directions of the flexural moments and torsional moments around the in-plane coordinate axes, shown in Figure 2b. We stipulate that positive Θ x refers to an anticlockwise angle around the x-axis, and likewise for the yaxis. Using the boundary conditions, we arrive at expressions for the flexural displacement of the X + beam and the torsion angle of the X + beam We also find the equivalent expressions for the X − beam, w(x, n) X − and τ x (x, n) X − . The shear force, the flexural bending moment and the torsional moment induced in the X ± beams are respectively, where C in M tors is the non-dimensionalised torsional stiffness coefficient and the negative sign corresponds to the direction of the moment about the x-axis [39 ]. The forces and moments that a beam applies to the node (m, n) can be expressed in terms of the generalised forcing vector F (m,n) = [F (m, n), Φ x (m, n), Φ y (m, n)] T ; where the moments Φ x (m, n) and Φ y (m, n) describe the total bending moment about their respective axes at node (m, n), and so each contains both flexural moments M f lex and torsion M tors . As such, the forces and moments that the X ± beams apply to the node (m, n) are (2.6) and we label the matrices for each beam X + A and X + B , X − A and X − B , as indicated in equation (2.6).
Following the same method, we apply boundary conditions to the ends of the Y ± beams; since flexion in the X ± beams induces torsion in the Y ± beams, this is reflected in the continuity condition at (m, n) such that w ,x (m, n) = Θ (0) y = τ y (m, n) and w ,y (m, n) = −Θ (0) x = τ x (m, n). Using the boundary conditions, the polynomials for the flexural displacement w(m, y) Y ± and torsion angles τ y (m, y) Y ± associated with the Y ± beams are found. The forces and bending moments that the Y ± beams apply to the node (m, n) are then found using the y-axis equivalent expressions to equation (2.5). In parallel with equation (2.6), we express the forces and moments from the Y ± beams as vectors, and we label the corresponding matrices Combining the forces and moments in both the x-and y-directions with Newton's second law for a time-harmonic system, the equation of motion of the (m, n) th node is then where u (m,n) is again the generalised displacement vector of the node (m, n) and the 3 × 3 matrix M = diag[1, µ, µ] describes the inertial properties of the lattice. The first component of M corresponds to the (unit) mass of the junction associated with the translational inertia, whilst the second and third components correspond to the rotational inertia associated with the flexural and torsional deformations, which we denote µ. We apply the discrete Fourier transformation T is the generalised displacement in reciprocal space in terms of the spectral parameters k 1 and k 2 . Following the Fourier transformation, we arrive at the lattice's equation of motion in reciprocal space, As such, the three components of u F (k 1 , k 2 ) are defined as the Fourier transformed real space displacements w(m, n), θ x (m, n) and θ y (m, n), respectively.

The dispersion equation
The dispersion equation σ(ω, k 1 , k 2 ) = 0 arises from the solvability condition of equation (2.9), where σ(ω, k 1 , k 2 ) = 144 sin 2 (k 1 )ζ(k 1 , k 2 ) + 144 sin 2 (k 2 )ζ(k 2 , k 1 ) and we define the repeated function ζ(P, Q) = −2C cos(P ) + 2C + 4 cos(Q) − µω 2 + 8 . As one would expect for a symmetric square lattice, the dispersion equation is symmetric in the spectral parameters k 1 and k 2 . We also see that the dispersion equation is cubic in ω 2 and therefore, closed form solutions can be found but are omitted here for brevity. In addition to the frequency and spectral parameters, σ(ω, k 1 , k 2 ) is dependent on the material constants µ and C, which describe the non-dimensionalised rotational inertia and torsional stiffness, respectively. In § 2.5 we investigate the effects of changing µ and C on the dispersion equation and show that it is possible to produce completely different dispersion diagrams, with and without the existence of a finite band gap. This control over the lattice's propagating frequencies (with their preferential directions and subsequent anisotropy) is a particularly interesting feature that provides applications for metamaterials designed to control the propagation of waves in structures.
In Figure 3, we plot the solutions to σ(ω, k 1 , k 2 ) = 0 for the illustrative values of µ = 0.01 and C = 0.1. It is shown that these values of rotational inertia and torsional stiffness produce a dispersion diagram with two distinct pass bands and two band gaps, one finite and also the semi-infinite band gap associated with discrete systems. Given that the closed form solutions to the dispersion equation are in-hand, we can determine the exact limits of the band gaps for any desired µ and C.

Constructing the Green's function
A notable feature of the flexural lattice is the ability to induce translational displacements, rotational moments and combinations thereof. To construct the Green's function, we consider the application of time-harmonic unit forces and unit moments to the central node of the lattice through the forcing vector f ∈ C 3 . In parallel with the components of F (m,n) , the forcing vector f = [f w , f θx , f θy ] T has components corresponding to the application of translational force f w out-of-plane (along the z-axis) and two moments f θx and f θy about the x-and y-axes respectively. For any chosen frequency, there are seven different forcing configurations. These modes are generated by the permutations of f w , f θx , f θy ∈ {0, −1} with the exclusion of f w = f θx = f θy = 0 corresponding to the absence of any applied forces or moments. We express the Green's function in reciprocal space in the form The inverse of the discrete Fourier transform is then applied to equation (2.12) to find the Green's function in real space, u (m,n) . The resulting components of u (m,n) are written in terms of Fourier integrals and, whilst useful, are cumbersome. Therefore, in the interests of brevity and clarity of presentation, we will restrict ourselves to consideration of the flexural displacement component W F (k 1 , k 2 ) of the Green's function u F (k 1 , k 2 ) for an arbitrary applied forcing. In this case, the spectral form of the flexural displacement is, where the functions σ and ζ are given in equations (2.10) and (2.11) respectively. The difference of sign for the coefficients of f θx and f θy in equation (2.13) is a consequence of the sign convention adopted for shear forces and moments in equation (2.5). As expected for a symmetric lattice, the coefficient of the linear force f w is symmetric in the two spectral parameters. We also note that the denominator of equation (2.13) coincides with the dispersion equation, which arises from the inverted operator matrix in equation (2.12).
The inverse Fourier transformation is applied to find the displacement in real space as follows, (2.14) Although explicit closed form representations do not exist in general, the integral form of the Green's function is amenable to numerical evaluation. In the following sections, we use numerical evaluation and finite element models to explore the effects of changing the forcing vector, and also the frequency with which the force is applied.

Anisotropic behaviour of localised modes
For a lattice with µ = 0.01 and C = 0.1, consistent with the dispersion diagram in Figure 3, the lower boundary of the finite band gap is ω = 4 √ 6. Using equation (2.14), the flexural displacement band gap Green's functions of the lattice nodes, w(m, n) is evaluated numerically for three examples of the seven forcing options, with each localised mode demonstrating a unique shape. Across each example in Figure 4, we chose the frequency of excitation to be ω = 9.8, such that it lies in the band gap, close to the band-edge; this demonstrates the rapid rate of decay and high degree of localisation associated with band gap Green's function in this regime.
Firstly, we consider out-of-plane forcing f = [−1, 0, 0] T and plot the flexural displacement in Figures 4a and 4b. In Figure 4a, the decay envelope of the displacement is clear, while Figure 4b displays the periodic oscillation of the nodes and the symmetry between the x-and y-directions. In mechanical lattices the preferential direction of wave propagation is most commonly seen along the principle axes of the lattice. Here, the localised mode is radially symmetric leading to isotropic behaviour which can be understood from the complex solutions (k 1 , k 2 ) ∈ C 2 to the dispersion equation, leading to a uniform decay rate in all directions.
When f = [0, −1, 0] T , a clockwise moment is induced about the x-axis at the 0 node; for this applied moment, the flexural displacement of the lattice nodes has been plotted in Figures 4c and 4d. While the displacement is localised radially, the displacement is not localised to the y-axis as one might expect for this type of forcing, this is again due to the complex solutions of the dispersion equation in this regime, meaning the mode has no preferential direction of travel. This is in particular comparison to Figure 5b where, for the same applied moment at a pass band frequency, the displacement field of the propagating mode is highly localised along the y-axis. In Figure 4d, we have indicated the dotted line (m, 0) (the x-axis) where the lattice nodes do not have out-of-plane displacement but do experience rotational displacement, as expected with this type of forcing. We note that a moment applied about the y-axis f = [0, 0, −1] T produces the same out-of-plane displacement in the lattice under a π/2 rotation about the z-axis. With f = [0, −1, −1] T , two unit moments are applied clockwise around the x-and y-axes at 0 simultaneously and the resulting flexural displacement is plotted in Figures 4e and 4f. As with the single moment forcing above, the displacement field generated by the double moments is radially localised but not confined to a single direction; in this case, the line of zero flexural displacement lies on π/4, as has been indicated by the dotted line.
In addition to evaluating the flexural displacement w(m, n), we can use the inverse Fourier transformation to evaluate the θ x (m, n) and θ y (m, n) components of the displacement vector. These components give a quantitative measure for the magnitude of the rota-tion each node experiences around the relevant axis. In the interests of clarity, the plots of the Green's functions for rotational displacements are omitted here because they do not lend themselves to convenient graphical representation. However, it is important to note that these components contribute to the overall lattice response, can be evaluated as described above, and play an important role in the construction of defect modes (cf. § 4.2).

Extreme dynamic anisotropy and non-reciprocity in the pass band
The dispersive behaviour of lattice systems and the associated extreme anisotropy has been well documented and leads to interesting effects such as highly localised waveforms, unidirectional waves, one-way edge modes, and DASERs [17 -19 , 25 , 35 ]. However these effects are often highly narrowband and usually require careful tuning of the material and geometric parameters through, for instance, the inclusion of resonant elements. In contrast, here we show how localised waveforms and uni-directional waves are easily achieved in broad frequency regimes as a result of incorporating the additional degrees of freedom associated with the interaction between flexural and torsional motion at the nodes. Furthermore, the lattice also supports the highly unusual phenomena of non-reciprocity, despite the lattice being symmetric.
In parallel to the numerical results derived from the analytical Green's function, we develop a novel finite element model to illustrate a range of interesting phenomena in the infinite lattice. Using COMSOL Multiphysics ® , we build a lattice model of 201 × 201 Euler-Bernoulli beams and impose frequency-independent damping, in the form of complex stiffness values, on the beams in the vicinity of the boundary of the computational window in order to simulate an infinite lattice and minimise any artificial reflections. In this section, we set the values of µ = 0.01 and C = 0.1 for the rotational inertia and torsional stiffness respectively across all figures. Figure 5 illustrates the use of the forcing vector to control the propagation of waves through the lattice for pass band frequencies. In particular, we demonstrate that either unidirectional or asymmetric waves can be achieved by choosing the forcing vector appropriately. Firstly, in Figure 5a we choose f = [−1, 0, 0] T corresponding to an out-of-plane force of frequency ω = 5. The corresponding slowness contour is provided in Figure 5f; as expected, the principle directions of anisotropy coincide with the normals to the slowness contours, which themselves define the directions of maximal group velocity. We note the stark contrast between Figures 5a and 5b and emphasise that the only difference between the two figures is in the choice of forcing; the frequency of excitation is identical in both figures. In particular, Figure 5b corresponds to f = [0, −1, 0] T , which is associated with the application of a clockwise moment about the x-axis. The symmetric slowness contour indicates that a wave should propagate along both of the principle axes of the lattice equally, however when the lattice is forced by a moment about the x-axis, zero flexural displacement is induced in the x-axis, and a uni-directional waveform is induced which propagates in the y-direction only. Likewise, a uni-directional wave along the x-axis can be induced using the forcing vector f = [0, 0, −1] T . In this way, the forcing vector can be used to select the direction of propagation, similar to the mechanism identified in [19 ] for in-plane elastic systems.
The forcing vector also allows the application of combined forces and moments, producing modes with rare asymmetric anisotropy which can be used for applications in wave control and energy harvesting. In Figure 5c, the lattice was subjected to out-of-plane forcing with a simultaneous moment applied about the x-axis, such that f = [−1, −1, 0] T , at a frequency of ω = 5. The resulting wave propagates along the principle axes as expected, however with remarkably higher amplitude of displacement along the positive y-axis. The combination of the applied translational force and rotational moment causes the Y + beam to displace more than the Y − beam, which induces the asymmetry. Such non-reciprocity is rarely observed and usually requires PT-symmetry breaking but can sometimes be induced by creating asymmetric eigenmodes as in, for example, [27 ].
It is of note that while these effects have been demonstrated for the frequency ω = 5, these displacements can be achieved in a broad frequency regime. As can be seen in Figure 5f, there is an extended interval of frequencies for which the slowness contours all have the same quasi-rectalinear shape. The waveforms at these frequencies all display similar localisation as Figure 5a, and as such the same anisotropy can be induced through the choice of forcing vector; we further show in § 2.5 that the associated shape of the slowness contours is stable over a wide range of values of µ and C. This is especially advantageous for the implementation of these effects in practical devices; in previous studies, generating uni-directional waves has required exact values of the material parameters and forcing frequency [35 ] and also has been associated with so-called 'parabolic metamaterials' in narrow frequency regimes around Dirac cones [17 ]. In contrast we have shown that, by fully taking into account the flexural and torsional interactions present in the lattice, it is possible to achieve similar effects without requiring highly precise tuning of material parameters. This is particularly important for the fabrication of such devices where a degree of tolerance in the manufacturing process is required.
In Figure 5, we also demonstrate that it is possible to control the dynamic anisotropy and shape of the localised waveforms propagating through lattice by altering the forcing frequency. In Figure 5d, we subject the lattice to an out-of-plane force of frequency ω = 8, producing a wave that travels equally along the diagonals of the lattice consistent with the principal directions predicted from the slowness contours in Figure 5f. Comparing Figures 5a  and 5d, it is observed that the principle directions of the two localised waveforms are rotated by π/4; the only difference between these two figures is the chosen forcing frequency.
A further example, illustrating non-reciprocity, is shown in Figure 5e, which gives the displacement under f = [−1, −1, −1] T , inducing simultaneous out-of-plane forcing and double moments applied clockwise around both the x-and y-axes, at frequency ω = 8. While there is a small displacement propagating at π/4 to the principle axes, this is dominated by a much larger displacement propagating along the 3π/4 line. We also remark that the displacement along 3π/4 is larger than along 7π/4; although this is difficult to discern from Figure 5e, it can be verified by direct evaluation of equation (2.14).

The effects of altering the rotational inertia and the torsional stiffness
We have shown above that the lattice provides significant control over the propagation of waves from the variety of forcing options alone, with rare examples of anisotropy. It is now shown that the dependence of the dispersion equation (2.10) on the rotational inertia and torsional stiffness, not only provides extensive control over the location (and existence) of band gaps for the lattice, but also the shape of the dispersion surfaces and in turn the preferential direction of wave propagation. The different dispersion diagrams can be roughly classified into three shapes. When µ is very small (around 0.01 or below) and the torsional stiffness is allowed to range from very small to very large values, the dispersion diagrams have the same overall configuration as Figure 6a, most often with a finite band gap. Generally the boundary of the acoustic band (the lower dispersion surface) does not change, although the shape of the surface does alter, in turn producing slowness contours with very different shapes. When µ is 1 or larger, we observe a conical shape which is more commonly associated with beam lattices, such as in Figure 6b. we also see that the boundary of the semiinfinite band gap usually occurs as much lower frequencies. For intermediate values of µ, the dispersion diagram can have a shape similar to Figure 6a, although a finite band gap is seen less often, as is the case in Figure 6d. Alternatively, we can see an entirely different shape, such as Figure 6c, depending on the torsional stiffness. Note that Figures 6c and 6d were produced using the same value of µ, but different C. This shows how amenable the dispersive properties of the lattice are to manipulation, indeed altering the values of µ and C allows us to tailor the dispersion surfaces to desired choices of band gaps and slowness contours.

A chain of Euler-Bernoulli beams
In this section, we construct a new class of Green's functions for an infinite one-dimensional chain of thin Euler-Bernoulli beams with connecting junctions that possess both mass and rotational inertia. Unlike the square lattice studied in § 2, there are no torsional interactions to account for in the 1D chain. This simplification, along with having only one Fourier parameter means that the Green's functions can be evaluated in closed form, as will be shown. In a right-handed coordinate system, we construct the chain along the y-axis, and define the translational motion along the z-axis, as shown in Figure 7. While there are no torsional interactions, each junction in the chain has two degrees of freedom, corresponding to the translational displacement w and rotational displacement (about the x-axis) θ x . As before, we introduce the generalised displacement vector u n = [w(n), θ x (n)] T of the n th junction to describe the translational and rotational displacements of the nodes, which are coupled through the flexural deformations of the beams. The mass of the nodes, flexural rigidity and length of the beams are once again chosen as natural units. To construct the equation of motion for the n th unit cell, we evaluate the forces acting on the node at y = n from the beam on n ≤ y ≤ n + 1 which we denote Y + , using the same method as previously.
The flexural deformation of the Euler-Bernoulli beams, w(y) is governed by the same fourth order differential equation used in § 2, equation (2.1). As before, w(y) is a cubic polynomial in y, which is found using boundary conditions for the ends of the Y + beam, where W i are the translational displacements and Θ i are the angles of rotation of the junctions. The positive flexural rotations of each junction are defined anticlockwise about the x-axis, shown in Figure 7. We use the same equations for the shear force F (y) and bending moment from the flexural rotations M f lex (y) as in § 2, equation (2.5). We introduce the generalised forcing vector F n = [F (n), M f lex (n)] T and as such, the force and the bending moment at y = n associated with the Y + beam can be expressed as where the matrices A = −12 −6 −6 −4 and B = 12 −6 6 −2 encapsulate the flexural rigidity of the beams connecting the n th node to the (n + 1) th node. Considering Figure 7, it can be seen that the translational forces at the y = n + 1, n − 1 nodes have opposite sign, while the bending moments have same direction. Hence we use the matrix R = diag[−1, 1], which describes a rotation by π of the coordinate system about the x-axis. Using this rotation, the forces exerted on the node at y = n by the beam that lies on n − 1 ≤ y ≤ n (which we denote Y − ), are With Newton's second law, the time-harmonic equation of motion of the n th node can be expressed where M = diag [1, µ] is the matrix that describes the inertial properties of the chain. The first component of M corresponds to the mass of the junction whilst the second component corresponds to the rotational inertia, denoted µ. In the chosen system of natural units µ is actually the moment of inertia per unit mass of the nodes and so, for junctions between thin beams, 0 < µ 1. After applying the discrete Fourier transformation where k is the Fourier variable and u F (k) = [W F (k), Θ F (k)] T is the generalised displacement in reciprocal space, formed of the Fourier transformation of the real space displacements w(n) and θ x (n).

The dispersion equation
The solvability condition of equation ( Given that equation (3.5) is quadratic in ω 2 , we can find analytical representations ω = ω(k) for the dispersion curves, such that σ(ω(k), k) = 0. Moreover, the exact values of the edges of the band gaps can be found easily. The dependence of the dispersion equation on µ illustrates the importance of the rotational inertia; in particular the rotational inertia can be used as a parameter to control the location and width of the band gaps.
In Figure 8 we plot the dispersion diagram for a chain with µ = 0.06, across the first Brillouin zone k ∈ [−π, π]. It is seen that, for µ = 0.06, the Euler-Bernoulli beam chain has two band gaps, one finite band gap for 4 √ 3 < ω < 10 2/3, and the semi-infinite band gap associated with discrete systems, for ω > 10 √ 2.  In Figure 9, we explore the evolution of the dispersion curves with changing rotational inertia. It is emphasised that Figure 9 does not show dispersion surfaces, rather planes have been formed by plotting the dispersion curves for a range µ. The planes are disjoint except for at one degenerate value of µ = 1/12 when ω = 4 √ 3. For any µ = 1/12, the chain always has two band gaps: a semi-infinite band gap for large ω, and also a finite band gap with non-zero width. Even as µ → ∞, a finite band gap still exists, although it becomes infinitesimally thin. Furthermore, regardless of the value of µ, the dispersion equation always satisfies σ(4 √ 3, π) = 0. This provides a connection to the square lattice (c.f. § 2) since σ(4 √ 3, 0, π) = σ(4 √ 3, π, 0) = 0 and, in this regime, the 2D problem becomes quasi-one-dimensional. When µ = 1/12 the semi-infinite band gap remains and, despite the fact there is no finite band gap, forcing the chain at this degeneracy nevertheless results in a localised mode, as discussed in § 3.8. In the vicinity of µ = 1/12, the critical point (ω, k) = (4 √ 3, π) is associated with the lower/upper boundaries of the finite band gap for µ ≶ 1/12, and for µ 1/12, ω = 4 √ 3 forms the lower boundary of the semi-infinite band gap.

The Green's function
To construct the Green's function, we consider the application of force f ∈ C 2 to the central node of the lattice. The forcing vector f = [f w , f θ ] T has components corresponding to the translational force f w and rotational moment f θ , and we let f w , f θ ∈ {0, −1} with the exclusion of f w = f θ = 0 which equates to the absence of an applied force. The Green's function in reciprocal space can then be expressed For convenience, we separate the flexural displacement W F (k) and rotation Θ F (k) components for evaluation. Explicitly, the flexural displacement of the chain in reciprocal space is The denominator of W F (k) coincides with the dispersion equation (3.5) as a result of the inverted matrix in equation (3.6). The flexural displacement in direct space is obtained by applying the inverse Fourier transform to the spectral representation equation (3.7), as follows

8)
A significant difference between the chain of beams and the square lattice from § 2 is the ability to evaluate the integral of the inverse Fourier transformation in closed form. To evaluate the integral we use the substitution z = e ik to map the line segment [−π, π] to the unit circle in the complex plane. Following the substitution, we denote the displacement W F (z), and apply Cauchy's Residue Theorem by determining where the poles of the integrand, z i , lie in respect to the unit circle and taking the residues as follows, The numerator of equation (3.7) is an entire function and, therefore, the poles of equation (3.8) are the zeros of the dispersion equation (3.5) and, consequently, are frequency dependent. Nevertheless, we can use the fact that the complex solutions of the dispersion equation only cross the unit circle when the frequency traverses the boundary of a pass band to construct closed form Green's functions for each frequency regime: finite band gap, semi-infinite band gap and pass band frequencies.

Finite band gap
As an example, we evaluate explicitly the flexural displacement Green's function in the finite band gap regime for translational forcing f = [−1, 0] T , using equation (3.9). For frequencies in the finite band gap, where the following repeated factors have been introduced In Figure 10, we plot w(n) for the chosen values of µ = 0.06 and ω = 7.5, it can be seen that the waves decay exponentially away from the forcing point at the origin as expected.

Rotational forcing
Choosing the forcing vector f = [0, −1] T in equation (3.6) has the effect of applying a clockwise moment about the x-axis, rather than the previous translational force. We then evaluate the real space Green's function in the same manner, using equation (3.9), to find the flexural displacement of the chain in response to the applied moment. This is plotted in Figure 11, using the same values for µ = 0.06 and ω = 7.5 as were used in Figure 10 for comparison.  Figure 11 clearly demonstrates the twist induced in the chain from the rotational forcing. There is notably zero translational motion at the n = 0 node, as expected for a purely rotational force, followed by a comparatively large displacement, anti-symmetric in each direction, which soon decays.

The semi-infinite band gap
Alongside studying the finite band gap, we also construct Green's functions for the semiinfinite band gap. Using the same value of the rotational inertia, the band gap Green's function for the flexural displacement in the case of translational forcing f = [−1, 0] T has been plotted in Figure 12 for ω = 14.3. Similarly, we evaluate the Green's function for the flexural displacement due to an applied moment f = [0, −1] T . This has been plotted in Figure 13 and exhibits the signature twist shape as expected.  Comparing Figures 10 and 12, it is seen that the displacements have significantly different shape despite having the same type of forcing. The same can be said comparing Figures 11  and 13. These differences can be explained by looking a the dispersion curves in Figure 8. For values of ω in the finite band gap, the complex solutions of the dispersion equation are of the form k = π + ik, wherek ∈ R, leading to oscillatory eigenmodes of the form e iπn−kn . In contrast, for the semi-infinite band gap, the complex solutions of the dispersion equation are purely imaginary, leading to eigenmodes of the form e −kn associated with non-oscillatory evanescent waves.

Pass band frequencies
While the Euler-Bernoulli beam chain has two pass bands, we note that in both the upper and lower pass bands, the position of the poles in relation to the unit circle, from equation (3.9), remains unchanged for fixed µ. Thus we can form one Green's function for both of the pass band regimes, which remains frequency dependent. The flexural displacement Green's function in response to a translational point force at n = 0 of frequency ω = 2, is plotted in Figure 14. As expected, we note that waves of constant amplitude are observed, indicating propagative modes.

Rotational displacement
Evaluating the inverse Fourier transformation of Θ F (k) component of u F (k) in equation (3.6) gives the rotational displacement of the chain θ x in response to the applied force. Until now, we have emphasised the flexural displacement component since the plots of w(n) can be interpreted intuitively as the deformed shape of the chain. However the rotational displacement Green's function is a useful tool that provides further insight into the lattice behaviour, and a quantitative measure of how much each mass rotates. In Figure 15 we plot θ x (n) in response to a translational point force f = [−1, 0] T for ω = 7.5 and µ = 0.06. These are the same conditions applied that were applied to Figure 10 when plotting the flexural displacement. In Figure 15, it can be seen that the n = 0 mass has zero rotational displacement under purely translational forcing, which is expected at the central node of the lattice where the force is applied. This agrees with Figure 10 which shows that while the beams either side of the n = 0 mass deform, the mass itself does not twist in place. It is also seen that the magnitude of the rotational displacement of each node decreases as n increases, which is characteristic of a decaying wave.

Degenerate point
As shown in § 3.1, the finite band gap collapses for only one value of the rotational inertia, µ = 1/12. When µ = 1/12, the dispersion equation degenerates and the two simple poles associated with the Green's function coalesce at ω = 4 √ 3. At this point of degeneracy, for a translational point force at n = 0, the inverse Fourier transformation of the flexural displacement has a convenient representation in terms of the PFQ regularised hypergeometric functionF as follows The displacement w(n) at the degeneracy results in an oscillatory evanescent wave, despite the frequency ω = 4 √ 3 being a solution to the dispersion equation, that is σ(4 √ 3, π) = 0, which would indicate a propagative mode. Other studies with detail on evaluating Green's functions using hypergeometric functions are found in [13 , 31 ] and references therein.

General Method
In this section we extend the previous work to form a general method for the construction of dynamic Green's functions on d-dimensional lattices with vector connections capable of describing a wide range of mechanical interactions, including elastic beams, rods, and springs. This method allows for the inclusion of interesting effects, such as the torsional stiffness and out-of-plane flexural deformations of the lattice connections, as was demonstrated in earlier sections.
The equations of motion are of the form where n ∈ Z d enumerates the lattice points in which the inertia of the system is concentrated, d ∈ N \ {0} is the dimension of space, u n ∈ C d denotes the generalised displacement which includes translation and rotation of the n th node, D is a discrete [tensor] operator which encapsulates the interaction (e.g. compressional, torsional, flexural, stiffness) of the lattice connections and M is a tensor describing the inertial properties of the system. The vector f n ∈ C d describes the [generalised] forcing -which incorporates, for example, linear forces, bending moments, and torsional moments -andü n is the second derivative of u n with respect to time. In the absence of forcing and for time-harmonic motion, of angular frequency ω, the equations of motion reduce to where 0 is the d-dimensional zero vector, and the common factor of e iωt has been omitted for convenience. We introduce the discrete Fourier transformation, where k is the d-dimensional For discrete systems this equation is polynomial in ω and the solutions of the dispersion equations can be plotted to produce dispersion diagrams and identify the regions where no real solutions to the dispersion equation exist, creating band gaps. In general, exciting the system at a frequency within a band gap will yield localised modes that decay exponentially with distance from the forcing point. We construct the Green's functions by considering the action of time-harmonic generalised forces at the origin of the d-dimensional lattices. In this case, the equations of motion are Du n − ω 2 Mu n = δ n,0 f n , (4.6) where δ n,m is the Krönecker delta, and we denote f 0 = f . Following the application of the Fourier transformation, equation (4.6) becomes The Green's functions describing the real space displacements u n are determined by applying the inverse of the discrete Fourier transformation as follows, Now, rather than applying external forces to the lattice systems, we consider creating localised modes through the introduction of an isolated defect in the inertial properties of the lattices at n = 0. In each system the defect is denoted by the d-dimensional tensor M which describes the altered inertial properties at the origin. The Fourier transformed equations of motion for lattices with a defect at the origin are then (4.9) By choosing an M such that it is possible to recover equation (4.7) and, in doing so, one can produce localised eigenmodes that coincide with band gap Green's functions. The displacement field corresponding to the localised defect modes can be found from equation (4.9), using the inverse Fourier transformation.

A mass defect for the Euler-Bernoulli beam chain
Using the above method, we derive a localised defect M for the Euler-Bernoulli beam chain such that the chain supports a localised mode for a chosen band gap frequency ω. The Green's functions w(n) and θ x (n), which form the vector u n can be evaluated in closed form for band gap frequencies using equation (3.9) and the analogous equation for θ x . Thus, u 0 = [w(0), θ x (0)] T from equation (4.10) is in-hand. It is noticed that in the case of a purely translational force, the central node experiences zero rotational displacement, that is, θ x (0) = 0 when f = [−1, 0] T , as can be seen explicitly in Figure 15. Likewise in the case of an applied moment f = [0, −1] T , we see that w(0) = 0, which is demonstrated in Figure 11. To avoid singular matrices, we define the functions α 1 and α 2 as As such, we define the inertial defect at the central node n = 0 using the matrix which is dependent on the applied forcing, and will result in a localised mode around the origin for the chosen band gap frequency. The function α 1 corresponds to a change in the mass of the central node n = 0, while α 2 corresponds to the changing the rotational inertia. Off diagonal elements in the matrices would correspond to coupling between the mass and rotational inertia; coupling of parameters is studied in other works such as [8 ] but is not relevant for this problem.

A mass defect for the square Euler-Bernoulli beam lattice
In the same manner as above, we can evaluate the inertial defect required for the 2D square lattice to support a localised mode for a chosen band gap frequency ω, coincident with the band gap Green's function for the desired forcing. Since the Green's function u (m,n) can be evaluated numerically at the 0 node, it remains possible for us to find the required values of w(0, 0), θ x (0, 0) and θ y (0, 0) that compose the u (0,0) vector. As above, to avoid singular matrices, we define the functions Following this, we define the inertial defect matrix as resulting in a localised defect mode about the 0 node. The component β 1 corresponds to a defect in the mass of the 0 node, while the β 2 and β 3 components correspond to defects in the rotational inertia for the x-and y-directions respectively.

Concluding Remarks
Firstly we studied a 2D square lattice of Euler-Bernoulli beams with junctions that possessed both mass and rotational inertia. The beams also possessed torsional stiffness, to account for the coupling between flexural and torsional waves at the junctions. The Green's function was used alongside finite element software to illustrate the extreme anisotropy (including remarkable asymmetric anisotropy, usually associated with PT-symmetry breaking) that can be achieved, and controlled in the lattice by altering the frequency and nature of the applied forcing. Analysis of the dispersion diagrams for various values of the rotational inertia and torsional stiffness was used to illustrate how the propagation of waves and dynamic behaviour of the lattice can be manipulated, giving an unprecedented level of control over the pass band frequencies and their preferential directions through the lattice. We also studied the related problem of a 1D chain of Euler-Bernoulli beams with junctions that possess mass and rotational inertia. Closed form analytical Green's functions were achieved for each of the pass band and band gap regimes of the chain, for different choices of forcing. We also demonstrated how the rotational inertia provides significant control over the propagating frequencies on the lattice. Lastly, we provided a method to construct Green's functions for a generalised d-dimensional lattice with discrete mass and arbitrary linear interactions between nodes. We then provided a method to design lattice defects such that the lattice possesses a localised mode of a chosen frequency. This method was used to derive inertial defect matrices for the 2D lattice and the 1D chain from earlier in the paper. The work in this paper has applications in many areas where the control of wave propagation through elastic structures is of particular interest -including metamaterials, seismic protection, energy dissipation, and others.