Two-Term Spectral Asymptotics in Linear Elasticity

We establish the two-term spectral asymptotics for boundary value problems of linear elasticity on a smooth compact Riemannian manifold of arbitrary dimension. We also present some illustrative examples and give a historical overview of the subject. In particular, we correct erroneous results published by Liu (J Geom Anal 31:10164–10193, 2021).

Let (Ω, g ) be a smooth compact connected d -dimensional (d ≥ 2) Riemannian manifold with boundary ∂Ω ̸ = .We consider the linear elasticity operator L acting on vector fields u and defined by1 (1.1) Here and further on ∇ is the Levi-Civita connection associated with g , Ric is Ricci curvature, and λ, µ are real constants called Lamé coefficients which are assumed to satisfy2 µ > 0, d λ + 2µ > 0. (1.2) We will also use the parameter α := µ λ + 2µ .
(1.3) Subject to (1.2), we have ⊆ (0, 1). (1.4) We also assume that the material density of the elastic medium, ρ mat , is constant.More precisely, we assume that ρ mat differs from the Riemannian density det g by a constant positive factor.We complement (1.1) with suitable boundary conditions, for example the Dirichlet condition u| ∂Ω = 0 (1.5) (sometimes called the clamped edge condition in the physics literature), or the free boundary condition T u| ∂Ω = 0 (1.6) (sometimes called the free edge or zero traction condition in the physics literature, and also the Neumann condition 3 ), where T is the boundary traction operator defined by (1.7) Here n is the exterior unit normal vector to the boundary ∂Ω. 4  It is easy to verify that subject to the restrictions (1.2) the operator L is elliptic.Its principal symbol has eigenvalues (λ + 2µ)∥ξ∥ 2 (with multiplicity one), µ∥ξ∥ 2 (with multiplicity d − 1).
Here and further on ∥ξ∥ denotes the Riemannian norm of the covector ξ.The quantities λ + 2µ and µ are known as the speeds of propagation of longitudinal and transverse elastic waves, respectively.It is also easy to verify that either of the boundary conditions (1.5) and (1.6) is of the Shapiro-Lopatinski type  for L , and therefore the corresponding boundary value problems are elliptic 5 .
The boundary conditions (1.5) and (1.6) are linked by Green's formula for the elasticity operator, (1.8) where the quadratic form equals twice the potential energy of elastic deformations associated with displacements u, and is nonnegative for all u ∈ H 1 (Ω) and strictly positive for all u ∈ H 1 0 (Ω).The structure of the quadratic functional (1.9) of linear elasticity is the result of certain geometric assumptions, see formula (8.28)], as well as [CaVa-20, Example 2.3 and formulae (2.5a), (2.5b) and (4.10e)].
Consider the Dirichlet eigenvalue problem (1.10) subject to the boundary condition (1.5), where Λ denotes the spectral parameter.The spectral parameter Λ has the physical meaning Λ = (ρ mat / det g ) ω 2 , where ρ mat is the material density, det g is the Riemannian density and ω is the angular natural frequency of oscillations of the elastic medium.With account of ellipticity and Green's formula (1.8), it is a standard exercise to show that one can associate with (1.10), (1.5), the spectral problem for a self-adjoint elliptic operator L Dir with form domain H 1 0 (Ω); we omit the details.The spectrum of the problem is discrete and consists of isolated eigenvalues (1.11) enumerated with account of multiplicities and accumulating to +∞.A similar statement holds for the free edge boundary problem (1.10), (1.6), which is associated with a self-adjoint operator L free whose form domain is H 1 (Ω); we denote its eigenvalues by (0 ≤)Λ free 1 ≤ Λ free 2 ≤ . . . .
We associate with the spectrum (1.11) of the Dirichlet elasticity problem on Ω the following functions.Firstly, we introduce the eigenvalue counting function N Dir (Λ) := # n : Λ Dir n < Λ , (1.12) defined for Λ ∈ R. Obviously, N Dir (Λ) is monotone increasing in Λ, takes integer values, and is identically zero for Λ ≤ Λ Dir 1 .An analogous eigenvalue counting function of the free boundary problem will be denoted N free (Λ) 6 .
Secondly, we introduce the partition function, or the trace of the heat semigroup, (1.13) defined for t > 0 and monotone decreasing in t .The free boundary partition function Z free (t ) is defined in the same manner 7 .The existence of asymptotic expansions of N (Λ) as λ → +∞ and of Z (t ) as t → 0 + , and precise expressions for the coefficients of these expansions in terms of the geometric invariants of Ω, for either the Dirichlet or the free boundary conditions, and similar questions for the Dirichlet and Neumann Laplacians, have been a topic of immense interest among mathematicians and physicists since the publication of the first edition of Lord Rayleigh's 8 The Theory of Sound in 1877 .A detailed historical review of the field is beyond the scope of this article; we refer the interested reader to , , and , and references therein.
Before stating our main results, we summarise below some known facts concerning the asymptotics of (1.12) and (1.13), and their free boundary analogues.Further on, we always assume that (Ω, g ) is a d -dimensional Riemannian manifold satisfying the conditions stated at the beginning of the article.
Fact 1.1.For any (Ω, g ) we have (1.15) is the Weyl constant for linear elasticity, and Vol d (Ω) denotes the Riemannian volume of Ω.
We note that the expansions (1.19) do not in themselves imply the existence of two-term asymptotic formulae for N ℵ (Λ).However, formula (1.18) implies the following Fact 1.4.Let ℵ ∈ {Dir, free}, and suppose that we have (1.21) In general, the validity of two-term asymptotic expansions (1.20) is still an open question (as it is for the scalar Dirichlet or Neumann Laplacian).However, similarly to the scalar case, there exist sufficient conditions which guarantee that (1.20) hold.These conditions are expressed in terms of the corresponding branching Hamiltonian billiards on the cotangent bundle T * Ω, see  for precise statements.
Fact 1.5.Suppose that (Ω, g ) is such that the corresponding billiards is neither dead-end nor absolutely periodic.Then (1.20) holds for both the Dirichlet and the free boundary conditions.Fact 1.5 is a re-statement of a more general [Va-84, Theorem 6.1] which is applicable to the elasticity operator L since the multiplicities of the eigenvalues of its principal symbol are constant on T * Ω, as we have mentioned previously.
We conclude this overview with the following observation, see also Remark 4.1(i)].
Fact 1.6.The coefficients b ℵ are numerical constants which do not contain any information on the geometry of the manifold Ω or its boundary ∂Ω.Therefore, to determine these coefficients it is enough to find them in the Euclidean case.
Fact 1.6 easily follows from a rescaling argument: stretch Ω by a linear factor κ > 0, note that the eigenvalues then rescale as κ −2 , and check the rescaling of the geometric invariants and of (1.19).
Fact 1.6 allows us to work from now on in the Euclidean setting, in which case (1.1) simplifies to Note that we define the curl of a planar vector field by embedding R 2 into R 3 ; thus, curl curl applied to a planar vector field is a planar vector field.
We are now in a position to state the main results of this paper.Before doing so, let us introduce some additional notation. Let (1.23) The cubic equation R α (w) = 0 has three roots w j , j = 1, 2, 3, over C, where w 1 is the distinguished real root in the interval (0, 1).We further define γ R := w 1 . (1.24) Remark 1.7.The subscript R in γ R stands for "Rayleigh".Indeed, the quantity has the physical meaning of velocity of the celebrated Rayleigh's surface wave .The cubic equation is often referred to as Rayleigh's equation, and it admits the equivalent formulation (1.26) Equation (1.25) can be obtained from (1.26) by multiplying through by 4 (1 − αw) (1 − w)+(w −2) 2 and dropping the common factor w corresponding to the spurious solution w = 0.It is well-known [RaBa-95, ViOg-04] that for all α ∈ (0, 1) equation (1.25) -or, equivalently, (1.26) -admits precisely one real root w 1 = γ 2 R ∈ (0, 1).The nature of the other two roots w j , j = 2, 3, depends on α; we will revisit this in §D.2.
As we shall see, the Rayleigh wave contributes to the second asymptotic term in the free boundary case. ◀ respectively, where α is given by (1.3) and γ R is given by (1.24).
Theorem 1.8 will be proved in § §3-5 by implementing the algorithm described in §2.
Remark 1.9.The bound (1.4) guarantees that (1.27) and (1.28) are well-defined and real.◀ We show the appropriately rescaled (for the ease of comparison and to remove the explicit dependence on µ) coefficients b Dir and b free as functions of α in Figure 1.As it turns out, in odd dimensions the integrals in formulae (1.27) and (1.28) can be evaluated explicitly.
(i) Integrals in formulae (1.27) and (1.28) can be evaluated explicitly in even dimensions as well, but in this case one ends up with complicated expressions involving elliptic integrals.Given that the outcome would not be much simpler or more elegant than the original formulae (1.27) and (1.28), we omit the explicit evaluation of the integrals in even dimensions.
(ii) In dimensions d = 2 and d = 3, formulae for the second Weyl coefficient for the operator of linear elasticity both for Dirichlet and free boundary conditions are given in [SaVa-97, Section 6.3].The formulae in  have been obtained by applying the algorithm described below in §2, but the level of detail therein is somewhat insufficient, with only the final expressions being provided, without any intermediate steps.Our results, when specialised to d = 2 and d = 3, agree with those of  and allow one to recover these results whilst providing the detailed derivation missing in .

