Elastodynamics on graphs—wave propagation on networks of plates

We consider the wave dynamics on networks of plates coupled along 1D joints. This set-up can be mapped onto an extension of wave graph systems studied in, for example, quantum graph theory. In the elastic case, different mode-types (flexural, longitudinal and shear waves) propagate in each plate and do so at different wave speeds. The flexural (or bending) modes are described in terms of fourth order equations introducing an always evanescent wave component into the system. Waves encounter plate intersections and can be transmitted, reflected or mode converted. The intersection or vertex scattering matrices mix different waves which can be propagating (open) and evanescent (closed). The local scattering matrices and the global transfer operator are no longer unitary; the consequences of this non-unitarity on secular equations and the Weyl law will be discussed. The findings are of relevance to describing complex engineering structures such as networks of beams and plates.


Introduction
We treat the problem of elastodynamic waves propagating on graph-like structures, or networks of plates, extending ideas developed in the context of quantum graphs [1][2][3]. The motiv ation for this investigation is two-fold. First, the elastodynamic extension of quantum graphs has novel and interesting theoretical features, emerging, for example, from the need to take account of evanescent modes which appear systematically in elastic problems. Second, the behaviour of elastic waves in complex structures is of practical importance in modelling noise and vibration in large structures such as ships or cars in the high frequency limit: in this context, the study of energy propagation on graphs can serve as a proxy for problems whose greater intrinsic complexity prevents explicit solutions, while sharing common characteristics such as a statistical description of the response of the system promising avenues to tackle otherwise intractable physical problems.
In the context of noise and vibration, numerical tools such as the finite element method (FEM) or the boundary element method (BEM) scale with the wavelength and lead to ever larger matrix problems in the short wavelength limit. Alternatively, ray-based [4][5][6] methods have been proposed and are used for modelling wave propagation within a structure. Wave methods based on the so-called wave-finite-element method (WFEM) have been useful in describing wave propagation within layered materials [7] and have found applications in mechanical engineering treating wave propagation problems within coupled beam and plate problems in terms of local scattering matrices at component interfaces and global operators connecting the local interface solutions [8,9].
Here, we consider elastodynamic waves propagating on graph-like structures such as networks of plates. We will restrict the analysis to homogenous and isotropic plates being extended to infinity in one spatial direction and being connected along 1D interfaces such as depicted in figure 1. Due to the translational invariance along the extended direction, we can project the dynamics onto graph networks consisting of 1D edges and vertices as shown in figure 1. Vertices correspond here to intersections between different plates or to plate boundaries (red dots in figure 1), while edges correspond to plate sections between intersections/ boundaries (blue lines in figure 1).
Scattering on graphs has so far mainly been studied for so-called quantum graphs [1][2][3] where the wave dynamics is given by the Helmholtz equation on 1D edges. These ideas have been extended to the Dirac operator and relativistic quantum mechanics on graphs propagating 2 or 4 spinor components [10] as well as to generic directed graphs with local unitary scattering matrices and averaging done over so-called unitary stochastic ensembles [11,12]. There is also a considerable literature on differential equations on graphs and their generalisations (see, for example, [13,14]) which has overlap with the underlying physical problems, although not with the scattering approach that is the focus here. Considering elastodynamics on graphs introduces novel features to the graph-scattering problem that are interesting both from a conceptual point of view and in terms of applications. First of all, plate equations are vector wave equations describing the linear dynamics of a displacement vector field. Using thin-plate theory, the vector field is split into a component perpendicular to the plate surface (bending modes) described by the fourth order biharmonic equation and in-plane modes (shear and pressure modes) described by the Helmholtz equation with mode specific wave velocities, Figure 1. Example of simple plate networks together with the corresponding graph description; the plates are considered to extend to ± infinity along the lines of intersection (marked as red dots in the cross section graph). see [15] for details. The biharmonic equation admits both propagating and evanescent solutions. We thus have in total four wave components describing the dynamics on each edge, of which at least one component is evanescent. This is in contrast to the vector wave dynamics on graphs for the Dirac operator considered in [10], for which the free wave solution on edges does not contain evanescent waves. Mode mixing can occur at the plate interfaces and reflection, transmission and mode-mixing coefficients will depend on the angle of incidence of the incoming wave front. The latter property has so far not been considered on wave graphs and can be introduced here due to the translational symmetry along the extended plates. By considering plates of finite thickness, the coupling coefficients at joints will not only depend on the number of plates which meet at the interface, but also on the angles between plates. (Wave dynamics on graphs depending on geometric features, such as angles between edges meeting at vertices, have also been termed 'fat graphs' in the literature [3,16].) In addition, there is a dependence on material parameters and the thickness of the plates as well as the driving frequency. The actual reflection/transmission coefficients can be obtained using methods described by Langley and Heron [17].
We will describe the general set-up for finding wave solutions on networks of plates in section 2. We will in particular focus on the role of evanescent contributions in the overall description of the wave propagation. We will introduce generalised scattering matrices at vertices or interfaces and derive a transfer matrix describing the wave dynamics on the network globally including evanescent contributions; the transfer matrices set up in this way are no longer unitary-in contrast to wave graph systems considered so far. It is often assumed that evanescent waves do not carry power and hence that their influence is negligible [17]. This is justified if the emphasis is on scattering into the far-field. In the presence of counterpropagating pairs of evanescent waves, these waves do, however, contribute to the power transmitted across plate sections [18] as will be shown. We will then give a brief introduction into elastic plate theory and describe the local interface scattering matrices in more detail. In section 3, we will present generalised scattering matrix conditions similar to those considered by Smilansky and co-workers in [19,20]. We will then investigate how the evanescent contributions influence the eigenvalue conditions and the so-called functional equations and we will construct a unitary transfer matrix for the overall graph. In section 4, we will derive the Weyl term describing the mean density of eigenmodes on the graph and discuss the influence of evanescent contributions on Weyl's law. Conclusions and an outlook on applications is given in section 5.