◀
Remark 1.12.Genquian Liu  claims to have obtained formulae for b Dir and b free .However, the strategy adopted in  is fundamentally flawed, because the "method of images" does not work for the operator of linear elasticity.Consequently, the main results from  are wrong.We postpone a more detailed discussions of , including the limitations of the method of images and a brief historical account of the development of the subject, until Appendix A. Below, we provide a preliminary "experimental" comparison of our results and those in .
Essentially,  aims to deduce the expression for the second asymptotic heat trace expansion coefficient b Dir in the Dirichlet case, as well as a corresponding expression in the case of the boundary conditions formula (1.5)] (called there the "Neumann" 12 conditions) which in our notation13 read n β ∇ β u α = 0. (1.31) We observe that the boundary conditions (1.31) are not self-adjoint for (1.10), as easily seen by simple integration by parts.Therefore, it is hard to assign a meaning to Liu's result in this case Theorem 1.1, the lower sign version of formula (1.10)].Nevertheless, even if one interprets the "Neumann" conditions (1.31) as our free boundary conditions (1.6), as the author suggests in a post-publication revision formula (1.3)], the result of  in the free boundary case remains wrong.
For the sake of clarity, let us compare the results in the case of Dirichlet boundary conditions only.The main result of  in the Dirichlet case is [Liu-21, Theorem 1.1, the upper sign version of formula (1.10)], which correctly states the coefficient a (cf.our formulae (1.15) and (1.17)), and also states, in our notation, that (1.32) This also implies, by (1.21), which differs from our expression (1.27) by a missing integral term.
For the reasons explained in Appendix A, formula (1.32) is incorrect.We illustrate this by first showing, in Figure 2  This ratio depends only on the dimension d and the parameter α.For each d , the ratio is monotone increasing in α (and is therefore monotone decreasing in λ/µ).
Dir /b Dir → 1 − in any dimension, see inset to Figure 2. Thus, for the smallest possible values of the Lamé coefficient λ, Liu's asymptotic formula would produce an almost correct result, however the error would become more and more noticeable as λ/µ gets large.
We illustrate this phenomenon "experimentally" in Figure 3 where we take Ω ⊂ R 2 to be the unit square.Neither the Dirichlet nor the free boundary problem in this case can be solved by separation of variables, so we find the eigenvalues using the finite element package FreeFEM [He-12].As Vol 2 (Ω) = 1 and Vol 1 (∂Ω) = 4, (1.20) in the Dirichlet case may be interpreted as N Dir (Λ) − aΛ ≈ 4b Dir Λ for sufficiently large Λ, and we compare the numerically computed left-hand sides with the right-hand sides given by our expression (1.27) and Liu's expression (1.33).As we have predicted, for λ = −1/2 both asymptotic formulae give a good agreement with the numerics, however for larger values of λ our formulae match the actual eigenvalue counting functions exceptionally well, whereas Liu's ones are obviously incorrect.Of course, the boundary of a square is not smooth, only piecewise smooth, but this does not cause problems because this case is covered by Theorem 1].Furthermore, Theorem 2] guarantees that sufficient conditions ensuring the validity of two-term asymptotic expansions (1.20) are satisfied. ◀ For an additional illustration of the validity of our asymptotics in the free boundary case, see  In this section we provide an algorithm for the determination of the second Weyl coefficient for more general elliptic systems.The algorithm given below is not new and appeared in  as well as in  for scalar operators, with [Va-84, §6] briefly outlining the changes needed to adapt the results to systems.However,  are not widely known and their English translations are somewhat unclear; therefore, we reproduce the algorithm here in a self-contained fashion and for matrix operators, for the reader's convenience.In the next section we will explicitly implement the algorithm for L Dir and L free .
Let A be a formally self-adjoint elliptic m×m differential operator of even order 2s, semibounded from below.Consider the spectral problem (2.1) where the B j 's are differential operators implementing self-adjoint boundary conditions of the Shapiro-Lopatinski type.
It is well-known that the spectrum of (2.1), (2.2) is discrete.Let us denote by the eigenvalues of (2.1), (2.2), with account of multiplicity, and let be the corresponding eigenvalue counting function.
In a neighbourhood of the boundary ∂Ω we introduce local coordinates (2.4) so that z > 0 for x ∈ Ω • , where Ω • is the interior of Ω.We will also adopt the notation (2.5) Let A prin (x, ξ) be the principal symbol of A and suppose that ξ ̸ = 0. Let h 1 (x, ξ), . . ., h m (x, ξ) be the distinct eigenvalues of A prin (x, ξ) enumerated in increasing order.Here m = m(x, ξ) is a positive integer smaller than or equal to m.
We will see in §3 that the above assumption is satisfied for the operator of linear elasticity L .
Theorem 2.2 ([Va-84, Theorem 6.1]).Suppose that (Ω, g ) is such that the corresponding billiards is neither dead-end nor absolutely periodic.Then the eigenvalue counting function (2.3) admits a twoterm asymptotic expansion for some real constants A and B B .Furthermore: (a) The first Weyl coefficient A is given by where n(x, ξ, Λ) is the eigenvalue counting function for the matrix-function A prin (x, ξ)14 .
(b) The second Weyl coefficient B B is given by where the spectral shift function is defined in accordance with and the phase shift ϕ(x ′ , ξ ′ , Λ) and the one-dimensional counting function N (x ′ , ξ ′ , Λ) are determined via the algorithm given below.
Step 1: One-dimensional spectral problem.Construct the ordinary differential operators A ′ and B ′ j from the partial differential operators A and B j as follows: (i) retain only the terms containing the derivatives of the highest order in A and B j ; (ii) replace partial derivatives along the boundary with i times the corresponding component of momentum: (iii) evaluate all coefficients at z = 0.
The operators ) are ordinary differential operators in the variable z with coefficients depending on x ′ and ξ ′ .
Consider the one-dimensional spectral problem (2.7) (2.8) Step 2: Thresholds and continuous spectrum.Suppose that ξ ′ ̸ = 0. Let h k (ζ), k = 1, . . ., m, be the distinct eigenvalues of (A ′ ) prin (ζ) enumerated in increasing order and let m k be their multiplicities, so that Clearly, for fixed (x ′ , ξ ′ ) we have In what follows, up to and including Step 6, we suppress, for the sake of brevity, the dependence on x ′ and ξ ′ .
Compute the thresholds of the continuous spectrum, namely, nonnegative real numbers Λ * such that the equation in the variable ζ has a multiple real root for at least one k ∈ {1, . . ., m}.We enumerate the m thresholds in increasing order The thresholds partition the continuous spectrum [Λ (1) * , +∞) of the problem (2.7), (2.8) into m intervals For Λ ∈ I (l ) , let k (l ) max be the largest k for which the equation (2.9) has real roots.Given a k ∈ {1, . . ., k (l ) max }, let 2q (l ) k be the number of real roots15 of equation (2.9).We define the multiplicity of the continuous spectrum in I (l ) as Step 3: Eigenfunctions of the continuous spectrum.At this step we suppress, for the sake of brevity, the dependence on l and write k max = k (l )  max , q k = q (l ) k , p = p (l ) .In each interval I (l ) denote the real roots of (2.9) for a given k = 1, . . ., k max by where the superscript ± is chosen in such a way that and the roots are ordered in accordance with be orthonormal eigenvectors of (A ′ ) prin (ζ) corresponding to the eigenvalues h k (ζ).Of course, these eigenvectors are not uniquely defined: there is a U(m k ) gauge freedom in their choice.
We seek eigenfunctions of the continuous spectrum (generalised eigenfunctions) for the one-dimensional spectral problem (2.7), (2.8) corresponding to Λ ∈ I (l ) in the form (2.10) where f 1 , . . ., f ms−p are linearly independent solutions of (2.7) tending to 0 as z → +∞, and the coefficients c @ k,q, j are not all zero.The coefficients c @ k,q, j are called incoming (@ = −) and outgoing (@ = +) complex wave amplitudes.
Step 4: The scattering matrix.Requiring that (2.10) satisfies the boundary conditions (2.8) allows one to express the outgoing amplitudes c + in terms of the incoming amplitudes c − .This defines the scattering matrix S (l ) (Λ), a p (l ) × p (l ) unitary matrix, via The order in which coefficients c @ k,q, j are arranged into p (l ) -dimensional columns c @ is unimportant.
Step 5: The phase shift.Compute the phase shift ϕ(Λ), defined in accordance with ϕ(Λ) := 0 for Λ ≤ Λ (1) * , arg det S (l ) (Λ) + s (l ) for Λ ∈ I (l ) . (2.11) The quantities s (l ) , l = 1, . . ., m, are some real constants whose role is to account for the fact that our construction of orthonormal bases for incoming and outgoing complex wave amplitudes involves a rigid (Λ-independent) unitary gauge degree of freedom, see Step 3 above.The branch of the multivalued function arg appearing in formula (2.11) is assumed to be chosen in such a way that the phase shift ϕ(Λ) is continuous in each interval I (l ) .
For each l , suppose that equation (2.9) with Λ = Λ (l ) * has a multiple real root for precisely one k = k (l ) , and that this multiple real root * is unique and is a double root 16 .Then the constants s (l ) in (2.11) are determined by requiring that the jumps of the phase shift at the thresholds satisfy (2.12) where j (l ) * is the number of linearly independent vectors v such that is a solution of the one-dimensional problem (2.7), (2.8), with f(z) = o(1) as z → +∞.
The threshold Λ (l ) * is called rigid if j (l ) * = 0 and soft if j (l ) * = m k (l ) .For rigid and soft thresholds formula (2.12) simplifies and reads minus for rigid and plus for soft.
Application of Steps 1-6 of the above algorithm to the elasticity operator with Dirichlet or free boundary conditions, which will be done in the next three sections, gives Theorem 2.3.

Second Weyl coefficients for linear elasticity: invariant subspaces
In this and the next two sections we will compute the spectral shift function for the operator of linear elasticity on a Riemannian manifold with boundary of arbitrary dimension d ≥ 2, both for Dirichlet and free boundary conditions, by explicitly implementing the algorithm from §2.This will establish Theorems 2.3 and 2.4.
In order to substantially simplify the calculations, we will turn some ideas of Dupuis-Mazo-Onsager [DuMaOn-60] into a rigorous mathematical argument, in the spirit of .Namely, we will introduce two invariant subspaces for the elasticity operator compatible with the boundary conditions, implement the algorithm in each invariant subspace separately, and combine the results in the end.
For starters, let us observe that the standard separation of variables leading to the one-dimensional problem (2.7), (2.8) can be achieved by seeking a solution of the form Next, suppose we have fixed ξ ′ ∈ R d −1 \{0}.Consider the pair of constant d -dimensional columns where 0 ′ stands for the (d − 1)-dimensional column of zeros.These define a two-dimensional plane Let us denote by Π the orthogonal projection onto P .Now, the principal symbol of the elasticity operator reads (3.1) Formula (3.1) immediately implies that the eigenvalues of the principal symbol are and of multiplicity m 2 = 1.It is easy to see that ξ ∈ P , P ⊥ ⊂ (I − ∥ξ∥ −2 ξξ T ) R d , and that P and P ⊥ are invariant subspaces of L prin .Furthermore, L prin P has two simple eigenvalues, (λ+2µ)∥ξ∥ 2 and µ∥ξ∥ 2 , whereas L prin P ⊥ has one eigenvalue µ∥ξ∥ 2 of multiplicity d − 2.
The above decomposition can be lifted to the space of vector fields.We define and be the one-dimensional operators associated with L and T , respectively; recall that the latter are defined by formulae (1.1) and (1.7).It turns out that the linear spaces V ∥ and V ⊥ are invariant subspaces of L ′ compatible with the boundary conditions.
Lemma 3.1 implies, via a standard density argument, that the operators L ′ ℵ , ℵ ∈ {Dir, free}, decompose as It is then a straightforward consequence of the Spectral Theorem that we can compute the spectral shift function for L ′ ℵ,⊥ and L ′ ℵ,∥ separately, and sum up the results in the end.More formally, we have Additional simplification: it suffices to implement our algorithm for the special case (3.10) The general case can then be recovered by rescaling the spectral parameter in the end, in accordance with In the next two sections, we assume (3.10).§4.First invariant subspace: normally polarised waves In this section we will compute the spectral shift functions shift ℵ,⊥ for the operators L ′ ℵ,⊥ , ℵ ∈ {Dir, free}.§4.1.Dirichlet boundary conditions.Consider the spectral problem The goal of this subsection is to prove the following result.
Lemma 4.1.We have and Here and further on χ A denotes the indicator function of a set A ⊂ R.
We shall prove Lemma 4.1 in several steps.The principal symbol (L ′ ⊥ ) prin , as a linear operator in P ⊥ , has only one eigenvalue The eigenvalue (4.3) determines the threshold which, in turn, yields exponents Therefore, the continuous spectrum of the operator L ′ ⊥ contains a single interval I (1) and the multiplicity of the continuous spectrum on this interval is p (1) ⊥ = 1.For Λ ∈ I (1) ⊥ , the eigenfunctions of the continuous spectrum read where (e j ) α = δ j α .Substituting (4.5) into (4.2) we obtain c 1 , . . . ,c d −2 ∈ C , other than the trivial one, from which the claim follows.
Proof.It is easy to see that the problem (4.1), (4.2) does not admit eigenvalues for Λ ≥ µ, i.e. eigenvalues embedded in the continuous spectrum.Furthermore, for Λ < µ a straightforward substitution shows that the only solution of (4.1), (4.2) of the form is the trivial one.This concludes the proof.
Combining formula (4.6), Lemma 4.2 and Lemma 4.3 one obtains Lemma 4.1.§4.2.Free boundary conditions.Consider the spectral problem (4.9) (4.10) The goal of this subsection is to prove the following result.
Proof.It is easy to see that the problem (4.9), (4.10) does not admit eigenvalues for Λ ≥ µ, i.e. eigenvalues embedded in the continuous spectrum.Furthermore, for Λ < µ a straightforward substitution shows that the only solution of (4.9), (4.10) of the form (4.8) is the trivial one.This concludes the proof.
Combining (4.11), Lemma 4.5, and Lemma 4.6, one obtains Lemma 4.4.§5.Second invariant subspace: reduction to the two-dimensional case In this section we will compute the spectral shift functions shift ℵ,∥ , ℵ ∈ {Dir, free}, for the L ′ ∥ .Calculations in the second invariant subspace are trickier, in that, unlike L ′ ⊥ , the operator L ′ ∥ is not diagonal.However, our decomposition into invariant subspaces implies the following (5.1)Fact (5.1) can be easily established by observing that, under assumption (3.10), elements in the domain of L ′ ∥ are of the form In the remainder of this section we will compute the spectral shift function for the operator of linear elasticity in dimension two.
The principal symbol L ′ plane prin has two simple eigenvalues18 These give us the two thresholds and the corresponding exponents so that the continuous spectrum [µ, +∞) is partitioned into the two intervals of multiplicities p (1) = 1 and p (2) = 2, respectively.
The normalised eigenvectors of L ′ plane prin are Hence, the eigenfunctions of the continuous spectrum for Λ in I (1) and (5.2) and (5.3) respectively.§5.1.Dirichlet boundary conditions.Consider the spectral problem (5.4) (5.5) The goal of this subsection is to prove the following result.
Lemma 5.4.The problem (5.4), (5.5) does not have eigenvalues, either below or embedded in the continuous spectrum.
Proof.Arguing as in the proof of Lemma 5.3, it is easy to see that thresholds are not eigenvalues.
For Λ ∈ I (1) we seek an eigenfunction in the form (5.2) with c ± 1 = 0, so that the solution is squareintegrable: (5.10) The latter satisfies (5.5) if and only if c = 0, which means there are no eigenfunctions for Λ ∈ I (1) , By looking at (5.3), it is easy to see that there are no (square-integrable) eigenfunctions corresponding to Λ ∈ I (2) , which completes the proof.
Lemma 5.5.We have and where γ R is given by (1.24), so that We shall prove Lemma 5.5 in several steps.Substituting (5.2) and (5.3) into (5.12)we get (5.13) We now have Lemma 5.6.
Proof.Arguing as in the proof of Lemma 5.6, it is easy to see that thresholds are not eigenvalues.
For Λ ∈ I (1) we seek an eigenfunction in the form (5.10).Substituting (5.10) into (5.12),one sees that the latter can only be satisfied if c = 0. Therefore, there are no eigenvalues in I (1) .
Lastly, by looking at (5.3), it is easy to see that there are no (square-integrable) eigenfunctions corresponding to Λ ∈ I (2) .This concludes the proof.
(i.e., the signs of the off-diagonal terms involving mixed derivatives change), and therefore the commutator [J , L ] does not vanish.Since the elasticity operator does not commute with reflections, the reflection method (or the method of images) is inapplicable.
The above argument shows that the principal symbol of the Laplacian (or the Laplace-Beltrami operator when working in curved space) is invariant under reflection, whereas the principal symbol of the operator of linear elasticity is not.This is what makes the method of images work for the Laplacian, but not for the operator of linear elasticity.The key difference between the Laplacian (or the Laplace-Beltrami operator) and the operator of linear elasticity is the presence of mixed derivatives in the leading term of the latter.
For the reader's convenience, let us spell out the precise points in  where the main mistakes occur, as a result of the operator L not being invariant under reflections.Below we use the notation from .
The author defines M := Ω ∪ ∂Ω ∪ Ω * to be the "double" of Ω, and P to be the "double" of the operator of linear elasticity in M .In the simplified setting of this appendix, M = R 2 and ) (the space of infinitely smooth functions with compact support), by a straightforward integration by parts one obtains Here (• , •) denotes the natural L 2 inner product and overline denotes complex conjugation.But (A.1) implies that P is not symmetric; therefore, it does not give rise to a heat operator.
As a result, the statement "P generates a strongly continuous semigroup e t P t ≥0 on L 2 (M ) with integral kernel K (t , x, y)" in [Liu-21, p. 10169, third line after (1.14)] is wrong, and all the analysis based on it breaks down (including formula (4.3)]).
In , the author states that they borrow their technique from .Indeed, the paragraph preceding formula (4.3) in [Liu-21, p. 10183] is taken, almost verbatim, from the beginning of Section 5 of p. 53].However, McKean and Singer applied the method of images to the Laplacian, for which the "double" operator is self-adjoint.
Let us conclude this appendix with a brief historical account.We note that the expression for b Dir was already found 19 in the 1960 paper by M. Dupuis, R. Mazo, and L. Onsager 20 [DuMaOn-60].Remarkably, this paper includes the critique of the 1950 paper by E. W. Montroll who presented exactly Liu's expression (1.32) for the second asymptotic coefficient, modulo some scaling, see [Mo-50, formulae (3)-( 5)].Dupuis, Mazo, and Onsager wrote, we quote: "Montroll . . .pointed out in 1950 a defect in the usual counting process of the normal modes of vibration and derived a corresponding correction term for the Debye frequency spectrum, . . ., proportional to the area of the solid.But it must be clearly realized that he used as a model a parallelopiped with perfectly reflecting faces, and that such boundary conditions 19 Not in a mathematically rigorous way, and for specific heat. 20The 1968 Nobel Laureate (Chemistry).
are not realistic.It is well known that in the case of a free surface as well as in the case of a clamped surface, one cannot satisfy the boundary conditions by the simple superposition of an incident wave and of a reflected wave of the same kind; one must add a transverse reflected wave if the incident is longitudinal and vice versa.The surface "scrambles" the waves so that one can no longer analyze the vibrations of the solid in terms of pure transverse and pure longitudinal modes." The results of Dupuis, Mazo, and Onsager for d = 3 were reproduced, for both the Dirichlet and the free boundary conditions, as rigorous theorems21 in [SaVa-97, Section 6.3], who also extended these results to the planar case d = 2.