Wave dynamics on plate networks-the general set-up
We will give here the general set-up for describing the wave dynamics in networks of plates. The plate segments between two intersections are considered isotropic and homogeneous, with plates being of constant width. We furthermore assume that the underlying elastic equations are linear and lossless in what follows. We work in the frequency domain with fixed frequency ω corresponding to a time dependence of the form e −iωt . The plates are of infinite extent in one direction, here the x direction, reducing the set-up to a graph structure with effectively 1D edges-the plates-and 0D vertices-the plate intersections or free plate boundaries. In what follows, we will briefly recapitulate the basic ingredients for the wave dynamics on a network, before discussing the wave modes describing the elastic deformations of plates and then the coupling coefficients for reflection/transmission at plate intersections.

Waves on networks
Wave propagation on networks is usually decomposed into free wave propagation on edges followed by scattering on vertices leading to coupling between different edges [1][2][3]. The wave field is described as a superposition of incoming and outgoing waves at vertices. The propagation along edges is written in terms of a shift matrix containing the phase shifts experienced by each mode, while coupling between different edges at vertices is given in terms of vertex scattering matrices derived from boundary conditions. A similar set-up has been proposed in terms of boundary integral equations in [21]. A semiclassical version valid in the short-wavelength limit has been presented by Bogomolny [22]. We follow the two step approach sketched above: first we determine the possible solutions of the wave equation on each edge, disregarding the boundary conditions at the vertices. The wave field on each edge is then split into incoming and outgoing solutions with (multicomponent) wave amplitudes at vertices collected in vectors ψ − and ψ + . This is sketched for a specific example, a T-junction (or T-beam), in figure 2. In the notation of figure 2, the components of the matrix on edge (ij) mapping outgoing waves at vertex i into incoming waves of the same type at vertex j is written in the form with 4D wave amplitude vectors ψ − ji and ψ + ij and S (ij) the so-called shift matrix of edge (ij) of the form Here, X = B 1 , B 2 , S or L denotes a specific mode type, such as one of two bending modes (B 1 and B 2 ), shear (S) or pressure/longitudinal (L) waves on plates, see section 2.2. Here, k X,⊥ is the associated wavevector component normal to the plate junction and L ij is the length of the plate segment perpendicular to the intersection. We have L ij = L ji . Note that the wave numbers can be either real and imaginary depending on whether we are dealing with propagating or evanescent modes. The shift matrices S (ij) describe wave propagation across individual plates between scattering events at plate junctions. At a plate junction i, (elastic) waves are partially reflected and partially transmitted as represented by a local scattering matrix denoted σ (i) , which transforms incoming waves into outgoing waves according to where the sum is over all vertices k connected to vertex i. The local scattering matrices σ (i) at a given intersection i are discussed in more detail in section 2.4. In a next step, we collect the local matrices S (ij) and σ (i) to obtain a global matrices S and σ describing connections across the entire graph, so that the global state vectors representing incoming and outgoing waves satisfy (when the system is undriven). Here, ψ ± denotes a vector of incoming or outgoing wave amplitudes of dimension dim ψ ± = Nm, where N is the number of edges or plate segments and m the number of modes on each segment. We have m = 4 in the case of plate elastodynamics, see section 2.2. We further define a transfer matrix [21,22] so that eigenmodes of the system are determined by the equation where T is a function of ω. The roots of the secular equation are then the resonant-or eigen-frequencies of the graph.