Appendix B. A two-dimensional example: the disk
Our aim in this (and in the next) appendix to give an experimental verification of the correctness of the second asymptotic coefficients (1.27) and (1.28) and to demonstrate the incorrectness of the second asymptotic coefficient (1.32)-(1.33).We work with counting functions rather than with partition functions since the former can, in some cases, be explicitly22 calculated for reasonably large values of the parameter, whereas computing the latter would require additional trickery.
Let Ω ⊂ R d be the unit disk, equipped with standard polar coordinates (r, φ).The equations of elasticity (1.10) with Dirichlet boundary conditions in the disk allow the separation of variables23 .To this end, we, in essence, separate variables in the three-dimensional cylinder Ω×R following [MoFe-53, Chapter XIII] and looking for solutions independent of the third coordinate, cf. also Supplementary materials].More precisely, we take where z is the third coordinate vector.Then it is easily seen that the scalar potentials ψ j (r, φ), j = 1, 2, should satisfy the Helmholtz equations where 3) The general solution of (B.2) regular at the origin is well-known, Every solution Λ > 0 of the secular equation (B.5) is an eigenvalue of multiplicity one of the Dirichlet elasticity operator L Dir on the unit disk, and every solution Λ > 0 of the secular equation (B.6) is an eigenvalue of multiplicity two.Note that for the case of the disk the branching Hamiltonian billiards associated with the operator of elasticity can be analysed explicitly, and one can check that the two-term asymptotics (1.20) is valid.
The numerical results are shown in Figure 5.The free boundary problem for the disk is treated in the same manner, the results are shown in Figure 6.

Appendix C. Three-dimensional examples: flat cylinders
We consider24 Ω = Ω 3,h := T 2 [0, h], where T 2 is a flat square torus with side 2π and h > 0 is the height of the cylinder, so that Vol 3 (Ω) = (2π) 2 h and Vol 2 (∂Ω) = 8π 2 .We can once again separate variables by first setting where the ψ j are scalar potentials, and z is the coordinate vector in the direction of x 3 .Once again, it is easy to see that each potential ψ j satisfies (B.2), (B.3), with ω 3,Λ := ω 2,Λ .

(C.2)
The general solutions of (B.2) are now Substitution of (C.1), (B.3), (C.2), and (C.3) into the boundary conditions at x 3 = 0 and x 3 = h leads to some secular equations25 which, as it turns out, depend only on the values of rather than on the values of k 1 , k 2 themselves; we therefore only need to consider the values of K with Σ 2 (K ) > 0 where is the sum of squares function.Each solution Λ > 0 of a secular equation corresponding to such a K will be an eigenvalue of multiplicity Σ 2 (K ).
The results of our computations in the Dirichlet case are collated in Figure 7 and in the free boundary case in Figure 8.

Appendix D. Second Weyl coefficients in odd dimensions: proof of Theorem 1.10
This appendix is devoted to the proof of Theorem 1.10.We will prove (1.29) and (1.30) separately, in §D.1 and §D.2 respectively, by explicitly evaluating the integrals in the right-hand sides of (1.27) and (1.28).
In this appendix we denote complex variables by w = t + is, t , s ∈ R. §D.1.Dirichlet case: proof of (1.29).We begin by observing that, by performing a change of variable t = τ −2 , one obtains Therefore, proving (1.29) reduces to establishing the following Lemma D.1.For k = 1, 2, . . .we have Proof.The task at hand is to evaluate the integral Observing that the inverse tangent turns to zero at the endpoints of the interval of integration and integrating by parts, one obtains It is easy to see that the function ] with poles at w = 0 (of order k + 1), w = 1 + α −1 (of order 1) .
We choose the branch of the square root so that it is positive above the branch cut and negative below the branch cut.Note that the poles are not on the branch cut.Let Γ ϵ be a negatively oriented (clockwise) dog-bone contour around [1, α −1 ], see Figure 9. Then (D.2) Let C r be a positively oriented (counterclockwise) circular curve of radius r , see Figure 9, with r > 1 + α −1 .Then by Cauchy's Residue Theorem we have so that, combining (D.2) and (D.3), we obtain (D.8) Proof.The task at hand is to evaluate the integral In what follows, we assume, for simplicity, that α ̸ = 1 2 (this corresponds to λ ̸ = 0).The case α = 1 2 can be handled in a similar fashion.
Observing that the inverse tangent tends to π 2 at the endpoints of the interval of integration and integrating by parts, one obtains dt , (D.10) where R α is defined in accordance with (1.23).Hence, evaluating (D.9) reduces to evaluating where we choose the branch of the square root in such a way that on the upper side of the branch cut [0, α −1 ] we have (1 − αw) (w − 1) > 0. It is easy to see that the function f k is holomorphic in with poles at w = 0 (of order k + 1), w = w j , j = 1, 2, 3, (of order 1) , where the w j , j = 1, 2, 3, are the roots R α (w), with 0 < w 1 < 1 -recall the discussion from Remark 1.7.The nature of the other two roots w j , j = 2, 3, depends on α.Let α * ∈ (0, 1) be the unique real root of α 3 − 107 64 α 2 + 31 32 α − 11 64 = 0.
We will carry out the proof for the case (i); the other two cases are analogous, and lead to the same final result.Let Γ ϵ be a dog-bone contour as in Figure 10.Then Hence, since the function f k is regular at infinity, Cauchy's Residue Theorem gives us Res( f k , w j ) .