Wave dynamics on plates
To describe the local vibrational dynamics on a given plate segment α = (ij), a local coordinate system (x, y, z) α is defined so that the x-axis is along the infinite extension of the plate, the y-axis points into the plate segment and the z-axis points in the direction normal to the plate. The x-axis is a global coordinate, which is the same for all plate segments. We furthermore introduce the vector of elastic deformations (u, v, w) α for each plate which are defined with respect to the local coordinates (x, y, z) α such that u and v are in-plane deformations parallel to x and y, respectively, and w is parallel to z. The equations governing deflections of a plate α are then [15] for the out-of-plane elastic deformations w and for the in-plane elastic deformations u and v. Here, ρ is the density (that is, the mass per unit area), ν is the Poisson ratio, E is Young's modulus and is the flexural rigidity, with h being the thickness of the plate. These parameters can be different for different plate elements α. Furthermore, incoming and outgoing waves meeting at a given junction have the same dependence exp (ik x x − iωt) along x and with t. Together with the dispersion relations given below, this fixes the y-component of the wavevector, k y , for each wave mode and thus also the angles of incidence/reflection at plate junctions.
The equations of motion (8)- (10) give rise to four types of wave solutions [15] with wave numbers The subscripts L, S and B refer again to longitudinal, shear and bending modes, respectively. The wavevector components perpendicular to interfaces and pointing into plates are denoted, for some fixed component k x along the junction. Note that one of the solutions for the bending deformations, k B2,y , is always purely imaginary, describing an evanescent mode. The other modes can be either propagating or evanescent depending on whether k x is smaller or larger than the wave number k X in (11) with X = B, S or L. We will use the symbol κ to denote the positive imaginary component of the evanescent wavevector components. For example if, in addition to the always-evanescent bending mode B 2 , the longitudinal mode L is also evanescent, we denote In the following, the discussion is valid for any possible combination of open and closed modes, but we refer to the case in equation (15) predominantly when considering specific examples.
Next, we set out our conventions for the wave amplitudes. The in-plane displacements u and v can be written as while the out-of-plane displacement takes the form w = a − B1 e −ikB 1 ,y y + a − B2 e κB 2 ,yy + a + B1 e ikB 1 ,yy + a + B2 e −κB 2 ,yy e ikxx .
The symbols + and −, respectively, label outgoing and incoming components at a junction, as in section 2.1, and we note that the conventions used for evanescent modes are chosen to align the sense of propagation with the direction of decay.

Power flow
When the problem is lossless, the scattering, shift and transfer operators exhibit symmetries associated with conservation of power. To state these concretely we relate the power flux incident on a junction to the wave amplitudes a ± X . Following the analysis in appendix A, this can be achieved more compactly by rescaling the amplitudes as follows: For evanescent modes we choose the branch √ i = e iπ/2 of the ensuing complex square root (and note that it is only for the mode B 2 that we alternate the sign before the square root when swapping between incoming and outgoing components-this leads to simpler expressions for power flux in the discussion below). As described in appendix A, the detailed form of the power flux depends on which modes are propagating and which are evanescent. We illustrate the calculation here for the case where the longitudinal mode is evanescent, along with the always-evanescent bending mode B 2 , while the shear wave and bending mode B 1 are propagating. Then, in terms of the scaled amplitudes, the power contributions arriving at a junction from the in-plane and out-of-plane modes in a given plate α = (ij) are respectively Equation (20) takes into account that counterpropagating evanescent waves can carry flux, the consequences of which are explored further in section 3. In open scattering systems there are no such contributions because evanescent waves can only propagate in one direction and flux is carried by open modes only. However in a scattering problem of finite extent such as this one, solutions which grow exponentially as we cross a plate from one junction to another are admissible, and we get a contribution to the current from the closed modes as described by the last terms in (20).
In the remainder of this paper we will assume that the solution vector on each plate α = (ij) is written in terms of the scaled amplitudes and, following [19,20], we will also order the modes on each plate into propagating 'open' and evanescent 'closed' components, using the notation Then the net power flux j α = j in-plane α + j bending α given by (20) can be written with the obvious interpretation of the inner product. With the conventions thus adopted, the form of the net power, equation (23), is found to apply generally, irrespective of which modes are evanescent or propagating, see also appendix A.

Scattering from plate junctions
Having defined the wave amplitudes as described in sections 2.2 and 2.3, the local shift matrix S α and the global shift matrix S are now obtained with respect to ψ ± α from (21) with entries defined in (2) using the wave numbers (12)- (14). To complete the construction of the transfer operator (5) and thus to obtain global solutions on the network, we need also to determine the local scattering matrices (3) at plate junctions in the basis (21). Reflection, transmission and mode-conversion coefficients at junctions with an arbitrary number of attached plates can be obtained following Langley and Heron [17], see also [23]. We will not reproduce the derivations here, but a typical set of coefficients is shown in figure 3 for an L-junction and an incoming bending wave (incoming from plate 1) using the line-junction approximation presented in [17]. (This approximation neglects geometric details of the junction itself). In figure 3, the relevant parameters, Young's modulus E, the area density ρ, the thickness h and the frequency ω, are chosen such that k L /k B = 0.3 with Poisson ratio ν = 0.3 reproducing one of the examples in [17]. The transmission coefficients τ XX ij with i = 1 and j = 1, 2 and X = B, X = B, S, L, are plotted here against cos φ, with φ the angle between the incoming wavevector and the x direction tangential to the interface. The coefficients τ XX ij shown in figure 3 are obtained by taking the absolute value square of the corresponding scattering matrix elements of the scattering matrix (3). Note that the scattering matrix σ is an 8 × 8 matrix in the example considered here and all matrix elements need to be computed including those between propagating and evanescent waves as well as between incoming and outgoing evanescent contributions. In figure 3, only the coefficients for propagating outgoing modes are shown-these form a unitary sub-matrix of the full scattering matrix, see section 3 for details. The corresponding transmission coefficients for a fixed incoming mode thus sum up to one. Note that the outgoing longitudinal mode becomes evanescent for cos φ > 0.3, the critical angle for the shear wave is cos φ ≈ 0.507.
In what follows, the exact form of the scattering matrix elements are not important and the main results in section 3 can be derived based on general principles such as energy and flux conservation. For all examples and the analysis of the generalised Weyl law in section 4, we use realistic plate parameters and associated scattering coefficients following [17].

Extended unitarity of the scattering matrix
The key conditions that extend the unitarity of scattering matrices to include evanescent channels or modes have been derived in [19,20] in the context of billiards coupled to waveguides. These conditions have also been set out separately in the electromagnetic context by Carminati et al in [24]. We review them here in the context of the transfer matrix T = σS, equation (5), with σ and S, the global interface scattering and shift matrix as defined in (4) offering a derivation based on flux conservation only.
First, in analogy with (22), we separate the global solution vectors into propagating 'open' and evanescent 'closed' components, and denote a corresponding block-decomposition of the scattering and shift matrices by Recall that there is at least one closed bending mode for each plate, and possibly more if k x is large enough. We also define matrices representing projection onto the open and closed subspaces, respectively. Let us now focus on scattering from plate junctions. If this is a lossless process, then the current arriving at each vertex should be balanced by the current leaving it, so that the net current vanishes: extending the notation of (23) in the obvious way. This condition of vanishing net current can be rearranged to give Asserting that this condition must hold independently of ψ − whenever ψ + = σψ − brings us to the condition Written out explicitly for each block, this is equivalent to the conditions being satisfied by individual blocks (noting that (29) and (30) are redundant and that σ † oc ≡ (σ oc ) † etc). Thus, even though the open-open block σ oo is unitary, the full matrix is not. Bearing in mind that the detailed forms of (30) and (31) depend on the phase convention used for closed modes, these are the conditions written for higher-dimensional scattering problems by Rouvinez and Smilansky in [19] and by Schanz and Smilansky in [20], and separately in the electromagnetic context by Carminati et al in [24].
It can be shown similarly that, under the assumption that propagation across plates is lossless, the shift matrix satisfies the same condition: However, the product T = σS is in general not of this form: That is, conditions (27) and (32) do not define a group property. We remark, however, that by redefining the conventions for decomposing closed modes into incoming and outgoing components, similar to an approach suggested for higher-dimensional problems in [25], it is possible to construct scattering and shift matrices that are unitary and remain so under composition, although space limitations inhibit us from describing this in detail here and its discussion is deferred to a future publication.

A functional equation for the nonunitary transfer matrix
Resonant frequencies of the elastic network are obtained as solutions of the secular equation (7). For lossless problems, it is physically clear that we should find real solutions of this equation, but since T(ω) is a complex, non-unitary matrix, it is not immediately obvious how this emerges from the formalism. Here, we demonstrate a functional equation which shows why zeros are found for real frequencies.
By exploiting the identity which follows from (27) and (32) after some manipulation, we arrive at the functional equation where N and M are the number of open and closed components, respectively and we write P o S − iP c and σ † P o + iP c in block form before taking determinants. Since σ oo and S oo are unitary, it follows that | det(−σ oo S oo )| = 1 and the function is real-valued when ω is real and therefore generically has zeros on the real axis corresponding to the eigen-frequencies of the plate ensemble. Note that the unitary matrix σ oo S oo appearing in the denominator here is not the open sub-block T oo of T. Nor is it the unitary matrix obtained after elimination of closed modes, as described in section 3.3. Nevertheless the phase of its determinant is directly relevant to the characteristic polynomial of T, as above, and will be important when considering the Weyl formula in section 4.

Reduction to a unitary problem
Alternatively, the eigenvalue condition (6) can be reduced to a unitary eigenvalue problem by elimination of the closed modes. To demonstrate how this is achieved, let us first consider the simplified problem where S is replaced by the identity matrix, so that we seek solutions of ψ = σψ.
By writing separate equations for open and closed components defined in (24), and eliminating the closed components, this can be restated as Then the matrix can be shown to be unitary following from equation (27). To see this use (30) to write and recall that σ oo is unitary, so it just remains to show that the matrix is unitary. This is easily seen using (31), that is, The full problem represented by equation (6) reduces on elimination of closed modes to It can be shown that V is also unitary using (27) and (32). The derivation is more involved because the conservation conditions satisfied by T are less simple (see (33)), so we omit the derivation here, see also [26] for details.
We can therefore also determine resonant frequencies as zeros of Although this function looks formally similar to (34) and has the same real zeros, the relationship between these two expressions is non-trivial. Off resonance, the function D V (ω) is derived from the characteristic polynomial of V which corresponds in the original setting to the generalised eigenvalue problem whereas the function D T (ω) defined in (34) corresponds in the original setting to the simple eigenvalue problem Although solutions to these problems coincide when the resonance condition λ = 1 is satisfied, they are otherwise distinct. An explicit relation between the functions D T (ω) and D V (ω) can be given exploiting the identity One obtains in particular and therefore where the relation V = σ oo V S oo defines V . In the high-frequency regime where evanescent components of the total wave solution decay rapidly, the components of the closed-closed block T cc are generally small and one has det(1 − T cc ) ≈ 1. Then the functions D T (ω) and D V (ω) would differ only by a phase. However, for moderate or low frequencies where evanescent terms contribute more significantly, or where evanescent waves incident on junctions can be scattered resonantly, we will see that D T (ω) and D V (ω) can behave quite differently despite necessarily having the same zeros. The differences are particularly large around frequencies supporting edge states, that is, states which are localised around junctions between plates and which decay exponentially away from the plates edges. These features and the ramifications of (38) for the Weyl law and a description of the spectrum in terms of periodic orbits are discussed further in section 4.