Figure 1 :
Figure 1: The appropriately rescaled coefficients b Dir (left image) and b free (right image) as functions of the parameter α for dimensions two to five.The inset in the right image shows the graph of γ R as a function of α.We note that lim α→0 + γ R ≈ 0.9553 and lim α→1 − γ R = 0.
, the ratio of the coefficient b Liu Dir and our coefficient b Dir .

Figure 2 :
Figure 2: The ratio b Liu Dir /b Dir for different dimensions, shown as functions of α.

Figure 3 :
Figure 3: The Dirichlet problem for the unit square.The second Weyl terms, both Liu's and ours, are compared to the actual numerically computed counting functions.In all three images we take µ = 1.

Figure 4 .
For further examples, both in the Dirichlet and the free boundary case, see Appendices B and C.

Figure 4 :
Figure 4: The free boundary problem for the unit square.The second Weyl terms are compared to the actual numerically computed counting functions.In all three images we take µ = 1.

Fact 5. 1 .
Let us denote by L plane the operator of linear elasticity for d = 2. Then the spectral shift function for the problem L ′ ∥ u ∥ = Λu ∥ , with Dirichlet/free boundary conditions coincides with the spectral shift function for the operator L ′ plane with the same boundary conditions.Namely, shift ℵ,∥ (Λ) = shift ℵ,plane (Λ), ℵ ∈ {Dir, free}.

Figure 5 :
Figure 5:The Dirichlet problem for the unit disk.In this case, (1.20) may be interpreted as N Dir (Λ)− aπΛ ≈ 2πb Dir Λ, and we compare the numerically computed left-hand sides with the right-hand sides computed using our and Liu's expressions for b Dir .In all three images we take µ = 1.

Figure 6 :
Figure 6: The free boundary counting function for the unit disk: in all three figures µ = 1.

Figure 7 :
Figure 7: The Dirichlet problem for flat cylinders.In all images µ = 1.

Figure 8 :
Figure 8: The free boundary problem for flat cylinders.In all images µ = 1.

Figure 9 :
Figure 9:The contour of integration for the Dirichlet problem.The two small circles centred at 1 and 1/α have radius ϵ.

Figure 10 :
Figure 10: The contour of integration for the free boundary problem.

Table 1 :
The coefficient b Dir for odd dimensions.