Weyl laws and trace formulas for elastic networks
In the previous section, we have presented two distinct secular equations for the lossless case given in equations (34) and (37). The eigenfrequencies are real in this case and can be expressed in terms of a density of states, or equivalently in terms of the staircase function N(ω; k x ) = #{ω n (k x ) < ω}.
Note that we consider here the spectrum as a function of ω for fixed wavenumber k x . We change the angle of incidence when changing the frequency and thus also the number of propagating versus evanescent modes. Figure 4 illustrates this, showing the number of propagating modes across the ω − k x plane for a plate with parameters given in the figure caption.
Here, white, yellow, green and blue correspond to the cases of no, one, two or three propagating modes, respectively. For fixed ω, there are no propagating modes if k x is large enough.
Increasing ω for fixed k x , one of the modes eventually switches on (bending mode B1 first for k x below about 23 m −1 for the plates considered in figure 4, and shear mode S first if k x is larger than this value). As ω increases further, a second and eventually a third propagating mode switches on, the order in which this happens depending on k x . An essential difference between the function D T (ω; k x ) in (34) and D V (ω; k x ) defined in (37) is that D T (ω; k x ) can have poles in addition to zeros on the real ω-axis, whereas D V (ω; k x ) has at most zeros due to the unitarity of V. The eigenvalues of the finite matrix V are then all on the unit circle and det(1 − V) is bounded. Because T is non-unitary, det(1 − T) is subject to no such restriction: we will find that poles can in fact arise which are associated with edge states and related to singularities of reflection of evanescent waves from plate edges or junctions. In what follows, we denote these poles by ν n (k x ) and we define a staircase function also for poles of the form P(ω; k x ) = #{ν n (k x ) < ω}.
We will now split these expressions into a smooth Weyl-like component and periodic-orbit contributions giving rise to an oscillatory part, see for example [2,3,22]. We will show that the smooth Weyl law obtained for D T (ω; k x ) differs from that obtained for D V (ω; k x ) in the way the pole-counting staircase function is included.

Decomposition of the staircase based on D T
We first describe a decomposition of the density of states using as a starting point the function D T (ω; k x ) defined in (34) with a corresponding decomposition based on D V (ω; k x ) discussed in section 4.2. We follow essentially the treatment described in Kottos and Smilansky [2].
To obtain the number of zeros and poles in the interval [ω 0 , ω 1 ] with ω 1 > ω 0 0 real, let C be a contour in the complex ω plane passing leftwards from ω 1 + i to ω 0 + i and ω 0 − i and then rightwards to ω 1 − i , before returning to its starting position (with > 0, small). From Cauchy's theorem, we have where N is the number of zeros on the real axis between ω 0 and ω 1 and P is the number of poles. Writing (34) in the form A Weyl formula is obtained by identifying the first term in brackets as providing the average growth rate of the staircase function in the form The remaining terms provide the fluctuations around this mean and can be written in terms of a trace formula [2]. Note, however, that N T (ω; k x ) corresponds to a staircase function that not only increases by one unit at each resonant frequency ω n , but also decreases by one unit at each pole ν n . This feature is seen in figure 5(a), where we show the staircase function for a T-junction configuration of three plates with parameter values similar to those of [27], but varied slightly on each plate so that we can resolve edge states supported at the tips of each plate and discussed at the end of this section. See the figure caption for the parameter values used. For quantum graphs we would typically substitute ω 0 = 0 in (40) and start counting states with N(ω 0 ; k x ) = 0. In the case of elastic graphs, however, care must be taken to properly account for changes in the numbers of propagating and evanescent modes as ω and k x are To leading order, the Weyl law is typically written in terms of volumes in 3D, areas in 2D or lengths in 1D-or more generally in terms of phase-space volumes-thus capturing the essential geometry of the problem; boundary conditions enter through next-to-leading contributions. In this spirit, we note that The contribution from arg(det(S oo )) gives the leading order term for the growth of N(ω; k x ) and corresponds to the total optical length of the graph taking account of all the open modes in the problem. In fact, from the description in section 2, we can write explicitly where the sum is taken over all plates and the open modes they support. The remaining terms derived from arg(− det(σ oo )) provide boundary corrections induced by the conditions imposed at junctions. These are nontrivial when the full physical boundary conditions are used and we cannot offer simple analytical expressions. As discussed in section 2.4, they can be approximated in terms of simple low-dimensional systems of equations using the treatment suggested in [17]. Denoting the fluctuating part of the staircase function related to D T as we can identify this part with periodic-orbit contributions after writing noting that Tr T t can be written as a sum over all contributions from periodic paths on the graph with period t. The expansion (41) is well-defined only if the eigenvalues of the matrix T lie entirely on or inside the unit circle: we show in the following sections that one can find examples of plate networks in which an eigenvalue of T associated predominantly with the closed-closed block T cc lies outside the unit circle and then the second form above cannot be used. When this expansion is admissable then Ñ T (ω) can be expressed as a sum over periodic orbits including orbits with evanescent segments or even being fully evanescent on all plates.
To illustrate the appearance of edge states and associated singularities in D T , we return to the staircase function shown in figure 5(a) obtained for a T-junction graph. The range of frequencies shown is chosen to coincide with an interval in which there is one propagating mode (bending in this case) on each plate while all other modes are evanescent-corresponding to a horizontal slice through the yellow region at the bottom left of figure 4. We have intentionally chosen a plate thickness here that is not particularly small relative to the plate lengths, despite the underlying physical model being based on thin-plate theory. This has been done so that features relating to poles of D T (ω) are easier to resolve. We emphasise that qualitatively similar features are seen also in figure 7, which uses smaller thicknesses.
We find that D T has three poles towards the end of the illustrated frequency interval, in which only the bending mode B 1 is propagating. These are associated with poles in the closed block σ cc of the junction scattering matrix, representing singular reflection of in-plane modes from the tips of each of the plates in turn. These in-plane modes (S and L) are evanescent and decouple from the bending mode while scattering from the plate tips. Because they are evanescent, we find that it is possible for reflection from the tips to be singular for these modes (see appendix C). This occurs in our illustration for a slightly different frequency on each of the three plate tips (because they each have slightly different values of Young's modulus E), thus giving rise to three poles nearby in frequency. We find systematically that for each such pole there is a nearby resonant frequency whose corresponding eigenvector is heavily localised at the associated plate tip, as an 'edge state'. Such edge states may be simply understood in the limit of large plate lengths and are illustrated schematically in figure 6. Then any inplane wave component approaching the plate tip evanescently must have vanishingly small amplitude when it reaches the tip itself. However, near the pole of the tip-scattering matrix, the amplitude of the corresponding reflected wave can be boosted arbitrarily strongly-note that this reflected wave is also in-plane and decays away from then plate tip, represented schematically by the black curve in figure 6. The boost provided by such singular reflection can be enough to overcome the smallness of the incoming wave and lead to a nontrivial overall solution that is strongly localised near a plate tip. Because the corresponding resonant frequency is extremely close to the pole frequency, the staircase function in (39) rises and falls again in a short frequency interval, leaving no net change to the staircase, and such edge states are essentially 'missed' by the Weyl formula defined by D T (ω; k x ). We will see in the next section that these edge states are positively counted by the Weyl formula defined through D V (ω; k x ).
We may quantify the frequency gap ∆ω between a pole and its corresponding edge resonance by noting that an evanescent wave leaving the central junction decreases in amplitude by a factor of O(e −κL ) when it arrives at a plate tip, where L is the length of the plate and κ the imaginary part of the wave number for the slowest-decaying in-plane mode. This wave is boosted by a factor of order σ = O(ω/∆ω) by reflection near the pole frequency, before suffering a further O(e −κL ) decay on its way back to the central junction. Assuming scattering amplitudes from the central junction are O(1), we then require O(e −2κL ω/∆ω) = O(1) or ∆ω/ω = O(e −2κL ) for a stationary state to arise. The parameters in figure 5 have been chosen so that this gap is not particularly small, but the gap may in practice be undetectably small if the plate lengths are large enough.

Decomposition of the staircase based on D V
An alternative decomposition of the density of states into a Weyl term and periodic orbit contributions can be obtained from D V (ω; k x ) defined in (37). In this case, equation (39) is replaced by where det(−V) = e iΘV (ω;kx)/2 defines the new phase Θ V (ω). Note that the corresponding Weyl law Figure 6. A schematic illustration is provided here of the mechanism producing edge states near poles of D T . The blue and black curves respectively represent evanescent waves approaching and leaving the plate edge on the right. Edge states occur when an evanescent wave approaches the plate edge with a very small amplitude (proportional to the blue curve). Singular scattering associated with a pole of σ can then boost the reflected amplitudes sufficiently that the wave decaying to the left (proportional to the black curve) is significant, leading to a global stationary state that is strongly localised on the corresponding edge.
is different to (40), where the difference manifests itself in the way periodic orbits made up entirely of evanescent segments are treated. This can be quantified by considering equation (38), from which the relation follows. The expansion in the last line is again only possible if all eigenvalues of T cc are inside or on the unit circle. In the D T representation, the fluctuating part of the staircase function contains all periodic orbits, including those periodic orbits made up entirely of evanescent contrib utions (contained in Tr T t cc ). In the V representation, these all-evanescent periodic orbits are counted as part of the smooth contribution N V . Given that these contributions are indeed non-oscillatory, an ambiguity of where to 'count' the Tr T t cc contributions is perhaps not surprising. It then becomes a matter of convention whether one prefers the reduction of the problem to a unitary matrix V, making the conservative nature of the dynamics immediately manifest, or rather prefers the full, nonunitary matrix T, resulting in a more equitable treatment of open and closed, i.e. evanescent, periodic orbits. The latter also behaves less pathologically when eigenvalues of the block T cc fall outside the unit circle, as discussed below. The difference between the two representations vanishes at high-enough frequencies (for fixed k x ) when the evanescent contributions are exponentially small in the lengths of the plates. Comparison between the staircase function N (blue) and the Weyl terms N T (yellow) and N V (red) for a T-junction using the same parameters as in figure 5 apart from h = 0.01 m. As in figure 5, we find edge states and associated poles, here at frequencies around 1.6 × 10 4 Hz leading to deviations between N T and N V . One also observes a sudden drop of N V compared to N T around 0.25 × 10 4 Hz: this is associated with a transition of an eigenvalue of T from outside to inside the unit circle, see the main text for details.
Significant differences between these two Weyl counting functions arise at low-to-medium frequencies in two ways. The first possibility has already been discussed in section 4.1, allowing for singularities in the T representation which are absent in the V representation. Eigenstates of the elastic graph system associated with edge states are linked to poles of the local scattering matrices σ and give no net contribution to the staircase function N T (ω; k x ) and are thus not contained in N T (ω; k x ); they are, however, counted in N V (ω; k x ) and appear in N V (ω; k x ). This difference is apparent when comparing figures 5(a) and (b). It can also be observed in figure 7 for frequencies just above 1.6 × 10 4 Hz where N T does not pick up a series of edge-states, while N V does. The computations for figure 7 are done for the same T-junction structure and plate parameters as in figure 5 except that the thickness of the plates is smaller, leading to significantly more resonant frequencies being found in the interval concerned.
In figure 7, a second scenario arises where the two Weyl formulas differ. The sudden drop seen in N V at a frequency ω c ≈ 0.25 × 10 4 Hz is not present in N T . A closer inspection shows that an eigenvalue of T cc lies outside the unit circle for frequencies below ω c in this case, moving inside the unit circle while passing close to unity as ω passes through ω c from below. Near ω c , there is a resonant mode which is predominantly evanescent and, as with the edge states discussed in section 4.1, supported near plate tips. There is, however, no associated pole in this case and N T (ω; k x ) counts this resonance positively. The eigenmode arises here because, low in the spectrum, exponential decay of evanescent waves along plates may be slow enough that the spectrum of T cc reaches outside the unit circle without singularities in scattering arising. A simple analogue calculation, using 2 × 2 matrices, is outlined in appendix D which shows qualitatively similar behaviour.
The Weyl counting function N T (ω; k x ) passes smoothly through this event: the decomposition (39) remains valid as long as we do not invoke the expansion in (41) and provided we track the branch of the logarithm continuously. There is no associated pole, so the staircase function N T (ω; k x ) counts the state even though it is supported predominantly on plate edges. We emphasise that (42) similarly remains valid throughout: there is simply a transfer of one unit from the Weyl term to the periodic-orbit contributions.
The main conclusion of this section is that comparison between the Weyl counting functions N T (ω; k x ) and N V (ω; k x ) allows us to distinguish between predominantly evanescent edge states and predominantly propagating bulk states. This makes it possible to count these two types of states separately without necessitating a detailed examination of the spectrum or the eigenfunctions.

Conclusions
We have established an extension of quantum graphs to the elastic case, motivated by the relevance of these structures to problems in noise and vibration and by the novelty of some of the theoretical features. Elastic graphs differ from quantum graphs due to the inclusion of the evanescent modes, which make scattering matrices and the transfer operator non-unitary in the most natural representation. We have demonstrated how underlying symmetries associated with flux conservation lead to an extension of unitarity which allows us to establish a functional equation satisfied by the secular equation for eigenstates. In particular, this functional equation can be used to establish a decomposition of the density of states (or corresponding staircase function) into a smooth Weyl-like part and periodic-orbit contributions accounting for periodic orbits that are fully or partially or not at all evanescent. Two alternative approaches have been explored which differ in the way in which they count edge states. In particular there is a difference between the corresponding Weyl counting functions which allows us to count edge states separately from bulk states. modes S and B 1 to be propagating. We find in this case, after some manipulation and using the dispersion relations (11), that the in-plane power flux becomes while the bending power flux can be written as Represented in terms of the rescaled amplitudes b ± X , this reduces to equation (20). In other cases of propagating and evanescent modes, similar calculations lead to the general expression for the current given in (23).

Appendix B. Time-reversal and parity
In addition to the conditions (28)-(31) imposed by flux-conservation on scattering and shift matrices, we get analogous but distinct conditions respectively imposed by time-reversal and parity symmetries, which are described in this section.

B.1. Time-reversal symmetry
Time-reversal symmetry emerges from the observation that, if equations (16)-(18) provide a solution to the junction-scattering or plate propagation problems, then so too do their complex conjugates T(u, v, w) = (u * , v * , w * ), written explicitly as Note that when the solution is conjugated, incoming components of propagating modes are mapped to outgoing components and vice versa; incoming components of evanescent waves, however, are mapped back to incoming components (and outgoing to outgoing). In addition, we need to take account of the different behaviour of in-plane and out-of-plane amplitudes when the wavevector is reversed. The displacements defined by the in-plane amplitudes (a ± L , a ± S ) rotate with the direction of propagation and we note in particular that, in contrast to the bending components, these change sign under a reversal (k x , k X,y ) → (−k x , −k X,y ) of the wavevector. This distinction motivates us to define the diagonal matrix where the negative sign is used for in-plane components and the positive sign for out-of-plane, bending components. We also decompose this matrix into blocks corresponding to open and closed modes J = J oo 0 0 J cc in analogy with (25). Let ψ ± and ϕ ± respectively denote decompositions into incoming and outgoing waves of the original, global solution and of its conjugate, using the convention defined in (19) and (21). Taking account of all of the effects just mentioned, along with the complex phase of closed modes in that convention, we find that Using the projection operators defined in (26), these can be combined in Let us now write the symmetry imposed on the junction-scattering matrix σ. Noting explicitly that σ depends on k x and that the sense of k x is reversed by conjugation, we can write by treating the conjugate solution afresh as a scattering solution in its own right. In order for these to be consistent for arbitrary ψ − the condition Then this condition can be rearranged in the form or written explicitly in block form as Note that this last equation generalises to elastic waves the standard condition for simple, quantum wave propagation that the scattering matrix is symmetric when there is time-reversal symmetry-see also the discussion at the end of appendix B.2.
Taking account of the transformation P : (k x , k X,y ) → (−k x , k X,y ), the effect on the amplitudes a ± X is to change the sign of the shear mode but to leave the longitudinal and bending modes unchanged. This motivates us to define the matrix K = diag(· · · , ±1, · · · ), (B.12) in analogy with (B.4), but with the difference that the minus sign applies to the shear amplitudes only, whereas the matrix J in (B.4) reverses the sign of both the shear and longitudinal amplitudes. We also define L = JK, which reverses the sign of the amplitudes of the longitudinal modes but leaves those of shear and bending modes unchanged. Then, by treating P(u, v, w) as a scattering solution for either the junction-scattering or shift matrices with k x → −k x we get the identities and KS(−k x )K = S(k x ) (although this symmetry is trivial for the shift matrix as it does not couple different modes). Note that this parity symmetry, combined with the time-reversal symmetry expressed in appendix B.1, allows us to write the following alternative to (B.11) σ T (k x ) = Lσ(k x )L, (B.13) in which there is no reversal of k x . Thus, in contrast to the treatment of scalar quantum mechanics in [19,20], the elastodynamic scattering matrix is not symmetric when presented in the basis (21). Although σ can be made symmetric by an alternative choice of basis, the conventions which achieve this lead to more complicated expressions of the symmetries presented above.

Appendix C. Examples: reflection from an edge
We provide here some simple special cases of elastic scattering matrices which exemplify the flux conservation conditions (28)-(31) and the symmetries discussed in appendix B. We consider specifically the case of reflection from the edge of a single plate, which can be written explicitly without needing the more complex formalism of [17]. In particular, in-plane and bending modes do not couple in these examples and can be treated independently.
Consider first the case of bending waves normally incident on an edge. Then there is precisely one open and one closed mode, corresponding respectively to the labels B 1 and B 2 . It is then straightforward to derive edge scattering matrices of the forms for a clamped edge (satisfying w = 0 = ∂w/∂n) and for a free edge under conditions of no traction (leading to ∂ 2 w/∂n 2 = 0 = ∂ 3 w/∂n 3 for normal incidence) and using the conventions of (19). Note that both of these examples are symmetric, confirming as a special case the time-reversal symmetry (B.11) for bending modes (for which we can set J = I) at normal incidence (for which k x = −k x ). Parity symmetry is trivial for these examples. Edge reflection of normally-incident, in-plane modes is trivial so we treat non-normal incidence to exemplify this scenario. We consider first the case where k x is in the range where the shear waves are propagating but longitudinal waves are evanescent: see (15). Then one can show for tractionless reflection that 2 (C.1) in the conventions of (19), where X and Y are defined X = 2 k S,y κ L,y k x and Y = k 2 S,y − k 2 x and k L,y = iκ L,y defines κ L,y . It is straightforward to check that any matrix of the form given in (C.1) satisfies conditions (28)-(31) as long as X and Y are real. Note also that time-reversal symmetry is less trivial for these in-plane modes: the matrix is not symmetric in the conventions used, as the off-diagonal elements have opposite sign. Invariance under transpose is therefore only achieved by, for example, also reversing the sign of k x , which reverses the sign of X (see (B.11) and note that we can set J = −I for in-plane modes), or by conjugating with matrix L (see (B.13)). Note that the matrix therefore also satisfies the parity symmetry expressed in appendix B.2. Finally, we note that the explicit calculations in section 4 were performed in a regime where both of the in-plane modes are evanescent. Reflection from the free ends of the three plates is still described by the matrix at the right of (C.1), except that X = 2e iπ/4 √ κ S,y κ L,y k x = e iπ/4 X is now complex (while Y remains real), where k S,y = iκ S,y defines κ S,y , so that The open blocks are empty here and σ cc is the entire 2 × 2 matrix. Conditions (28)-(31) reduce to the condition that σ in-plane free be Hermitian, which is manifestly the case. Furthermore, (C.1) has poles when (X ) 2 = Y 2 , the condition for Rayleigh waves to be supported on the edge [28], which underlies the behaviour illustrated in figures 5 and 7.