---
abstract: |
    We derive and analyze a hybridizable discontinuous Galerkin (HDG) method for approximating weak solutions to the equations of time-harmonic linear elasticity on a bounded Lipschitz domain in three dimensions. The real symmetry of the stress tensor is strongly enforced and its coefficients as well as those of the displacement vector field are approximated simultaneously at optimal convergence with respect to the choice of approximating spaces, wavenumber, and mesh size. Sufficient conditions are given so that the system is indeed transferable onto a global hybrid variable that, for larger polynomial degrees, may be approximated via a smaller-dimensional space than the original variables. We construct several variants of this method and discuss their advantages and disadvantages, and give a systematic approach to the error analysis for these methods. We touch briefly on the application of this error analysis to the time-dependent problem, and finally, we examine two different implementations of the method over various polynomial degrees and numerically demonstrate the convergence properties proven herein.\
    [**AMS Subject classification.**]{} 65N30, 65M50\
    [**Keywords.**]{} Hybridizable DG methods, elastic wave equation, time-harmonic solutions, optimal convergence.
author:
- |
    Allan Hungria${}^{(1)}$, Daniele Prada${}^{(2)}$, & Francisco–Javier Sayas${}^{(1)}$[^1]\
    (1) Dept of Mathematical Sciences, University of Delaware, USA\
    (2) Dept of Mathematical Sciences, Indiana University-Purdue University Indianapolis, USA\
    [allanh@udel.edu, dprada@umail.iu.edu, fjsayas@udel.edu ]{}
bibliography:
- 'HDGpaperBibliography.bib'
title: HDG methods for elastodynamics
---

Dedicated to Peter Monk on his 60th birthday

Introduction
============

We are concerned with numerical methods for the evolution of elastic waves on general (heterogeneous anisotropic) linearly elastic solids. It is well known that elastodynamics, in the time and frequency domains, has multiple applications in the fields of geophysics, material science, structural engineering, oil exploration, aerospace, etc. This work is a first contribution on the use of the Hybridizable Discontinuous Galerkin (HDG) method to the three-dimensional linear elastic wave equation in the time-harmonic regime, including some insights on energy conservation properties when these ideas are applied to evolutionary cases.

Mathematical literature contains a plethora of numerical methods for dealing with the elastic wave equation, each with its own virtue and applications: spectral elements [@Cohen:2001], particle-based methods such as the Hamiltonian Particle method (HPM) [@GoMaMiTa:2012], as well as the more finite element styled Continuous Galerkin (CG) methods [@Joly10:2008] and Discontinuous Galerkin (DG) methods [@CoNgPe:2011]. CG is well-known for its accuracy and reliability with smooth data and simple meshes. The DG framework is praised for its capacity to handle all sorts of complicated meshes and discontinuous data, to provide high-order accurate solutions, to perform $h/p$ adaptivity, and to retain very good scalability properties. However, DG methods have been criticized because, for the same mesh and the same polynomial degree, the number of globally coupled degrees of freedom is much larger than those of CG methods.

Certain DG methods, however, including the ones we shall explore here, have the key property of being *hybridizable*, i.e., the global system can be recast in terms of (statically condensed onto) a single “hybrid" variable that represents the trace of the solution on the boundaries of the elements. These form a family of methods that are, naturally, called the Hybridizable Discontinuous Galerkin (HDG) methods [@CoGoLa:2009]. The main idea for devising these methods is to: (i) use a characterization of the exact solution in terms of solutions of local problems and transmission conditions; (ii) use discontinuous approximations for both the solution inside each element and its trace on the element boundary; (iii) define the local solvers with the DG method; (iv) define a global problem by weakly imposing the transmission conditions. This creates a global linear system for only the hybrid variable, which is solved, after which the unknowns are recovered locally, again in parallel. This is similar to the hybridized implementation of mixed methods such as the Raviart-Thomas elements (see [@CoGoLa:2009], [@MoScSi:2010] and [@Sayas:2013] for more on this), except that the HDG method uses different polynomial spaces and a stabilization function instead of a stable mixed finite element pair. In certain cases, the hybrid space is smaller than that of the displacement/stress spaces [@KiMoShYa2015], and this has resulted in renewed interest in HDG.

The HDG methodology was successfully applied to time-harmonic acoustic waves by Roland Griesmaier and Peter Monk [@GrMo:2011]. Their analysis involves first rephrasing the classical system as first order in frequency before moving to the weak formulation. Testing the equations with the projected errors leads to a G[å]{}rding-type identity, and, combined with the dual equations to the classical system, the projected errors of both amplitude and its gradient can be bounded. This last bound requires a rather involved bootstrapping argument which is indispensable within our argument here.

Hybridizable DG methods have lately enjoyed further exposure in time-domain wave problems. For example, Cockburn and Quenneville-B[é]{}lair’s work on HDG for the acoustic wave equation [@CoQu:2014] provides much of the framework for the insights on the time-domain elastic problem at the end of this work. Nguyen, Peraire, and Cockburn implement an implicit HDG numerical scheme for both time-dependent acoustic and elastic equations [@CoNgPe:2011], and more recently, Stanglmeier, Nguyen, Peraire, and Cockburn explore an explicit scheme for the acoustic case [@CoNgPeSt:2016].

The vector field formulation of elasticity introduces several distinct complications in both the analysis and the implementation of HDG. Cockburn, Soon, and Stolarski give a numerical implementation of HDG for planar elasticity, along with a proof of existence and uniqueness of a solution to their particular HDG formulation [@CoSoSt:2009]. Cockburn, Fu, and Stolarski go on to analyze the convergence of this last method, which uses degree $k$ polynomial bases for displacement, stress, and hybrid spaces. They prove convergence at an order of $k+1$ for displacement and $k+1/2$ for the stress [@CoFuSt:2014], which is suboptimal; this has prompted the exploration of optimally convergent HDG methods.

One issue is that the tailored HDG projections often used in the analysis may not play well with the symmetry of the approximate stress tensor. Another is that using bases of the same polynomial degree for displacement, stress, and hybrid spaces leads to a suboptimal method. One method for addressing both of these issues is to introduce special divergence-free symmetric “bubble matrices" as in [@CoGoGu:2010], providing an extra control on the stress-associated approximation space. This yields a weakly-enforced symmetry of the approximate stress as well as optimal convergence of a postprocessed solution, taking advantage of some superconvergent quantities.

Another approach entirely is that of Weifeng Qiu, Jiguang Shen, and Ke Shi for the steady-state elasticity problem [@QiSh:2014]. The special tailored HDG projections are left behind for simpler $L^2$ projections, and the displacement-associated approximation space is expanded by one polynomial degree. While this does then require some extra terms to be bounded in the analysis, the net result is shown to achieve optimal convergence directly. (Note that due to the disparity of polynomial degrees for the stress and displacement, optimal convergence of this method yields the same quality of the solution as a postprocessed method based on a superconvergent scheme.) An important feature of this approach is the strong symmetry of the approximate stress. See the introduction of [@QiSh:2014] for more on this. Expanding the polynomial degree of the primal unknown by one is an idea that can be traced back to Lehrenfeld and Schoberl [@LeSc:2015], but Qiu, Shen and Shi compensate by adjusting the order of the primal unknown piece of the stabilization function to $O(h^{-1})$ as well as a projection operator from primal approximation space onto hybrid space.

Our choice of polynomial approximating spaces and projections is that of Qiu, Shen, and Shi in order to be able to work on the most general polyhedral mesh possible. However, the frequency-domain problem, unlike the steady-state problem, is not coercive, so we wind up with a Gårding-type identity similar to that of Griesmaier and Monk’s [@GrMo:2011], after following their example and first phrasing the classical system as first order in both frequency and space. The two analytical recipes from [@QiSh:2014] and [@GrMo:2011] are here carefully blended to approach the time-harmonic elasticity case, which has implications on the choice of numerical flux and its dependence on the wavenumber.

The following treatment of HDG for time-harmonic elasticity, however, comes with its own complications, not only with regard to the hybridization of the DG scheme, but also in consideration of the dependence on the wavenumber. We have also developed a simplified system for dealing with the double-bootstrapping process, which is now even messier due to the use of $L^2$ projections rather than tailored HDG projections. By varying the numerical flux, we wind up with several different HDG methods for the time-harmonic linear elastic problem. We proceed to show how some of these methods can be used to produce semi-discretizations in the time domain and that one of them is actually conservative.

What follows is a rigorous treatment of the error analysis and well-posedness of HDG methodology as applied to the problem of three-dimensional time-harmonic elasticity on a polyhedron with mixed boundary conditions and strong symmetric stresses. We explore how this analysis can shape the stability mechanism for a method of numerically integrating the time-dependent system, in particular for the second-order-in-frequency case. Numerical experiments are carried out to demonstrate convergence of both the first-order method and a second-order variant. We then compare, using various polynomial degrees and tetrahedrizations, the sizes of the global linear systems involved in HDG and Lagrange element CG, demonstrating an advantage of HDG at large polynomial degrees.

#### Notation.

Given an open set $\Omega$, we will write $$(u,v)_\Omega:=\int_\Omega u\, v,
    \quad
(\mathbf u,\mathbf v)_\Omega:=\int_\Omega \mathbf u\cdot\mathbf v,
    \quad
({\boldsymbol}\xi,{\boldsymbol}\chi)_\Omega:=\int_\Omega {\boldsymbol}\xi:{\boldsymbol}\chi:=\int_\Omega \mathrm{trace}({\boldsymbol}\chi^\top {\boldsymbol}\xi),$$ for real square-integrable scalar, vector-valued, and matrix-valued functions. The symbol $\top$ will be used for real transposition of matrices. When used for complex-valued fields, all brackets will still be defined in the same way, and will be therefore bilinear and not sesquilinear. In the same spirit, $\top$ will denote transposition and the colon will be defined as above, even when applied to complex matrices. The set of symmetric real $3\times 3$ matrices will be denoted ${{\mathbb R}}^{3\times 3}_{\mathrm{sym}}$ and the set of symmetric (not Hermitian) complex $3\times 3$ matrices will be denoted ${{\mathbb C}}^{3\times 3}_{\mathrm{sym}}$. Similarly, the set of antisymmetric real $3\times 3$ matrices will be denoted ${{\mathbb R}}^{3\times 3}_{\mathrm{skw}}$ and the set of antisymmetric complex $3\times 3$ matrices will be denoted ${{\mathbb C}}^{3\times 3}_{\mathrm{skw}}$.

Problem setting and HDG discretization {#sec:2}
======================================

Let $\Omega\subset \mathbb R^3$ be a polyhedron with Lipschitz boundary $\Gamma$. We assume that $\Gamma$ is divided into two non-overlapping parts $\Gamma_D$ and $\Gamma_N$, where we will respectively impose Dirichlet and Neumann boundary conditions. For simplicity we will assume that each of $\Gamma_D$ and $\Gamma_N$ are made up of full faces of $\Gamma$. We will be looking for a displacement field $\mathbf u:\Omega\to {{\mathbb C}}^3$ and for the associated stress tensor $\widetilde{{\boldsymbol}\sigma}:\Omega\to {{\mathbb C}}^{3\times 3}_{\mathrm{sym}}$. The stress field is given by a general linear non-homogeneous anisotropic law: $$\widetilde{{\boldsymbol}\sigma}=\mathcal C{\boldsymbol}\varepsilon(\mathbf u),
\qquad{\boldsymbol}\varepsilon(\mathbf u):=\tfrac12 (\nabla\mathbf u+(\nabla \mathbf u)^\top),$$ where for almost every $\mathbf x\in \Omega$, the linear operator $\mathcal C(\mathbf x)$ transforms real symmetric matrices into real symmetric matrices, satisfies the symmetry condition $$(\mathcal C(\mathbf x){\boldsymbol}\xi):{\boldsymbol}\chi 
={\boldsymbol}\xi : (\mathcal C(\mathbf x){\boldsymbol}\chi) \quad \forall {\boldsymbol}\xi,{\boldsymbol}\chi\in \mathbb R^{3\times 3}_{\mathrm{sym}}$$ and the positivity condition $$(\mathcal C(\mathbf x) {\boldsymbol}\xi):{\boldsymbol}\xi
\ge C_0\,  {\boldsymbol}\xi:{\boldsymbol}\xi \qquad \forall {\boldsymbol}\xi\in {{\mathbb R}}^{3\times 3}_{\mathrm{sym}},$$ for some $C_0>0$. Moreover, we assume that the components of the tensor $\mathcal C$ with respect to the canonical basis of $\mathbb R^{3\times 3}_{\mathrm{sym}}$ are $L^\infty(\Omega)$ functions. The other physical parameter in the equations to follow is the strictly positive bounded density $\rho:\Omega \to {{\mathbb R}}$, so that the weighted norm is equivalent $$C_1 \|\mathbf u\|_\Omega^2 \le \|\mathbf u\|_\rho^2:=(\rho\mathbf u,\overline{\mathbf u})_\Omega
\le C_2 \|\mathbf u\|_\Omega^2.$$ In strong primal form, our problem is the search of $\mathbf u$ such that

\[eq:2.1\] $$\begin{aligned}
{6}
\nabla\cdot\widetilde{{\boldsymbol}\sigma} + \kappa^2\rho\,\mathbf u
     &= \widetilde{\mathbf f} &\qquad &\mbox{in $\Omega$},\\
\mathbf u &=\mathbf g_D & & \mbox{on $\Gamma_D$},\\
\widetilde{{\boldsymbol}\sigma}\,\mathbf n &=\widetilde{\mathbf g}_N & & \mbox{on $\Gamma_N$},\end{aligned}$$

where: $\widetilde{{\boldsymbol}\sigma}=\mathcal C{\boldsymbol}\varepsilon(\mathbf u)$, the divergence operator is applied to the rows of $\widetilde{{\boldsymbol}\sigma}$, $\widetilde{\mathbf f}\in \mathbf L^2(\Omega)=L^2(\Omega)^3$, $\mathbf g_D\in \mathbf H^{1/2}(\Gamma_D)=H^{1/2}(\Gamma_D)^3$, $\widetilde{\mathbf g}_N \in \mathbf L^2(\Gamma_N)=L^2(\Gamma_N)^3$, $\mathbf n$ is the unit outward-pointing normal vector field on $\Gamma_N$, and $\kappa > 0$ is the wave number. In the weak primal formulation (where problem is typically studied), the Dirichlet condition is understood in the sense of traces, while the Neumann condition holds in a dual space with negative Sobolev index on the boundary $\Gamma_N$. We will assume that $\kappa^2$ is not an eigenvalue for the associated Navier-Lamé operator $\mathbf u\mapsto -\nabla\cdot (\mathcal C{\boldsymbol}\varepsilon(\mathbf u))$ with the given boundary conditions, i.e., we assume that the only solution of with zero right-hand side is the trivial solution.

The discretization will be done for a first order (in space and frequency) reformulation of . We first need to invert the elastic law $\mathcal C$. With the hypothesis given for $\mathcal C$, we can assert that for almost every $\mathbf x\in \Omega$, there exists a linear operator $\mathcal A(\mathbf x)=\mathcal C(\mathbf x)^{-1}$, transforming real symmetric matrices into real symmetric matrices. On the set of matrix-valued ${\boldsymbol}\xi:\Omega \to \mathbb C^{3\times 3}_{\mathrm{sym}}$ functions with $L^2(\Omega)$ components, we define the elastic potential norm $$C_1\| {\boldsymbol}\xi\|_\Omega^2 \le \| {\boldsymbol}\xi\|_{\mathcal A}^2 :=
(\mathcal A {\boldsymbol}\xi,\overline{{\boldsymbol}\xi})_{\Omega}=(\mathcal A{\boldsymbol}\xi_{\mathrm{re}}+\imath \mathcal A{\boldsymbol}\xi_{\mathrm{im}},\overline{{\boldsymbol}\xi})_\Omega
\le C_2 \| {\boldsymbol}\xi\|_\Omega^2.$$ We emphasize that we will always use bilinear (not sesquilinear) brackets for $L^2$-type products, and that the symbol used for the Frobenius product of matrices (the colon) will not include conjugation. We then introduce the new unknown and data $${\boldsymbol}\sigma:=\frac\imath\kappa \widetilde{{\boldsymbol}\sigma}=\frac\imath\kappa \mathcal C {\boldsymbol}\varepsilon(\mathbf u),
    \qquad
\mathbf f:=\frac\imath\kappa \widetilde{\mathbf f},
    \qquad
\mathbf g_N:=\frac\imath\kappa \widetilde{\mathbf g}_N,$$ and write as the equivalent first order system

\[eq:2.2\] $$\begin{aligned}
{6}
\imath\kappa\mathcal A{\boldsymbol}\sigma +{\boldsymbol}\varepsilon(\mathbf u) &=\mathbf 0
    &\qquad &\mbox{in $\Omega$},\\
\nabla\cdot{\boldsymbol}\sigma +\imath\kappa\,\rho\,\mathbf u &=\mathbf f
    &\qquad &\mbox{in $\Omega$},\\
\mathbf u &=\mathbf g_D
    &\qquad &\mbox{on $\Gamma_D$},\\
{\boldsymbol}\sigma\mathbf n &=\mathbf g_N
    &\qquad &\mbox{on $\Gamma_N$}.\end{aligned}$$

A similar formulation can be found using the original data and the stress tensor $\widetilde{\bf\sigma}$, so that the equations are $\mathcal A\widetilde{{\boldsymbol}\sigma}-{\boldsymbol}\varepsilon(\mathbf u)=\mathbf 0$ and $\nabla\cdot\widetilde{{\boldsymbol}\sigma}+\kappa^2 \,\rho\,\mathbf u=\widetilde{\mathbf f}$. (This will be discussed in Section \[sec:6\].)

We next introduce the HDG discretization of . Since the method we use is Qiu & Shi’s [@QiSh:2014], we will not repeat the derivation. We start with a shape-regular conforming tetrahedrization ${{\mathcal T_h}}$ of the domain $\Omega$. The set of all faces of elements of ${{\mathcal T_h}}$ is denoted $\mathcal E_h$, and we will sometimes understand that $\mathcal E_h$ is the geometric skeleton of the triangulation, i.e., the union of all the faces of all elements. The method involves three discrete spaces

\[eq:2.3\] $$\begin{aligned}
{{\mathbb V_h}}&:=& \{ {\boldsymbol}\xi:\Omega\to {{\mathbb C}}^{3\times 3}_{\mathrm{sym}}\,:\, {\boldsymbol}\xi|_K \in \mathcal P_k(K;{{\mathbb C}}^{3\times 3}_{\mathrm{sym}}) \quad \forall K\in {{\mathcal T_h}}\},\\
{{\mathbf W_h}}&:=& \{ \mathbf u:\Omega\to {{\mathbb C}}^3\,:\, \mathbf u|_K \in \mathcal P_{k+1}(K;{{\mathbb C}}^3)\quad
\forall K\in {{\mathcal T_h}}\},\\
{{\mathbf M_h}}&:=& \{ {\boldsymbol}\mu:\mathcal E_h \to {{\mathbb C}}^3\,:\, {\boldsymbol}\mu|_F\in \mathcal P_k(F;{{\mathbb C}}^3)\quad
\forall K\in \mathcal E_h\}.\end{aligned}$$

In , $\mathcal P_r(K;S)$ is the set of polynomials of total degree up to $r$ defined on $K$ and with values in $S\in \{{{\mathbb C}}^{3\times 3}_{\mathrm{sym}},{{\mathbb C}}^3\}$, while $\mathcal P_k(F;{{\mathbb C}}^3)$ are vector valued polynomials on the tangential coordinates defined on the face $F$ and of degree not greater than $k$. We will also use the orthogonal projector $$\label{eq:2.4new}
{{\mathbf P_M}}: \prod_{K\in {{\mathcal T_h}}} L^2(\partial K) \longrightarrow \prod_{K\in {{\mathcal T_h}}} \prod_{F\in \mathcal E(K)} \mathcal P_k(F;{{\mathbb C}}^3),$$ where $\mathcal E(K)$ is the set of faces of $\partial K$. Note that ${{\mathbf M_h}}$ can be identified with the subspace of the set of the right-hand side of consisting of functions that are single-valued on internal faces.

Stabilization is carried out through a function ${\boldsymbol}\tau$ defined as follows: for each element $K\in {{\mathcal T_h}}$, a function ${\boldsymbol}\tau_K:\partial K \to \mathbb R^{3\times 3}_{\mathrm{sym}}$ satisfying (a) ${\boldsymbol}\tau_K|_{F}$ is constant on each $F\in \mathcal E(K)$; (b) there exist two positive constants such that $$\label{eq:tau}
C_1 h_K^{-1} \|{\boldsymbol}\mu\|_{\partial K}^2 \le \langle{\boldsymbol}\tau_K {\boldsymbol}\mu,\overline{{\boldsymbol}\mu}\rangle_{\partial K}
    \le C_2 h_K^{-1} \|{\boldsymbol}\mu\|_{\partial K}^2
    \quad \forall {\boldsymbol}\mu\in \mathbf L^2(\partial K), \quad \forall K\in {{\mathcal T_h}},$$ where $h_K$ is the diameter of $K$. The symbol ${\boldsymbol}\tau$ will be used to denote the function defined on the set of boundaries of all elements as above, understanding that ${\boldsymbol}\tau$ can be double-valued on interior faces. The numerical fluxes follow the pattern of HDG methods: the one corresponding to the displacement will be an unknown $\widehat{\mathbf u}_h \in {{\mathbf M_h}}$, while the one related to the (normal) stress is given elementwise with a formula in terms of all the unknowns $$\label{eq:2.4}
\widehat{{\boldsymbol}\sigma}_h\mathbf n:={\boldsymbol}\sigma_h\mathbf n+{\boldsymbol}\tau_K ({{\mathbf P_M}}\mathbf u_h-\widehat{\mathbf u}_h) 
\,:\, \partial K\to {{\mathbb C}}^3.$$ Here the normal vector field $\mathbf n:\partial K\to {{\mathbb R}}^3$ is unitary and points to the exterior of $K$. At this time, we can write the HDG discrete equations for . We look for $({\boldsymbol}\sigma_h,\mathbf u_h,\widehat{\mathbf u}_h)\in {{\mathbb V_h}}\times {{\mathbf W_h}}\times {{\mathbf M_h}}$ satisfying

\[eq:2.5\] $$\begin{aligned}
{6}
\label{eq:2.5a}
\imath\kappa (\mathcal A{\boldsymbol}\sigma_h,{\boldsymbol}\xi)_{{{\mathcal T_h}}}
    -(\mathbf u_h,\nabla\cdot{\boldsymbol}\xi)_{{\mathcal T_h}}+\langle\widehat{\mathbf u}_h,{\boldsymbol}\xi\mathbf n\rangle_{\partial{{\mathcal T_h}}}
     & = 0 & \qquad & \forall {\boldsymbol}\xi\in {{\mathbb V_h}},\\
\label{eq:2.5b} 
-({\boldsymbol}\sigma_h,\nabla\mathbf w)_{{{\mathcal T_h}}}
    +\langle \widehat{{\boldsymbol}\sigma}_h\mathbf n,\mathbf w\rangle_{\partial{{\mathcal T_h}}}
    +\imath\kappa (\rho\,\mathbf u_h,\mathbf w)_{{\mathcal T_h}}&=(\mathbf f,\mathbf w)_{{\mathcal T_h}}& & \forall \mathbf w\in {{\mathbf W_h}},\\
\label{eq:2.5c}
\langle\widehat{{\boldsymbol}\sigma}_h\mathbf n,{\boldsymbol}\mu\rangle_{\partial{{\mathcal T_h}}\setminus\Gamma_D}
    &=\langle\mathbf g_N,{\boldsymbol}\mu\rangle_{\Gamma_N} & & \forall {\boldsymbol}\mu \in {{\mathbf M_h}},\\
\label{eq:2.5d}
\langle\widehat{\mathbf u}_h,{\boldsymbol}\mu\rangle_{\Gamma_D} &=\langle\mathbf g_D,{\boldsymbol}\mu\rangle_{\Gamma_D}
    & &\forall {\boldsymbol}\mu\in {{\mathbf M_h}},\end{aligned}$$

with as the definition of $\widehat{{\boldsymbol}\sigma}_h\mathbf n$ and brackets defined as follows: $$(\mathbf u,\mathbf v)_{{\mathcal T_h}}:=\sum_{K\in {{\mathcal T_h}}} (\mathbf u,\mathbf v)_K, \qquad
\langle \mathbf u,\mathbf v\rangle_{\partial{{\mathcal T_h}}}:=\sum_{K\in {{\mathcal T_h}}} \langle \mathbf u,\mathbf v\rangle_{\partial K}
    :=\sum_{K\in {{\mathcal T_h}}}\int_{\partial K}\mathbf u\cdot\mathbf v,$$ and $$\langle \mathbf u,\mathbf v\rangle_{\partial{{\mathcal T_h}}\setminus\Gamma_D}
:=\sum_{K\in {{\mathcal T_h}}} \langle \mathbf u,\mathbf v\rangle_{\partial K\setminus\Gamma_D}.$$ Equations and can be added together as a single equation tested against ${{\mathbf M_h}}$, which shows that is a square system of linear equations. The discrete momentum equation can be equivalently written as $$\label{eq:2.5bnew}
(\nabla\cdot{\boldsymbol}\sigma_h,\mathbf w)_{{\mathcal T_h}}+\langle{\boldsymbol}\tau ({{\mathbf P_M}}\mathbf u_h-\widehat{\mathbf u}_h),{{\mathbf P_M}}\mathbf w\rangle_{\partial{{\mathcal T_h}}}
    +\imath\kappa (\rho\,\mathbf u_h,\mathbf w)_{{\mathcal T_h}}=(\mathbf f,\mathbf w)_{{\mathcal T_h}}\quad \forall\mathbf w\in {{\mathbf W_h}}.$$ We note that the degree of the polynomial space used for $\mathbf u_h$ is one higher than the one used for the other unknowns and the fact that ${{\mathbf P_M}}$ has been introduced in the definition of the flux so that $\widehat{{\boldsymbol}\sigma}_h\mathbf n \in 
\prod_{F\in \mathcal E(K)} \mathcal P_k(F;{{\mathbb C}}^3).$

Main results {#sec:3}
============

#### Regularity assumptions.

From now on we will assume that $\rho$ and the coefficients of $\mathcal C$ are in $W^{1,\infty}({{\mathcal T_h}})$. Let us now consider the coercive problem

\[eq:3.1\] $$\begin{aligned}
{6}
\nabla\cdot(\mathcal C {\boldsymbol}\varepsilon(\mathbf w))-\rho\,\mathbf w
     &= \mathbf r &\qquad &\mbox{in $\Omega$},\\
\mathbf w &={\boldsymbol}0 & & \mbox{on $\Gamma_D$},\\
(\mathcal C {\boldsymbol}\varepsilon(\mathbf w))\,\mathbf n 
&={\boldsymbol}0 & & \mbox{on $\Gamma_N$}.\end{aligned}$$

We will also assume that the solution of for arbitrary $\mathbf r\in L^2(\Omega;\mathbb R^3)$ is in $H^2(\Omega;\mathbb R^3)$ and that there exists a constant $C>0$ such that $$\label{eq:3.2}
\| \mathbf w\|_{2,\Omega} \le C \|\mathbf r\|_\Omega.$$ For the time-harmonic problem, we will denote by $C_\kappa>0$ the constant such that the solution of

\[eq:3.3\] $$\begin{aligned}
{6}
\nabla\cdot(\mathcal C {\boldsymbol}\varepsilon(\mathbf w))+\kappa^2\rho\,\mathbf w
     &= \mathbf r &\qquad &\mbox{in $\Omega$},\\
\mathbf w &={\boldsymbol}0 & & \mbox{on $\Gamma_D$},\\
(\mathcal C {\boldsymbol}\varepsilon(\mathbf w))\,\mathbf n 
&={\boldsymbol}0 & & \mbox{on $\Gamma_N$}.\end{aligned}$$

can be bounded by $$\label{eq:3.4}
\| \mathbf w\|_{1,\Omega}\le C_\kappa \|\mathbf r\|_\Omega.$$ Note that we have assumed the unique solvability of .

#### Error quantities.

The error analysis will be carried out by comparing numerical solutions and orthogonal projections. Let ${{\Pi_V}}: L^2(\Omega;{{\mathbb C}}^{3\times 3}_{\mathrm{sym}})\to {{\mathbb V_h}}$ and ${{\Pi_W}}:L^2(\Omega;{{\mathbb C}}^3)\to{{\mathbf W_h}}$ be the orthogonal projections onto the discrete spaces. Consider the errors $${{e_h^\sigma}}:={{\Pi_V}}{\boldsymbol}\sigma-{\boldsymbol}\sigma_h,
    \qquad
{{\mathbf e_h^u}}:={{\Pi_W}}\mathbf u-\mathbf u_h,
    \qquad
{{\widehat{\mathbf e}_h^u}}:={{\mathbf P_M}}\mathbf u-\widehat{\mathbf u}_h,$$ and the best approximation errors $${{\varepsilon_h^\sigma}}:={{\Pi_V}}{\boldsymbol}\sigma-{\boldsymbol}\sigma,
    \qquad
{{{\boldsymbol}\varepsilon_h^u}}:={{\Pi_W}}\mathbf u-\mathbf u.$$ For convenience, we introduce the skeleton norm $$\| {\boldsymbol}\mu\|_\tau:=\langle {\boldsymbol}\tau {\boldsymbol}\mu,\overline{{\boldsymbol}\mu}\rangle_{\partial{{\mathcal T_h}}}^{1/2}.$$ For the error analysis we will allow constants to depend on the density $\rho$ and on the coefficients of $\mathcal A$. While the influence of these physical coefficients in the inequalities can be tracked with careful arguments, the results seem to be too involved to obtain precise conclusions on how $h$ and $\kappa$ interact with them. However, we will pay attention to the maximum spectral value of the inverse compliance tensor, i.e., to the positive bounded function such that for almost every $\mathbf x\in \Omega$ $$\label{eq:5.3new}
(\mathcal A(\mathbf x){\boldsymbol}\xi):{\boldsymbol}\xi\le c_{\mathcal A}(\mathbf x) \, 
     {\boldsymbol}\xi :{\boldsymbol}\xi  \qquad \forall {\boldsymbol}\xi \in {{\mathbb R}}^{3\times 3}_{\mathrm{sym}}.$$

\[the:3.1\] There exist $C_1, C_2>0$, dependent only on the shape-regularity of ${{\mathcal T_h}}$, the density $\rho$ and the coefficients of the inverse compliance tensor $\mathcal A$ such that if $h (1+\kappa)^{3/2}(1+\kappa C_\kappa+C_\kappa)$ is small enough, then the errors can be bounded by $$\| {{e_h^\sigma}}\|_{\mathcal A}+\kappa^{-1/2} \| {{\mathbf P_M}}{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}}\|_\tau
    \le C_1 (1+\kappa^{-1/2}) \big( h^t |{\boldsymbol}\sigma|_{t,\Omega}+h^{s-1} |\mathbf u|_{s,\Omega}\big)$$ and $$\| {{\mathbf e_h^u}}\|_\Omega\le C_2 (1+\kappa C_\kappa) \kappa^{-1/2}(1+\kappa)^2 
\big( h^{t+1} |{\boldsymbol}\sigma|_{t,\Omega}+h^{s} |\mathbf u|_{s,\Omega}\big),$$ if $k\ge 1$, $\mathbf u\in H^s(\Omega;{{\mathbb C}}^3)$ with $1\le s\le k+2$, and ${\boldsymbol}\sigma\in H^t(\Omega;{{\mathbb C}}^{3\times 3})$ with $1\le t\le k+1$.

Optimal error estimates are $$\| {\boldsymbol}\sigma-{\boldsymbol}\sigma_h\|_\Omega=\mathcal O(h^{k+1}), 
    \qquad
\|  \mathbf u-\mathbf u_h\|_\Omega=\mathcal O(h^{k+2}).$$ With some additional scaling inequalities, keeping in mind that ${\boldsymbol}\tau$ scales like $h^{-1}$ elementwise, it is possible to show that $$\| \mathbf u-\widehat{\mathbf u}_h\|_\tau =\mathcal O(h^{k+1}).$$ The estimates of Theorem \[the:3.1\] can also be written in terms of the original physical variables. If we denote $\widetilde{{\boldsymbol}\sigma}_h:=-\imath\kappa{\boldsymbol}\sigma_h$, then $$\begin{aligned}
{6}
\| {{\Pi_V}}\widetilde{{\boldsymbol}\sigma}-\widetilde{{\boldsymbol}\sigma}_h\|_\Omega+
    \kappa^{1/2} \| {{\mathbf P_M}}{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}}\|_\tau 
    & \le
    C_1 (1+\kappa^{-1/2}) \big( h^t |\widetilde{{\boldsymbol}\sigma}|_{t,\Omega}+h^{s-1}\kappa |\mathbf u|_{s,\Omega}\big),\\
\| {{\mathbf e_h^u}}\|_\Omega
    &\le C_2 (1+\kappa C_\kappa) \kappa^{-3/2}(1+\kappa)^2 
       \big( h^{t+1} |\widetilde{{\boldsymbol}\sigma}|_{t,\Omega}+h^{s} \kappa|\mathbf u|_{s,\Omega}\big).\end{aligned}$$

#### Unique solvability.

Theorem \[the:3.1\] can be used to prove existence and uniqueness of solution of for $h$ small enough (depending on the wave number $\kappa$). The argument is as follows. Consider the system with homogeneous data: $\mathbf f={\boldsymbol}0$, $\mathbf g_N={\boldsymbol}0$, and $\mathbf g_D={\boldsymbol}0$. Let $({\boldsymbol}\sigma_h,\mathbf u_h,\widehat{\mathbf u}_h)$ be any solution of this homogeneous set of linear equations. Theorem \[the:3.1\] applied to this solution and the exact zero solution shows that $({\boldsymbol}\sigma_h,\mathbf u_h,\widehat{\mathbf u}_h)$ has to vanish. Therefore, the linear system (with as many equations as unknowns) is uniquely solvable for any right-hand side. The logic of the use of Theorem \[the:3.1\] is slightly warped: it assumes the existence of a discrete solution, which we know to happen at least for the homogeneous case, and then it uses the error estimates to show that the system is actually uniquely solvable.

Local solvability and energy identity
=====================================

\[lemma:5.1\] There exists $C>0$, depending only on the shape regularity of the grid, such that $$\|\mathbf v\|_K \le C h_K \|{\boldsymbol}\varepsilon(\mathbf v)\|_K$$ for all $\mathbf v\in H^1(K;{{\mathbb C}}^3)$ satisfying $$\label{eq:5.1}
\langle \mathbf v,{\boldsymbol}\mu\rangle_{\partial K}=0 \qquad \forall {\boldsymbol}\mu\in \prod_{F\in \mathcal E(K)} \mathcal P_0(F;{{\mathbb C}}^3).$$

A scaling argument, using only that $\int_{\partial K}\mathbf v=\mathbf 0$ and a Poincaré inequality on the reference element prove that $$\|\mathbf v\|_K \le C h_K \|\nabla\mathbf v\|_K.$$ On the other hand, by a straightforward extension of [@QiSh:2014 Lemma 4.1] to our complex-valued fields, we have the local Korn inequality $$\label{eq:5.2}
\inf_{{\boldsymbol}\xi\in \mathbb C^{3\times 3}_{\mathrm{skw}}}\| \nabla\mathbf v+{\boldsymbol}\xi\|_K
    \le C \|{\boldsymbol}\varepsilon(\mathbf v)\|_K \quad\forall\mathbf v\in H^1(K;{{\mathbb C}}^3).$$ The constant in depends only on the shape-regularity constant of the mesh. Finally, if ${\boldsymbol}\xi\in {{\mathbb C}}^{3\times 3}$, then $$(\nabla\mathbf v,{\boldsymbol}\xi)_K=\langle \mathbf v,{\boldsymbol}\xi\mathbf n\rangle_{\partial K}=0,$$ since $\mathbf v$ satisfies . Therefore $$\inf_{{\boldsymbol}\xi\in \mathbb C^{3\times 3}_{\mathrm{skw}}}\| \nabla\mathbf v+{\boldsymbol}\xi\|_K=\|\nabla\mathbf v\|_K$$ and the proof is finished.

The following result shows that the local equations associated to - are uniquely solvable, i.e., given the data functions and $\widehat{\mathbf u}_h$, we can compute ${\boldsymbol}\sigma_h$ and $\mathbf u_h$ element by element. This is the key ingredient to show that the HDG method is actually hybridizable, that is, it can be recast as a linear system where $\widehat{\mathbf u}_h\in {{\mathbf M_h}}$ is the only variable. To simplify the proof, we introduce the weighted norms $$\| {\boldsymbol}\xi\|_{\mathcal A,K}^2:=(\mathcal A{\boldsymbol}\xi,\overline{{\boldsymbol}\xi})_K,
    \qquad
\| \mathbf v\|_{\rho,K}^2:=(\rho\,\mathbf v,\overline{\mathbf v})_K.$$

\[prop:5.2\] If $C>0$ is the constant of Lemma \[lemma:5.1\] and $$\label{eq:5.4new}
\kappa\, h_K < \frac1{C\,\| c_{\mathcal A}\|_{L^\infty(K)}^{1/2} \| \rho\|_{L^\infty(K)}^{1/2}},$$ then the local solver associated to the element $K\in {{\mathcal T_h}}$ is well defined. In other words, if $({\boldsymbol}\sigma,\mathbf u)\in \mathcal P_k(K;{{\mathbb C}}^{3\times3}_{\mathrm{sym}})\times \mathcal P_{k+1}(K;{{\mathbb C}}^3)$ satisfies

$$\begin{aligned}
{6}
\label{eq:5.3a}
\imath\kappa (\mathcal A{\boldsymbol}\sigma,{\boldsymbol}\xi)_K -(\mathbf u,\nabla\cdot{\boldsymbol}\xi)_K 
    &=0
    &\qquad &\forall {\boldsymbol}\xi\in \mathcal P_k(K;{{\mathbb C}}^{3\times3}_{\mathrm{sym}}),\\
\label{eq:5.3b}
(\nabla\cdot{\boldsymbol}\sigma,\mathbf w)_K+\langle{\boldsymbol}\tau{{\mathbf P_M}}\mathbf u,\mathbf w\rangle_{\partial K}
    +\imath\kappa(\rho\mathbf u,\mathbf w)_K
    &=0
    &\qquad & \forall\mathbf w\in \mathcal P_{k+1}(K;{{\mathbb C}}^3),
    \end{aligned}$$

then $({\boldsymbol}\sigma,\mathbf u)=(\mathbf 0,\mathbf 0)$.

Note that we only need to prove that $\mathbf u=\mathbf 0$. Testing with $\overline{{\boldsymbol}\sigma}$, conjugating and testing it with $\overline{\mathbf u}$, and adding the result of these two equations, it follows that $$\imath\kappa\left(\|{\boldsymbol}\sigma\|_{\mathcal A,K}^2-\|\mathbf u\|_{\rho,K}^2\right)
    +\langle{\boldsymbol}\tau{{\mathbf P_M}}\overline{\mathbf u},{{\mathbf P_M}}\mathbf u\rangle_{\partial K}=0.$$ By it follows that ${{\mathbf P_M}}\mathbf u=\mathbf 0$ and $\|{\boldsymbol}\sigma\|_{\mathcal A,K}=\|\mathbf u\|_{\rho,K}$. Going back to , integrating by parts, and using that ${{\mathbf P_M}}\mathbf u=\mathbf 0$, it follows that $$\label{eq:5.4}
\imath\kappa (\mathcal A{\boldsymbol}\sigma,{\boldsymbol}\xi)_K+(\nabla\mathbf u,{\boldsymbol}\xi)_K=0 \quad \forall  {\boldsymbol}\xi\in \mathcal P_k(K;{{\mathbb C}}^{3\times3}_{\mathrm{sym}}).$$ Testing with ${\boldsymbol}\xi={\boldsymbol}\varepsilon(\overline{\mathbf u})$, it follows that $$\| {\boldsymbol}\varepsilon( \mathbf u)\|_K^2
    = (\nabla\mathbf u,{\boldsymbol}\varepsilon(\overline{\mathbf u}))_K
            = \kappa | (\mathcal A{\boldsymbol}\sigma,{\boldsymbol}\varepsilon(\overline{\mathbf u}))_K|  \le 
            \kappa \| c_{\mathcal A}\|_{L^\infty(K)}^{1/2} \|{\boldsymbol}\sigma\|_{\mathcal A,K} \| {\boldsymbol}\varepsilon( \mathbf u)\|_K,$$ where we have used . Note that $\mathbf u$ satisfies , given the fact that ${{\mathbf P_M}}\mathbf u=\mathbf 0$. Therefore, by Lemma \[lemma:5.1\], if ${\boldsymbol}\varepsilon(\mathbf u)=\mathbf 0$, then $\mathbf u=\mathbf 0$ and the proof is finished. Otherwise $\mathbf u\neq \mathbf 0$ and, by Lemma \[lemma:5.1\] and the equality $\|{\boldsymbol}\sigma\|_{\mathcal A,K}=\|\mathbf u\|_{\rho,K}$, we can bound $$\begin{aligned}
\| \mathbf u\|_K & \le & C h_K \| {\boldsymbol}\varepsilon(\mathbf u)\|_K 
            \le C \kappa\, h_K \| c_{\mathcal A}\|_{L^\infty(K)}^{1/2} \|{\boldsymbol}\sigma\|_{\mathcal A,K}
            = C \kappa\, h_K \| c_{\mathcal A}\|_{L^\infty(K)}^{1/2} \|\mathbf u\|_{\rho,K} \\
            &\le & C \kappa\, h_K \| c_{\mathcal A}\|_{L^\infty(K)}^{1/2} \| \rho\|_{L^\infty(K)}^{1/2}
                    \|\mathbf u\|_K\end{aligned}$$ and we arrive at a contradiction if holds.

\[prop:5.3\] $$\begin{aligned}
{6}
& \imath\kappa (\|{{e_h^\sigma}}\|_{\mathcal A}^2-\|{{\mathbf e_h^u}}\|_\rho^2)+\|{{\mathbf P_M}}{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}}\|_\tau^2 \\
& \hspace{5pt} =\imath\kappa\left( (\mathcal A{{\varepsilon_h^\sigma}},\overline{{e_h^\sigma}})_{{\mathcal T_h}}-(\rho \overline{{{\boldsymbol}\varepsilon_h^u}},{{\mathbf e_h^u}})_{{\mathcal T_h}}\right)
    +\langle\overline{{\varepsilon_h^\sigma}}\mathbf n,{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}}\rangle_{\partial{{\mathcal T_h}}}
    +\langle{\boldsymbol}\tau \, \overline{{{\boldsymbol}\varepsilon_h^u}}, {{\mathbf P_M}}{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}}\rangle_{\partial{{\mathcal T_h}}}.\end{aligned}$$

Substituting $({{\Pi_V}}{\boldsymbol}\sigma,{{\Pi_W}}\mathbf u,{{\mathbf P_M}}\mathbf u)$, where $({\boldsymbol}\sigma,\mathbf u)$ is the exact solution of in the left-hand side of discrete equations, and subtracting the actual discrete equations (with better written in the form ), it is simple to prove that the errors $({{e_h^\sigma}},{{\mathbf e_h^u}},{{\widehat{\mathbf e}_h^u}}) \in {{\mathbb V_h}}\times {{\mathbf W_h}}\times {{\mathbf M_h}}$ satisfy

\[eq:5.7\] $$\begin{aligned}
{6}
\label{eq:5.7a}
\imath\kappa (\mathcal A {{e_h^\sigma}},{\boldsymbol}\xi)_{{\mathcal T_h}}-({{\mathbf e_h^u}},\nabla\cdot{\boldsymbol}\xi)_{{\mathcal T_h}}+\langle {{\widehat{\mathbf e}_h^u}},{\boldsymbol}\xi\mathbf n\rangle_{\partial{{\mathcal T_h}}}
    &=\imath\kappa (\mathcal A {{\varepsilon_h^\sigma}},{\boldsymbol}\xi)_{{\mathcal T_h}},\\
\nonumber
(\nabla\cdot{{e_h^\sigma}},\mathbf w)_{{\mathcal T_h}}+\imath\kappa (\rho{{\mathbf e_h^u}},\mathbf w)_{{\mathcal T_h}}\\
\label{eq:5.7b}
    +\langle{\boldsymbol}\tau({{\mathbf P_M}}{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}}),{{\mathbf P_M}}\mathbf w\rangle_{\partial{{\mathcal T_h}}} 
    &=\imath\kappa(\rho{{{\boldsymbol}\varepsilon_h^u}},\mathbf w)_{{\mathcal T_h}}\\
\nonumber
    &\phantom{=}+\langle {{\varepsilon_h^\sigma}}\mathbf n,\mathbf w\rangle_{\partial{{\mathcal T_h}}}
        +\langle{\boldsymbol}\tau{{{\boldsymbol}\varepsilon_h^u}},{{\mathbf P_M}}\mathbf w\rangle_{\partial{{\mathcal T_h}}},\\
\label{eq:5.7c}
\langle {{e_h^\sigma}}\mathbf n+{\boldsymbol}\tau({{\mathbf P_M}}{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}}),{\boldsymbol}\mu\rangle_{\partial{{\mathcal T_h}}\setminus\Gamma_D}
    &=\langle{{\varepsilon_h^\sigma}}\mathbf n,{\boldsymbol}\mu\rangle_{\partial{{\mathcal T_h}}\setminus\Gamma_D}
        +\langle{\boldsymbol}\tau {{\mathbf P_M}}{{{\boldsymbol}\varepsilon_h^u}},{\boldsymbol}\mu\rangle_{\partial{{\mathcal T_h}}\setminus\Gamma_D}\\
\label{eq:5.7d}
\langle{{\widehat{\mathbf e}_h^u}},{\boldsymbol}\mu\rangle_{\Gamma_D} &=0,\end{aligned}$$

for all $({\boldsymbol}\xi,\mathbf w,{\boldsymbol}\mu)\in{{\mathbb V_h}}\times {{\mathbf W_h}}\times {{\mathbf M_h}}$. Clearly is equivalent to ${{\widehat{\mathbf e}_h^u}}=\mathbf 0$ on $\Gamma_D$. We now (i) test with ${\boldsymbol}\xi=\overline{{e_h^\sigma}}$, (ii) conjugate and test the result with $\mathbf w=\overline{{\mathbf e_h^u}}$, (iii) conjugate and test the result with ${\boldsymbol}\mu=-\overline{{\widehat{\mathbf e}_h^u}}$. The results of these three steps are added and reorganized ( is used for this) to prove the proposition.

We next bound the last two terms in the right-hand side of the identity in Proposition \[prop:5.3\]. From this moment on, we will frequently use, without additional warning, approximation properties of the $L^2$ projections onto the space of piecewise polynomial functions.

\[prop:5.4\] If $k\ge 1$, $$\begin{aligned}
{4}
 \left|\langle\overline{{\varepsilon_h^\sigma}}\mathbf n,{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}}\rangle_{\partial{{\mathcal T_h}}}
    +\langle{\boldsymbol}\tau \, \overline{{{\boldsymbol}\varepsilon_h^u}}, {{\mathbf P_M}}{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}}\rangle_{\partial{{\mathcal T_h}}}\right|
& \le C_1 
    \left( h^t |{\boldsymbol}\sigma |_{t,\Omega} + h^{s-1} |\mathbf u|_{s,\Omega}\right) \| {{\mathbf P_M}}{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}}\|_\tau \\
&  \phantom{\le}
    +C_2 \kappa   h^t |{\boldsymbol}\sigma|_{t,\Omega}  
        \left(h^t |{\boldsymbol}\sigma|_{t,\Omega} +\|{{e_h^\sigma}}\|_{\mathcal A}\right),\end{aligned}$$ if ${\boldsymbol}\sigma\in H^t(\Omega;{{\mathbb C}}^{3\times 3})$ for $1\le t \le k+1$ and $\mathbf u\in H^s(\Omega;{{\mathbb C}}^3)$ for $1\le s\le k+2$.

Following [@QiSh:2014 Lemma 4.3], it is possible to show that

\[eq:5.8\] $$\begin{aligned}
    \label{eq:5.8a}
\left|\langle\overline{{\varepsilon_h^\sigma}}\mathbf n,{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}}\rangle_{\partial{{\mathcal T_h}}}\right|
    & \le & C h^t |{\boldsymbol}\sigma |_{t,\Omega} 
        \left( \| {{\mathbf P_M}}{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}}\|_\tau + \| {\boldsymbol}\varepsilon({{\mathbf e_h^u}})\|_{{\mathcal T_h}}\right)  \\
\left| \langle{\boldsymbol}\tau \, \overline{{{\boldsymbol}\varepsilon_h^u}}, {{\mathbf P_M}}{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}}\rangle_{\partial{{\mathcal T_h}}}\right|
    &\le & C h^{s-1}|\mathbf u|_{s,\Omega}\| {{\mathbf P_M}}{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}}\|_\tau.
    \label{eq:5.8b}\end{aligned}$$

We recall that the argument leading to the proof of needs the traces of rigid motions to be in ${{\mathbf M_h}}$, which is where the additional hypothesis $k\ge 1$ is used.

We next test the first error equation with ${\boldsymbol}\xi={\boldsymbol}\varepsilon(\overline{{\mathbf e_h^u}})$ restricted to $K$ to obtain $$\label{eq:5.9}
\| {\boldsymbol}\varepsilon({{\mathbf e_h^u}})\|_K^2 = (\nabla {{\mathbf e_h^u}},{\boldsymbol}\varepsilon(\overline{{\mathbf e_h^u}}))_K
    =\imath\kappa ( \mathcal A({{\varepsilon_h^\sigma}}-{{e_h^\sigma}}),{\boldsymbol}\varepsilon(\overline{{\mathbf e_h^u}}))_K
        +\langle {{\mathbf P_M}}{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}},{\boldsymbol}\varepsilon(\overline{{\mathbf e_h^u}})\mathbf n\rangle_{\partial K}.$$ The scaling hypothesis on ${\boldsymbol}\tau$ given in and a scaling argument using the fact that ${\boldsymbol}\varepsilon({{\mathbf e_h^u}})$ is a polynomial on $K$ show then $$\begin{aligned}
\left|\langle {{\mathbf P_M}}{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}},{\boldsymbol}\varepsilon(\overline{{\mathbf e_h^u}})\mathbf n\rangle_{\partial K}\right|
    &\le & \| {{\mathbf P_M}}{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}}\|_{\partial K} \| {\boldsymbol}\varepsilon({{\mathbf e_h^u}})\|_{\partial K} \\
    &\le & C \| {\boldsymbol}\tau_K^{1/2} ({{\mathbf P_M}}{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}})\|_{\partial K} 
        h_K^{1/2} \| {\boldsymbol}\varepsilon({{\mathbf e_h^u}})\|_{\partial K} \\
    &\le & C' \| {\boldsymbol}\tau_K^{1/2} ({{\mathbf P_M}}{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}})\|_{\partial K}  \| {\boldsymbol}\varepsilon({{\mathbf e_h^u}})\|_K.\end{aligned}$$ Substituting these bounds in the right-hand side of , and adding over all elements, using (the spectral bounds on the inverse compliance tensor) it follows that $$\label{eq:5.10}
\| {\boldsymbol}\varepsilon({{\mathbf e_h^u}})\|_{{\mathcal T_h}}\le C \kappa \left( \| {{\varepsilon_h^\sigma}}\|_{\mathcal A}+\| {{e_h^\sigma}}\|_{\mathcal A}\right)
            + C \| {{\mathbf P_M}}{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}}\|_\tau.$$ Plugging in , the proposition is proved.

Dual problem and bootstrapping process
======================================

We consider the adjoint system to :

\[eq:6.1\] $$\begin{aligned}
{6}
-\imath\kappa\mathcal A{\boldsymbol}\psi -{\boldsymbol}\varepsilon({\boldsymbol}\phi) &=\mathbf 0
    &\qquad &\mbox{in $\Omega$},\\
-\nabla\cdot{\boldsymbol}\psi -\imath\kappa\,\rho\,{\boldsymbol}\phi &={{\mathbf e_h^u}}&\qquad &\mbox{in $\Omega$},\\
{\boldsymbol}\phi &=\mathbf 0
    &\qquad &\mbox{on $\Gamma_D$},\\
{\boldsymbol}\psi\mathbf n &=\mathbf 0
    &\qquad &\mbox{on $\Gamma_N$}.\end{aligned}$$

This problem is uniquely solvable if , or equivalently, , is. If we assume regularity for the coercive problem , expressed in the bound , and assume the local smoothness of the coefficients given at the beginning of Section \[sec:3\], then it is easy to see that ${\boldsymbol}\phi\in H^2(\Omega;\mathbb C^3)$ and ${\boldsymbol}\psi\in H^1(\Omega;\mathbb C^{3\times3})$. Morever, we can bound $$\label{eq:6.2}
\| {\boldsymbol}\phi\|_{1,\Omega} + \| \rho{\boldsymbol}\phi\|_{1,{{\mathcal T_h}}}+\| \mathcal A {\boldsymbol}\psi\|_{1,\Omega}+ \|{\boldsymbol}\psi\|_{1,{{\mathcal T_h}}}
+\kappa^{-1} \|{\boldsymbol}\phi\|_{2,\Omega} \le D_\kappa\|{{\mathbf e_h^u}}\|_\Omega,$$ where $D_\kappa:= C\,(1+\kappa C_\kappa+C_\kappa)$, $C_\kappa$ being the constant in and $C$ being allowed to depend on the regularity constant as well as on the physical coefficients. Here the norm $\|\cdot\|_{1,{{\mathcal T_h}}}$ is the natural norm of the broken Sobolev space $\prod_{K\in {{\mathcal T_h}}} H^1(K)$. We note that, while the regularity requirement can be somewhat relaxed, the analysis in this paper (see also [@QiSh:2014] and [@GrMo:2011]) needs a certain amount of regularity for the solution of the dual problem (which can be translated to regularity of the solution of ) due to the need of having square integrable normal traces of ${\boldsymbol}\psi$ on the faces of the elements.

\[prop:6.1\] $$\begin{aligned}
\| {{\mathbf e_h^u}}\|_\Omega^2
    &=& \imath\kappa \left( (\mathcal A{{e_h^\sigma}},\overline{{\boldsymbol}\psi}-{{\Pi_V}}\overline{{\boldsymbol}\psi})_{{\mathcal T_h}}+(\mathcal A{{\varepsilon_h^\sigma}},{{\Pi_V}}\overline{{\boldsymbol}\psi})_{{\mathcal T_h}}\right)\\    
    & & +\imath\kappa \left((\rho\,{{\mathbf e_h^u}},\overline{{\boldsymbol}\phi}-{{\Pi_W}}\overline{{\boldsymbol}\phi})_{{\mathcal T_h}}+(\rho\,{{{\boldsymbol}\varepsilon_h^u}},{{\Pi_W}}\overline{{\boldsymbol}\phi})_{{\mathcal T_h}}\right)\\
    & & +\langle {{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}},(\overline{{\boldsymbol}\psi}-{{\Pi_V}}\overline{{\boldsymbol}\psi})\mathbf n\rangle_{\partial{{\mathcal T_h}}}\\
    & & +\langle{\boldsymbol}\tau({{\mathbf P_M}}{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}})-{\boldsymbol}\tau{{\mathbf P_M}}{{{\boldsymbol}\varepsilon_h^u}},\overline{{\boldsymbol}\phi}-{{\Pi_W}}\overline{{\boldsymbol}\phi}\rangle_{\partial{{\mathcal T_h}}}
           + \langle {{\varepsilon_h^\sigma}}\mathbf n,{{\Pi_W}}\overline{{\boldsymbol}\phi}-{{\mathbf P_M}}\overline{{\boldsymbol}\phi}\rangle_{\partial{{\mathcal T_h}}}.\end{aligned}$$

The proof is similar to duality arguments in [@QiSh:2014] (and related references). We give here a very systematic approach to help understand the logic of the argument. We first conjugate equations and then test them with the discrete errors $({{e_h^\sigma}},{{\mathbf e_h^u}},{{\widehat{\mathbf e}_h^u}})$:

\[eq:6.2D\] $$\begin{aligned}
{6}
\label{eq:6.2a}
\imath\kappa(\mathcal A{{e_h^\sigma}},\overline{{\boldsymbol}\psi})_{{\mathcal T_h}}+(\nabla\cdot{{e_h^\sigma}},{{\Pi_W}}\overline{{\boldsymbol}\phi})_{{\mathcal T_h}}-\langle{{e_h^\sigma}}\mathbf n,\overline{{\boldsymbol}\phi}\rangle_{\partial{{\mathcal T_h}}} &=0,\\
\label{eq:6.2b}
-({{\mathbf e_h^u}},\nabla\cdot{{\Pi_V}}\overline{{\boldsymbol}\psi})_{{\mathcal T_h}}+\langle {{\mathbf e_h^u}},({{\Pi_V}}\overline{{\boldsymbol}\psi}-\overline{{\boldsymbol}\psi})\mathbf n\rangle_{\partial{{\mathcal T_h}}}
    +\imath\kappa (\rho\,{{\mathbf e_h^u}},\overline{{\boldsymbol}\phi})_{{\mathcal T_h}}& = \| {{\mathbf e_h^u}}\|_\Omega^2,\\
\label{eq:6.2c} 
\langle {{\widehat{\mathbf e}_h^u}},\overline{{\boldsymbol}\psi}\mathbf n\rangle_{\partial{{\mathcal T_h}}} &=0.\end{aligned}$$

Note that to reach and we need to integrate by parts and introduce projections wherever possible. Also, reflects the fact that ${\boldsymbol}\psi$ does not jump across interelement faces as well as the equality ${{\widehat{\mathbf e}_h^u}}=\mathbf 0$ on $\Gamma_D$. The second ingredient for the proof is the set of error equations tested with $({{\Pi_V}}\overline{{\boldsymbol}\psi},{{\Pi_W}}\overline{{\boldsymbol}\phi},{{\mathbf P_M}}\overline{{\boldsymbol}\phi})$, to yield

\[eq:6.3\] $$\begin{aligned}
{6}
\imath\kappa (\mathcal A {{e_h^\sigma}},\overline{{\boldsymbol}\psi})_{{\mathcal T_h}}-({{\mathbf e_h^u}},\nabla\cdot {{\Pi_V}}\overline{{\boldsymbol}\psi})_{{\mathcal T_h}}+\langle {{\widehat{\mathbf e}_h^u}},\overline{{\boldsymbol}\psi}\mathbf n\rangle_{\partial{{\mathcal T_h}}}
    &=\ell_1({\boldsymbol}\psi),\\
(\nabla\cdot{{e_h^\sigma}},{{\Pi_W}}\overline{{\boldsymbol}\phi})_{{\mathcal T_h}}+\imath\kappa (\rho{{\mathbf e_h^u}},\overline{{\boldsymbol}\phi})_{{\mathcal T_h}}&=\ell_2({\boldsymbol}\phi),\\
-\langle {{e_h^\sigma}}\mathbf n,\overline{{\boldsymbol}\phi}\rangle_{\partial{{\mathcal T_h}}}
    &=\ell_3({\boldsymbol}\phi),\end{aligned}$$

where $$\begin{aligned}
\ell_1({\boldsymbol}\psi)&:=&\imath\kappa \left(
    (\mathcal A{{e_h^\sigma}},\overline{{\boldsymbol}\psi}-{{\Pi_V}}\overline{{\boldsymbol}\psi})_{{\mathcal T_h}}+(\mathcal A {{\varepsilon_h^\sigma}},{{\Pi_V}}\overline{{\boldsymbol}\psi})_{{\mathcal T_h}}\right)
    +\langle{{\widehat{\mathbf e}_h^u}},(\overline{{\boldsymbol}\psi}-{{\Pi_V}}\overline{{\boldsymbol}\psi})\mathbf n\rangle_{\partial{{\mathcal T_h}}},\\
\ell_2({\boldsymbol}\phi) &:=& \imath\kappa\left(
    (\rho{{\mathbf e_h^u}},\overline{{\boldsymbol}\phi}-{{\Pi_W}}\overline{{\boldsymbol}\phi})_{{\mathcal T_h}}+(\rho{{{\boldsymbol}\varepsilon_h^u}},{{\Pi_W}}\overline{{\boldsymbol}\phi})_{{\mathcal T_h}}\right)\\
     & & +\langle {{\varepsilon_h^\sigma}}\mathbf n,{{\Pi_W}}\overline{{\boldsymbol}\phi}\rangle_{\partial{{\mathcal T_h}}}
        +\langle{\boldsymbol}\tau{{\mathbf P_M}}{{{\boldsymbol}\varepsilon_h^u}},{{\Pi_W}}\overline{{\boldsymbol}\phi}\rangle_{\partial{{\mathcal T_h}}}
         -\langle{\boldsymbol}\tau({{\mathbf P_M}}{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}}),{{\Pi_W}}\overline{{\boldsymbol}\phi}\rangle_{\partial{{\mathcal T_h}}} \\
\ell_3({\boldsymbol}\phi) &:=& -\langle{{\varepsilon_h^\sigma}}\mathbf n,{{\mathbf P_M}}\overline{{\boldsymbol}\phi}\rangle_{\partial{{\mathcal T_h}}}
        -\langle{\boldsymbol}\tau {{\mathbf P_M}}{{{\boldsymbol}\varepsilon_h^u}},\overline{{\boldsymbol}\phi}\rangle_{\partial{{\mathcal T_h}}}
        +\langle {\boldsymbol}\tau({{\mathbf P_M}}{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}}),\overline{{\boldsymbol}\phi}\rangle_{\partial{{\mathcal T_h}}}.\end{aligned}$$ Note that in we have kept in the left-hand side of the error equations only those terms that appear in the left-hand side of . We have also eliminated some redundant projections and applied that ${\boldsymbol}\phi=\mathbf 0$ on $\Gamma_D$. The proof of the result is now straightforward: add equations and substitute equations in the result.

The next step in the proof of the error estimates is a bound for $\|{{\mathbf e_h^u}}\|_\Omega$ obtained by carefully working on the right-hand side of the duality identity in Proposition \[prop:6.1\]. To alleviate the proof from an excess of constants, we will use the convention that $a\lesssim b$, whenever there exists a positive constant $C$ independent of $h$ and $\kappa$ such that $a\le C \, b$.

\[prop:6.2\] If $h \kappa D_\kappa$ is small enough if ${\boldsymbol}\sigma\in H^t(\Omega;{{\mathbb C}}^{3\times 3})$ for $1\le t \le k+1$ and if $\mathbf u\in H^s(\Omega;{{\mathbb C}}^3)$ for $1\le s\le k+2$, then $$\| {{\mathbf e_h^u}}\|_\Omega \lesssim  h (\kappa+1) D_\kappa 
\left( \| {{e_h^\sigma}}\|_{\mathcal A}+\| {{\mathbf P_M}}{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}}\|_\tau
    + h^t |{\boldsymbol}\sigma|_{t,\Omega}
    +h^{s-1} |\mathbf u|_{s,\Omega}\right).$$

As already mentioned, we first estimate the right-hand side of the equality in Proposition \[prop:6.1\]. Let ${[ f ]_h}$ be the best $L^2(\Omega)$ projection of $f$ on the space of piecewise constant functions, i.e., ${[ f ]_h}=\frac1{|K|}\int_K f$ in $K$ for every $K$.

Notice that $$(\mathcal A{{e_h^\sigma}},\overline{{\boldsymbol}\psi}-{{\Pi_V}}\overline{{\boldsymbol}\psi})_{{\mathcal T_h}}+(\mathcal A{{\varepsilon_h^\sigma}},{{\Pi_V}}\overline{{\boldsymbol}\psi})_{{\mathcal T_h}}= (\mathcal A({\boldsymbol}\sigma-{\boldsymbol}\sigma_h),\overline{{\boldsymbol}\psi}-{{\Pi_V}}\overline{{\boldsymbol}\psi})_{{\mathcal T_h}}+ ({{\varepsilon_h^\sigma}},\mathcal A\overline{{\boldsymbol}\psi}-{[ \mathcal A\overline{{\boldsymbol}\psi} ]_h})_{{\mathcal T_h}},$$ and then we can bound $$\label{eq:6.5}
|(\mathcal A{{e_h^\sigma}},\overline{{\boldsymbol}\psi}-{{\Pi_V}}\overline{{\boldsymbol}\psi})_{{\mathcal T_h}}+(\mathcal A{{\varepsilon_h^\sigma}},{{\Pi_V}}\overline{{\boldsymbol}\psi})_{{\mathcal T_h}}|
    \lesssim h \|{\boldsymbol}\sigma-{\boldsymbol}\sigma_h\|_{\mathcal A}  |{\boldsymbol}\psi|_{1,{{\mathcal T_h}}}
        + h \|{{\varepsilon_h^\sigma}}\|_\Omega  |\mathcal A {\boldsymbol}\psi|_{1,\Omega}.$$ Similarly, the equality $$(\rho\,{{\mathbf e_h^u}},\overline{{\boldsymbol}\phi}-{{\Pi_W}}\overline{{\boldsymbol}\phi})_{{\mathcal T_h}}+(\rho\,{{{\boldsymbol}\varepsilon_h^u}},{{\Pi_W}}\overline{{\boldsymbol}\phi})_{{\mathcal T_h}}=(\rho(\mathbf u-\mathbf u_h),\overline{{\boldsymbol}\phi}-{{\Pi_W}}\overline{{\boldsymbol}\phi})_{{\mathcal T_h}}+({{{\boldsymbol}\varepsilon_h^u}},\rho\overline{{\boldsymbol}\phi}-{[ \rho\overline{{\boldsymbol}\phi} ]_h})_{{\mathcal T_h}}$$ can be used to estimate $$\label{eq:6.6}
|(\rho\,{{\mathbf e_h^u}},\overline{{\boldsymbol}\phi}-{{\Pi_W}}\overline{{\boldsymbol}\phi})_{{\mathcal T_h}}+(\rho\,{{{\boldsymbol}\varepsilon_h^u}},{{\Pi_W}}\overline{{\boldsymbol}\phi})_{{\mathcal T_h}}|
\lesssim h \|\mathbf u-\mathbf u_h\|_\Omega | {\boldsymbol}\phi|_{1,\Omega}
    + h \|{{{\boldsymbol}\varepsilon_h^u}}\|_\Omega |\rho\,{\boldsymbol}\phi|_{1,{{\mathcal T_h}}}$$ Using with ${\boldsymbol}\psi$ in place of ${\boldsymbol}\sigma$ and $t=1$ and , we can estimate $$\label{eq:6.7}
|\langle {{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}},(\overline{{\boldsymbol}\psi}-{{\Pi_V}}\overline{{\boldsymbol}\psi})\mathbf n\rangle_{\partial{{\mathcal T_h}}}|
\lesssim h \left( \| {{\varepsilon_h^\sigma}}\|_{\mathcal A}+\| {{e_h^\sigma}}\|_{\mathcal A}+ \| {{\mathbf P_M}}{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}}\|_\tau\right)
|{\boldsymbol}\psi|_{1,{{\mathcal T_h}}}.$$ Using with ${\boldsymbol}\phi$ in place of $\mathbf u$ and $s=2$, we bound $$\label{eq:6.8}
|\langle{\boldsymbol}\tau({{\mathbf P_M}}{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}}),\overline{{\boldsymbol}\phi}-{{\Pi_W}}\overline{{\boldsymbol}\phi}\rangle_{\partial{{\mathcal T_h}}}|
\lesssim h \|{{\mathbf P_M}}{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}}\|_\tau | {\boldsymbol}\phi|_{2,\Omega}.$$ With a scaling argument and the bound (stating the size of the stabilization paraneter), we can bound on every $K$ $$\begin{aligned}
|\langle {\boldsymbol}\tau {{\mathbf P_M}}(\mathbf u-{{\Pi_W}}\mathbf u),\overline{{\boldsymbol}\phi}-{{\Pi_W}}\overline{{\boldsymbol}\phi}\rangle_{\partial K}|
    &\lesssim & h_K^{-1} \| \mathbf u-{{\Pi_W}}\mathbf u\|_{\partial K}
                        \| {\boldsymbol}\phi-{{\Pi_W}}{\boldsymbol}\phi\|_{\partial K}\\
    & \lesssim & h_K^s |\mathbf u|_{s,K}|{\boldsymbol}\phi|_{2,K}\end{aligned}$$ and therefore $$\label{eq:6.9}
|\langle {\boldsymbol}\tau{{\mathbf P_M}}{{{\boldsymbol}\varepsilon_h^u}},\overline{{\boldsymbol}\phi}-{{\Pi_W}}\overline{{\boldsymbol}\phi}\rangle_{\partial{{\mathcal T_h}}}|
\lesssim h^s |\mathbf u|_{s,\Omega} |{\boldsymbol}\phi|_{2,\Omega}.$$ Similarly $$\label{eq:6.10}
|\langle {{\varepsilon_h^\sigma}}\mathbf n,{{\mathbf P_M}}\overline{{\boldsymbol}\phi}-{{\Pi_W}}\overline{{\boldsymbol}\phi}\rangle_{\partial{{\mathcal T_h}}}|
\lesssim h^{t+1} |{\boldsymbol}\sigma|_{t,\Omega} |{\boldsymbol}\phi|_{2,\Omega}.$$ Collecting the estimates - to bound the right-hand side of the identity in Proposition \[prop:6.1\], and using the regularity bound , we can bound $$\begin{aligned}
\| {{\mathbf e_h^u}}\|_\Omega^2
    \lesssim  \, & h \kappa D_\kappa \| {{\mathbf e_h^u}}\|_\Omega 
    \left( \| {{\mathbf e_h^u}}\|_\Omega + \|{{e_h^\sigma}}\|_{\mathcal A} +\| {{\mathbf P_M}}{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}}\|_\tau \right) \\
    & + h (1+\kappa ) D_\kappa\left(h^t |{\boldsymbol}\sigma|_{t,\Omega}  + h^{s-1} |\mathbf u|_{s,\Omega} \right)\end{aligned}$$ The proposition is now a simple consequence of the latter inequality.

The proof of Theorem \[the:3.1\] follows from the energy identity (Proposition \[prop:5.3\]) and the estimates of Propositions \[prop:5.4\] and \[prop:6.2\] by a careful bootstrapping process.

To simplify the algebra involved in this final step, let us give symbols for the quantities we want to bound $$\Sigma:=\| {{e_h^\sigma}}\|_{\mathcal A}, \qquad 
    \mathrm T:=\kappa^{-1/2}\| {{\mathbf P_M}}{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}}\|_\tau, \qquad
    \mathrm U:=\| {{\mathbf e_h^u}}\|_\Omega,$$ and for the approximation terms $$\Sigma_h:= h^t |{\boldsymbol}\sigma|_{t,\Omega},\qquad
    \mathrm U_h:=h^{s-1} |\mathbf u|_{s,\Omega}.$$ With this shorthand, Propositions \[prop:5.3\] and \[prop:5.4\] yield $$\label{eq:6.11}
|\Sigma^2-\imath \mathrm T^2|
    \lesssim \Sigma\,\Sigma_h + \mathrm U\,\mathrm U_h 
    + \kappa^{-1/2} (\Sigma_h+\mathrm U_h) \mathrm T +\Sigma_h^2+\mathrm U^2.$$ If $\alpha:=h (1+\kappa) (1+\kappa C_\kappa)=D_\kappa h (1+\kappa)$, Proposition \[prop:6.2\] can then be rephrased as $$\label{eq:6.12}
\mathrm U \lesssim \alpha (\Sigma+\Sigma_h+\mathrm U_h +\kappa^{1/2}\mathrm T).$$ Substituting in the right-hand side of , and reordering terms, we have $$\begin{aligned}
\label{eq:6.13}
\nonumber
\Sigma^2+\mathrm T^2 
    & \lesssim &
    \alpha^2\Sigma^2+\Sigma(\Sigma_h+\alpha \mathrm U_h)\\
    && +\alpha^2 \kappa \mathrm T^2
    +\mathrm T(\kappa^{-1/2}\mathrm U_h+\alpha \kappa^{1/2}\mathrm U_h
                +\kappa^{-1/2} \Sigma_h)\\
\nonumber
    & &
    +(1+ \alpha^2) \Sigma_h^2
     + (\alpha+\alpha^2)\mathrm U_h^2.\end{aligned}$$ Let now $C$ be the constant that is hidden in the symbol $\lesssim$, and let us assume that $$\label{eq:hypo}
C \alpha^2\le \tfrac14 \qquad \mbox{and} \qquad C \alpha^2\kappa\le \tfrac14.$$ We now use Young’s inequality $ab\le \frac14 a^2+ b^2$ in to get $$\begin{aligned}
\Sigma^2+\mathrm T^2  
        & \le & \tfrac12\Sigma^2+\tfrac12 \mathrm T^2 \\
        & & + C^2(\Sigma_h+\alpha \mathrm U_h)^2
            +C^2(\kappa^{-1/2}\mathrm U_h+\alpha \kappa^{1/2}\mathrm U_h+\kappa^{-1/2} \Sigma_h)^2\\
        & & +C(1+ \alpha^2) \Sigma_h^2
            + C(\alpha+\alpha^2)\mathrm U_h^2.\end{aligned}$$ We can now simplify this expression using to obtain $\Sigma^2+\mathrm T^2 \lesssim (1+\kappa^{-1})(\Sigma_h^2+\mathrm U_h^2)$, or equivalently $$\label{eq:6.14}
\Sigma+\mathrm T 
    \lesssim 
    (1+\kappa^{-1/2})(\Sigma_h+\mathrm U_h).$$ Using in , we can finally prove that $$\label{eq:6.15}
\mathrm U \lesssim \alpha(1+\kappa^{1/2}+\kappa^{-1/2})(\Sigma_h+\mathrm U_h).$$ This finishes the proof.

Variants and insights {#sec:6}
=====================

#### Matrix form.

We first give a matrix representation of the method of Section \[sec:3\]. Equations and suggest the following orthogonal decomposition ${{\mathbf M_h}}={{\mathbf M_h^{nD}}}\oplus {{\mathbf M_h^{D}}}$, where ${{\mathbf M_h^{nD}}}=\{{\boldsymbol}\mu\,:\,{\boldsymbol}\mu|_{\Gamma_D}=0\}\equiv \{ {\boldsymbol}\mu|_{\partial{{\mathcal T_h}}\setminus\Gamma_D}\,:\,{\boldsymbol}\mu\in {{\mathbf M_h}}\}$. We now take real-valued bases for the spaces ${{\mathbb V_h}}$, ${{\mathbf W_h}}$, ${{\mathbf M_h^{nD}}}$, and ${{\mathbf M_h^{D}}}$ and identify the unknowns ${\boldsymbol}\sigma_h\in{{\mathbb V_h}}$, ${\mathbf}u_h\in{{\mathbf W_h}}$, $\widehat{{\mathbf}u}_h|_{\partial{{\mathcal T_h}}\setminus\Gamma_D}\in {{\mathbf M_h^{nD}}}$, and $\widehat{{\mathbf}u}_h|_{\Gamma_D}\in {{\mathbf M_h^{D}}}$, with respective complex column vectors $\underline\sigma$, $\underline u$, $\underline{\widehat u}^{{nD}}$, and $\underline{\widehat u}^{{D}}$. We then consider [*real matrices*]{} associated to the following [*bilinear forms*]{} that are understood as functionals acting on the unknowns: $$\begin{aligned}
{6}
(\mathcal A {\boldsymbol}\sigma_h,{\boldsymbol}\xi)_{{\mathcal T_h}}& \qquad  & {\boldsymbol}\xi\in {{\mathbb V_h}}& \qquad & & \mathrm A \underline\sigma, \\
(\nabla\cdot{\boldsymbol}\sigma_h,{\mathbf}w)_{{\mathcal T_h}}& & {\mathbf}w\in {{\mathbf W_h}}& & & \mathrm D\underline\sigma, \\
\langle{\boldsymbol}\sigma_h{\mathbf}n,{\boldsymbol}\mu\rangle_{\partial{{\mathcal T_h}}\setminus\Gamma_D} & & {\boldsymbol}\mu\in {{\mathbf M_h^{nD}}} 
        & & &\mathrm N \underline\sigma,\\
\langle{\boldsymbol}\sigma_h{\mathbf}n,{\boldsymbol}\mu\rangle_{\Gamma_D} & & {\boldsymbol}\mu\in {{\mathbf M_h^{D}}} 
        & & &\mathrm N_{{D}} \underline\sigma,\\
(\rho{\mathbf}u_h,{\mathbf}w)_{{\mathcal T_h}}& & {\mathbf}w\in {{\mathbf W_h}}& & & \mathrm M\underline u,\\
\langle {\boldsymbol}\tau{{\mathbf P_M}}{\mathbf}u_h,{{\mathbf P_M}}{\mathbf}w\rangle_{{{\mathcal T_h}}} & & {\mathbf}w\in {{\mathbf W_h}}& & & \mathrm T_{11} \underline u,\\
\langle{\boldsymbol}\tau\widehat{{\mathbf}u}_h,{\mathbf}w_h\rangle_{\partial{{\mathcal T_h}}\setminus\Gamma_D} & & {\mathbf}w\in {{\mathbf W_h}}& & & \mathrm T_{12}\underline u^{{nD}},\\
\langle{\boldsymbol}\tau\widehat{{\mathbf}u}_h,{\mathbf}w_h\rangle_{\Gamma_D} & & {\mathbf}w\in {{\mathbf W_h}}& & & \mathrm T_{D}\underline u^{{D}},\\
\langle{\boldsymbol}\tau\widehat{{\mathbf}u}_h,{\boldsymbol}\mu\rangle_{\partial{{\mathcal T_h}}\setminus\Gamma_D} & & {\boldsymbol}\mu\in {{\mathbf M_h^{nD}}} 
        & & & \mathrm T_{22}\underline u^{{nD}},\\
\langle\widehat{{\mathbf}u}_h,{\boldsymbol}\mu\rangle_{\Gamma_D} & &{\boldsymbol}\mu\in {{\mathbf M_h^{D}}}
        & & & \mathrm M_{{D}} \underline{\widehat u}^{{D}}.\end{aligned}$$Note that the matrices $\mathrm A$, $\mathrm M$, and $\mathrm T_{22}$ are symmetric and positive definite, while $\mathrm T_{11}$ is symmetric and positive semidefinite. The method given by equations is then equivalent to the linear system $$\label{eq:7.0}
\left[\begin{array}{cccc}
    \imath\kappa \mathrm A & -\mathrm D^\top & \mathrm N^\top & \mathrm N_{{D}}^\top \\
    \mathrm D & \imath\kappa\mathrm M+\mathrm T_{11} & -\mathrm T_{12} & -\mathrm T_{D} \\
    -\mathrm N & -\mathrm T_{12}^\top & \mathrm T_{22} & \mathrm O \\
    \mathrm O & \mathrm O &\mathrm O & \mathrm M_{{D}}
\end{array}\right]
\left[\begin{array}{l} 
    \underline\sigma \\ \underline u\\ \underline{\widehat u}^{{nD}} \\ \underline{\widehat u}^{{D}}
\end{array}\right]
=
\left[\begin{array}{c}
    \underline 0 \\ \underline f \\ -\underline g_N \\ \underline g_D
\end{array}\right],$$ where the definition of the vectors $\underline f$, $\underline g_N$ and $\underline g_D$ is self-evident. Equations are equivalent to the following system $$\left[\begin{array}{cccc}
    \mathrm A & \mathrm D^\top & -\mathrm N^\top & -\mathrm N_{{D}}^\top \\
    -\mathrm D & -\kappa^2\mathrm M-\imath\kappa\mathrm T_{11} & 
            \imath\kappa\mathrm T_{12} & \imath\kappa\mathrm T_{D} \\
    \mathrm N & -\imath\kappa\mathrm T_{12}^\top & \imath\kappa\mathrm T_{22} & \mathrm O \\
    \mathrm O & \mathrm O &\mathrm O & \mathrm M_{{D}}
\end{array}\right]
\left[\begin{array}{l} 
    \widetilde{\underline\sigma} \\ \underline u\\ \underline{\widehat u}^{{nD}} \\ \underline{\widehat u}^{{D}}
\end{array}\right]
=
\left[\begin{array}{c}
    \underline 0 \\ \widetilde{\underline f} \\ \widetilde{\underline g}_N \\ \underline g_D
\end{array}\right],$$ where $\widetilde{\underline\sigma}=-\imath\kappa\underline\sigma$, $\widetilde{\underline f}=-\imath\kappa \underline f$, and $\widetilde{\underline g}_N=-\imath\kappa\underline g_N$. This change of variables in the first unknown and in the right-hand side reverts the system to the original physical variables (the ones with a tilde in Section \[sec:2\]) so that the equations are second-order-in-frequency. It is clear how the stabilization terms are the only complex-valued ones in the system.

#### Hybridization.

The four matrices in the upper left $2\times 2$ block of the matrix of are elementwise block diagonal. The hybridization process consists of solving the system $$\begin{aligned}
{6}
\left( \left[\begin{array}{cc} \mathrm N & \mathrm T_{12}^\top \end{array}\right] \mathrm C^{-1} 
\left[\begin{array}{c} \mathrm N^\top \\ - \mathrm T_{12} \end{array}\right]+\mathrm T_{22} \right) \underline{\widehat u}^{{nD}}
 = & -\underline g_N + 
 \left[\begin{array}{cc} \mathrm N & \mathrm T_{12}^\top \end{array}\right] \mathrm C^{-1} 
\left[\begin{array}{c} \underline 0 \\ \underline f \end{array}\right] \\
 & - \left[\begin{array}{cc} \mathrm N & \mathrm T_{12}^\top \end{array}\right] \mathrm C^{-1} 
\left[\begin{array}{c} \mathrm N_D^\top \\ - \mathrm T_D \end{array}\right]\mathrm M_D^{-1} \underline g_D,\end{aligned}$$ where the invertibility of $$\mathrm C:= 
    \left[\begin{array}{cc} 
        \imath \kappa \mathrm A & -\mathrm D^\top \\
        \mathrm D & \imath \kappa\mathrm M+\mathrm T_{11}
    \end{array}\right]$$ was the object of Proposition \[prop:5.2\].

#### Variant \# 1: time reversal.

While we are not making any claims about the behavior of the method for high frequency problems, we have kept $\kappa$ visible everywhere. We next explore some variants of the method that can be obtained by changing to second-order-in-frequency form and exploring different choices of the stabilization parameter. The energy identity of Proposition \[prop:5.3\] is the trigger for the analysis of the method. It is there clear that the sign of the boundary term is not relevant and a method based on the numerical flux $$\widehat{{\boldsymbol}\sigma}_h\mathbf n:={\boldsymbol}\sigma_h\mathbf n-{\boldsymbol}\tau_K ({{\mathbf P_M}}\mathbf u_h-\widehat{\mathbf u}_h) 
\,:\, \partial K\to {{\mathbb C}}^3$$ has the same convergence properties as the method presented in Section \[sec:2\]. As we will see later in this section, this method corresponds to time reversal.

#### Variant \# 2: $\kappa$-scaled stabilization.

The factor $\kappa^{-1/2}$ in the error estimate of Theorem \[the:3.1\] suggests the following variant of the numerical method: we still use equations by change the definition of the numerical flux to be $$\label{eq:7.1}
\widehat{{\boldsymbol}\sigma}_h\mathbf n:={\boldsymbol}\sigma_h\mathbf n+\kappa {\boldsymbol}\tau_K ({{\mathbf P_M}}\mathbf u_h-\widehat{\mathbf u}_h) 
\,:\, \partial K\to {{\mathbb C}}^3.$$ (Note that, as shown in [@GoLaNiPe:2015] for the acoustic wave equation, making the stabilization parameter depend on $\kappa$ is a must when we want to deal with complex frequencies. This dependence has also some desirable properties.) The proof of Theorem \[the:3.1\] can be easily adapted to deal with the method whose stabilization term is given by . The error estimate is given in the following theorem.

There exist $C_1, C_2>0$, dependent only on the shape-regularity of ${{\mathcal T_h}}$, the density $\rho$ and the coefficients of the inverse compliance tensor $\mathcal A$ such that if $h (1+\kappa)^{3/2}(1+\kappa C_\kappa)$ is small enough, then the errors can be bounded by $$\| {{e_h^\sigma}}\|_{\mathcal A}+ \| {{\mathbf P_M}}{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}}\|_\tau
    \le C_1 \big((1+\kappa^{-1})  h^t |{\boldsymbol}\sigma|_{t,\Omega}+h^{s-1} |\mathbf u|_{s,\Omega}\big)$$ and $$\| {{\mathbf e_h^u}}\|_\Omega\le C_2 (1+\kappa C_\kappa) (1+\kappa) 
\big((\kappa+\kappa^{-1}) h^{t+1} |{\boldsymbol}\sigma|_{t,\Omega}+h^{s} (1+\kappa) |\mathbf u|_{s,\Omega}\big),$$ if $k\ge 1$, $\mathbf u\in H^s(\Omega;{{\mathbb C}}^3)$ with $1\le s\le k+2$, and ${\boldsymbol}\sigma\in H^t(\Omega;{{\mathbb C}}^{3\times 3})$ with $1\le t\le k+1$.

#### Second-order-in-frequency formulations.

Since all methods presented above are based in a first-order-in-frequence formulation, defining ${\boldsymbol}\sigma:=(\imath/\kappa)\widetilde{{\boldsymbol}\sigma}$, where $\widetilde{{\boldsymbol}\sigma}$ is the physical stress for the displacement field ${\mathbf}u$. Consider now the following family of HDG schemes based on a second-order-in-frequency formulation: the spaces are unchanged and $\alpha_\kappa$ is a fixed parameter that is allowed to depend on the frequency:

\[eq:7.2\] $$\begin{aligned}
{6}
 (\mathcal A\widetilde{{\boldsymbol}\sigma}_h,{\boldsymbol}\xi)_{{{\mathcal T_h}}}
    +(\mathbf u_h,\nabla\cdot{\boldsymbol}\xi)_{{\mathcal T_h}}-\langle\widehat{\mathbf u}_h,{\boldsymbol}\xi\mathbf n\rangle_{\partial{{\mathcal T_h}}}
     & = 0 & \qquad & \forall {\boldsymbol}\xi\in {{\mathbb V_h}},\\
(\widetilde{{\boldsymbol}\sigma}_h,\nabla\mathbf w)_{{{\mathcal T_h}}}
    - \langle \widehat{{\boldsymbol}\sigma}_h{\mathbf}n,\mathbf w\rangle_{{\boldsymbol}\tau}
    -\kappa^2 (\rho\,\mathbf u_h,\mathbf w)_{{\mathcal T_h}}&=-(\widetilde{\mathbf f},\mathbf w)_{{\mathcal T_h}}& & \forall \mathbf w\in {{\mathbf W_h}},\\
\langle\widehat{{\boldsymbol}\sigma}_h{\mathbf}n,{\boldsymbol}\mu\rangle_{\partial{{\mathcal T_h}}\setminus\Gamma_D}
    &=\langle\widetilde{\mathbf g}_N,{\boldsymbol}\mu\rangle_{\Gamma_N} & & \forall {\boldsymbol}\mu \in {{\mathbf M_h}},\\
\langle\widehat{\mathbf u}_h,{\boldsymbol}\mu\rangle_{\Gamma_D} &=\langle\mathbf g_D,{\boldsymbol}\mu\rangle_{\Gamma_D}
    & &\forall {\boldsymbol}\mu\in {{\mathbf M_h}},\end{aligned}$$ where $$\widehat{{\boldsymbol}\sigma}_h{\mathbf}n:=\widetilde{{\boldsymbol}\sigma}_h{\mathbf}n-\alpha_\kappa {\boldsymbol}\tau({{\mathbf P_M}}{\mathbf}u_h-\widehat{{\mathbf}u}_h).$$

Note that the equations are written in terms of the original data in . The choice $\alpha_\kappa=1$ is the direct application of the method in [@QiSh:2014] to the equation $\nabla\cdot\widetilde{{\boldsymbol}\sigma}+\kappa^2\rho{\mathbf}u=\widetilde{{\mathbf}f}$. This choice of the parameter $\alpha_\kappa$ yields a method that transitions smoothly (analytically) to the zero-frequency limit. Methods based on the first-order-in-frequency formulation can be rewritten in the form with the relation $\widetilde{{\boldsymbol}\sigma}_h=-\imath\kappa{\boldsymbol}\sigma_h$ and the parameter $\alpha_\kappa:=\imath\kappa$ (for the method of Section \[sec:3\]), $\alpha_\kappa:=-\imath\kappa$ (time reversed method) or $\alpha_\kappa=\imath\kappa^2$ (for the method with the flux defined in .

#### Variant \# 3: conservative method.

The error estimate for the method in with the choice $\alpha_\kappa=1$ is given in the next theorem. Not surprisingly the estimates hold as $\kappa\to 0$, since we end up with a smooth perturbation of the discretization of the steady-state equations. Note that when $\kappa\to 0$ and we are not dealing with the pure Neumann problem, the quantity $C_\kappa$ converges to a finite value. Later in this section we will see that this choice corresponds to a conservative method in the time domain.

\[the:7.2\] There exist $C_1, C_2>0$, dependent only on the shape-regularity of ${{\mathcal T_h}}$, the density $\rho$ and the coefficients of the inverse compliance tensor $\mathcal A$ such that if $h (1+\kappa) E_\kappa$ is small enough, then the errors can be bounded by $$\| {{\Pi_V}}\widetilde{{\boldsymbol}\sigma}-\widetilde{{\boldsymbol}\sigma}_h \|_{\mathcal A}+ \| {{\mathbf P_M}}{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}}\|_\tau
    \le C_1 \big(  h^t |\widetilde{{\boldsymbol}\sigma}|_{t,\Omega}+(1+\kappa) h^{s-1} |\mathbf u|_{s,\Omega}\big)$$ and $$\| {{\mathbf e_h^u}}\|_\Omega\le C_2 E_\kappa (1+\kappa)
\big( h^{t+1} |\widetilde{{\boldsymbol}\sigma}|_{t,\Omega}+h^{s} (1+\kappa) |\mathbf u|_{s,\Omega}\big),$$ if $k\ge 1$, $\mathbf u\in H^s(\Omega;{{\mathbb C}}^3)$ with $1\le s\le k+2$, and $\widetilde{{\boldsymbol}\sigma}\in H^t(\Omega;{{\mathbb C}}^{3\times 3})$ with $1\le t\le k+1$. Here $E_\kappa\le C (1+\kappa C_\kappa+ C_\kappa+\kappa C_\kappa^{1/2})$.

See Section \[sec:A\].

#### Methods in the time domain.

Let us now write some HDG-semidiscrete methods for the transient elastic wave equation. The data functions are $\widetilde{{\mathbf}f}:[0,\infty)\to L^2(\Omega;{{\mathbb R}}^3)$, ${\mathbf}g_D:[0,\infty)\to H^{1/2}(\Gamma_D;{{\mathbb R}}^3)$, and $\widetilde{{\mathbf}g}_N:[0,\infty)\to L^2(\Gamma_N;{{\mathbb R}}^3)$. We then look for ${\mathbf}u:[0,\infty)\to H^1(\Omega;{{\mathbb R}}^3)$ and ${\boldsymbol}\sigma:[0,\infty)\to H(\mathrm{div},\Omega;{{\mathbb R}}^{3\times 3}_{\mathrm{sym}})$ satisfying

\[eq:6.4\] $$\begin{aligned}
{6}
\mathcal A \widetilde{{\boldsymbol}\sigma}(t)-{\boldsymbol}\varepsilon (\mathbf u)(t)&={\boldsymbol}0&\qquad &\mbox{in $\Omega$}, \quad\forall t\ge 0,\\
\nabla\cdot\widetilde{{\boldsymbol}\sigma} -\rho\,\ddot{\mathbf u}(t)
     &= \widetilde{\mathbf f}(t) &\qquad &\mbox{in $\Omega$}, \quad\forall t\ge 0,\\
\mathbf u(t) &=\mathbf g_D(t) & & \mbox{on $\Gamma_D$}, \quad\forall t\ge 0,\\
\widetilde{{\boldsymbol}\sigma}(t)\,\mathbf n &=\widetilde{\mathbf g}_N(t) & & \mbox{on $\Gamma_N$}, \quad\forall t\ge 0,\end{aligned}$$

and initial conditions ${\mathbf}u(0)={\mathbf}u_0$, $\dot{{\mathbf}u}(0)={\mathbf}v_0$. The HDG semidiscretization consists of the search for $\widetilde{{\boldsymbol}\sigma}_h:[0,\infty)\to {{\mathbb V_h}}$, ${\mathbf}u_h:[0,\infty)\to {{\mathbf W_h}}$, and $\widehat{{\mathbf}u}_h:[0,\infty)\to {{\mathbf M_h}}$ satisfying

\[eq:7.5\] $$\begin{aligned}
{6}
\label{eq:7.5a}
 (\mathcal A \widetilde{{\boldsymbol}\sigma}_h(t),{\boldsymbol}\xi)_{{{\mathcal T_h}}}
    +(\mathbf u_h(t),\nabla\cdot{\boldsymbol}\xi)_{{\mathcal T_h}}-\langle\widehat{\mathbf u}_h(t),{\boldsymbol}\xi\mathbf n\rangle_{\partial{{\mathcal T_h}}}
     & = 0 & \\
\label{eq:7.5b}
(\widetilde{{\boldsymbol}\sigma}_h(t),\nabla\mathbf w)_{{{\mathcal T_h}}}
    - \langle \widehat{{\boldsymbol}\sigma}_h(t){\mathbf}n,\mathbf w\rangle_{\partial{{\mathcal T_h}}}
    + (\rho\,\ddot{\mathbf u}_h(t),\mathbf w)_{{\mathcal T_h}}&=-(\widetilde{\mathbf f}(t),\mathbf w)_{{\mathcal T_h}}& & \\
\label{eq:7.5c}
\langle\widehat{{\boldsymbol}\sigma}_h(t){\mathbf}n,{\boldsymbol}\mu\rangle_{\partial{{\mathcal T_h}}\setminus\Gamma_D}
    &=\langle\widetilde{\mathbf g}_N(t),{\boldsymbol}\mu\rangle_{\Gamma_N} & & \\
\label{eq:7.5d}
\langle\widehat{\mathbf u}_h(t),{\boldsymbol}\mu\rangle_{\Gamma_D} &=\langle\mathbf g_D(t),{\boldsymbol}\mu\rangle_{\Gamma_D}
    & &\end{aligned}$$

for all $({\boldsymbol}\xi,{\mathbf}w,{\boldsymbol}\mu)\in {{\mathbb V_h}}\times {{\mathbf W_h}}\times {{\mathbf M_h}}$ amd $t\ge 0$. The numerical flux $\widehat{{\boldsymbol}\sigma}_h$ can be defined in different ways which influence the choice of initial conditions. If we take the inverse Fourier transform of equations with $\alpha_\kappa=\pm \imath\kappa$, we obtain the following proposals for the numerical flux $$\label{eq:7.6}
\widehat{{\boldsymbol}\sigma}_h(t){\mathbf}n:=\widetilde{{\boldsymbol}\sigma}_h(t){\mathbf}n\pm {\boldsymbol}\tau({{\mathbf P_M}}\dot{{\mathbf}u}_h(t)-\dot{\widehat{{\mathbf}u}}_h(t)).$$ Note the positive sign corresponds to the method of Section \[sec:3\] while the negative sign is the one obtained by time reversal. (It is clear from this why the sign change in the parameter $\alpha_\kappa$ corresponded to time reversal.) For equations with $\alpha_\kappa=1$ we obtain the flux $$\label{eq:7.7}
\widehat{{\boldsymbol}\sigma}_h(t){\mathbf}n:=\widetilde{{\boldsymbol}\sigma}_h(t){\mathbf}n+{\boldsymbol}\tau({{\mathbf P_M}}{\mathbf}u_h(t)-\widehat{{\mathbf}u}_h(t)).$$ The following result shows that: the method with flux given by with positive sign accumulates energy over time, the method with flux with negative sign is dissipative, and the method with flux is conservative. The build-up or dissipation of energy happens at the interfaces, while the conservative method needs to add a potential energy term in the interfaces.

Assume that the problem is unforced ($\widetilde{\mathbf f}=\mathbf 0$ and $\widetilde{\mathbf g}|_N=\mathbf 0$) and the Dirichlet boundary conditions are static $(\dot{\mathbf g}_D=\mathbf 0$). Then the solution to the HDG-semidiscrete equations with flux defined by satisfy $$\frac{\mathrm d}{\mathrm dt}
\left(\frac12 \| \widetilde{{\boldsymbol}\sigma}_h(t)\|_{\mathcal A}^2+\frac12 \| \dot{{\mathbf}u}_h(t)\|_\rho^2 \right)
=\pm \| {{\mathbf P_M}}\dot{{\mathbf}u}_h(t)-\dot{\widehat{{\mathbf}u}}_h(t)\|_\tau^2 \qquad \forall t.$$ The solution to the HDG equations with flux defined by satisfy $$\frac{\mathrm d}{\mathrm dt}
\left(\frac12 \| \widetilde{{\boldsymbol}\sigma}_h(t)\|_{\mathcal A}^2+\frac12 \| \dot{{\mathbf}u}_h(t)\|_\rho^2 
+\frac12 \|{{\mathbf P_M}}{\mathbf}u_h(t)-\widehat{{\mathbf}u}_h(t)\|_\tau^2\right)=0 \qquad \forall t.$$

It follows from a simple argument: (a) differentiate and test with $\widetilde{{\boldsymbol}\sigma}_h(t)$, (b) test with $\dot{{\mathbf}u}_h(t)$, (c) test with $\dot{\widehat{{\mathbf}u}}_h(t)$; finally add the result of (a)-(c) using the fact that $\widehat{{\mathbf}u}_h(t)$ is constant in time on the Dirichlet faces.

Numerical experiments
=====================

In order to have a good collection of meshes of the domain, we will conduct all the experiments in the unit cube $\Omega=(0,1)^3$. The cube will be partitioned in a sequence of tetrahedrizations obtained as follows: we first divide the cube along the coordinate planes into $n^3$ equally sized cubes (for $n=1,\ldots,7$), and then partition the cubes into six tetrahedra. We will tag the triangulation as ${{\mathcal T_h}}$, where we can consider that $h=1/n$. The partitions corresponding to $n_1$ and $n_2$, where $n_2$ is a multiple of $n_1$, are nested.

Tests will be run for constant and variable coefficients. In all cases we will consider isotropic materials, with $\mathcal C {\boldsymbol}\varepsilon=2\mu{\boldsymbol}\varepsilon+\lambda({\mathrm{tr}\,}{\boldsymbol}\varepsilon) I_{3\times3}$, for constant or variable Lamé parameters $\lambda$ and $\mu$. In the case of constant coefficients we will take $\rho=\lambda=\mu=1$. In the variable case we take $$\begin{aligned}
{6}\label{eq:coeffsvar}
& \rho(\mathbf x) = 1+|\mathbf x|^2,\\
& \lambda(\mathbf x) = 2 + 0.2x_1^2 + 0.3x_2^2 + 0.04x_3^2,\\
& \mu(\mathbf x)= 3 + 0.5x_2^2 + 0.03x_3^2.\end{aligned}$$

The code expands the three dimensional HDG Matlab package of [@FuGaSa:2013]. The code is written for the second-order-in-frequency form given in , which allows for easy comparisons with the physical quantities of interest and for almost straightforward changes to have all the methods we have considered in this paper. All the element-by-element operations (needed for building the local solvers and recovering displacement and stress after the hybridized system has been solved) are carried out in parallel.

#### A problem with variable coefficients. {#sec:7.2}

We imbue the domain with the variable density and material coefficients given by . Dirichlet boundary conditions are imposed on the top and bottom faces of the cube, while Neumann conditions are given on the remaining faces of $\Omega$. Data are built so that $$\mathbf u(\mathbf x) = 
    \begin{bmatrix}
        \cos(\pi x_1) \sin (\pi x_2) \cos (\pi x_3) \\
        5x_1^2x_2x_3+4x_1x_2x_3+3x_1x_2x_3^2 + 17       \\
        \cos(2x_2)\cos(3x_2)\cos(x_3)
    \end{bmatrix}$$ is the exact solution. Figure \[fig:test1-fo\] shows approximation errors for the method that is equivalent to the one applied to first-order-in-frequency formulation, i.e., we set $\alpha_\kappa := \imath \kappa, \tau = h^{-1}$ in . We measure errors for the displacement and the ‘physical’ stress $$\label{eq:errors}
\|\mathbf u-\mathbf u_h\|_\Omega, 
    \qquad
\|\widetilde{{\boldsymbol}\sigma}-\widetilde{{\boldsymbol}\sigma}_h\|_\Omega.$$ Error plots in Figure \[fig:test1-fo\] are shown for the method for different values of the polynomial degree (recall that ${\boldsymbol}\sigma$ is approximated with piecewise $\mathcal P_k$ functions, while $\mathbf u$ is approximated with piecewise $\mathcal P_{k+1}$ functions). The lowest order method ($k=0$), for which we do not have theoretical support shows a quite bad performance in practice and the corresponding results are not shown.

\
\

We next explore the method based on the pure second-order in frequency formulation ($\alpha_\kappa = 1, \tau = h^{-1}$ in ). Approximation errors for several values of the polynomial degree $k$ are shown in Figure \[fig:test1-so\].

\
\

#### Plane wave solutions. {#pwave}

We now test a problem with constant coefficients and plane-wave solutions as exact data, with Dirichlet conditions imposed on the entire boundary. Plane waves for the time-harmonic elastic wave equation are functions of the form: $$\mathbf u(\mathbf x) = a\, \mathbf e \exp \left( (-\imath\kappa/c) \mathbf x \cdot \mathbf d \right),$$ where $\mathbf d$ (direction of propagation) and $\mathbf e$ (direction of elastic displacement) are unit vectors, $a$ is a constant and: (a) either $\mathbf d = \mathbf e$ and the wavespeed is $c = \sqrt{(\lambda + 2 \mu)/\rho}$ (pressure wave), or (b) $\mathbf d \cdot \mathbf e = 0$ and $c = \sqrt{\mu/\rho}$ (shear wave).

We use amplitude $a=0.3$ and frequency $\kappa = 1$ and tetrahedrizations of the cube each with $6n^3$ elements. Convergence results for pressure and shear waves with the first-order-in-frequency formulation ($\alpha_\kappa = i\kappa, \tau = h^{-1}$ in ) are shown in Figure \[fig:pressurewave-1\]. Both show the same pattern of optimal order $\mathcal O(h^{k+2})$ for the displacement variable and $\mathcal O(h^{k+1})$ for the stress variable.

\
\

#### Increasing the frequency $\kappa$. {#sec:7.4}

In a final test, we look for the behavior of the method as the frequency increases $h\to 0$ while $h \kappa$ remains constant. We will take $\kappa$-dependent plane waves as exact solutions. In order to play with different values of $\kappa$ without hitting elastic modes, we impose an impedance boundary condition $$\widetilde{{\boldsymbol}\sigma}\,\mathbf n -\imath \kappa \mathbf u =\widetilde{\mathbf g}_R \qquad \mbox{on $\Gamma$}.$$ Note that we have not carried out the analysis for these boundary conditions, but that convergence estimates can be obtained with uncomplicated modifications of what we have done in the previous sections. The HDG equation corresponding to the boundary conditions and interelement stress balance become $$\langle\widehat{{\boldsymbol}\sigma}_h\mathbf n - {\boldsymbol}\tau({{\mathbf P_M}}\mathbf u_h - \widehat{\mathbf u}_h),{\boldsymbol}\mu\rangle_{\partial{{\mathcal T_h}}}-\imath\kappa\langle \widehat{\mathbf u}_h, {\boldsymbol}\mu \rangle_\Gamma
= \langle \widetilde{\mathbf g}_R, {\boldsymbol}\mu \rangle_\Gamma \qquad \forall {\boldsymbol}\mu \in {{\mathbf M_h}}.$$ This is just a slight perturbation of the pure Neumann problem, so the linear system (with coefficient representations as in Section \[sec:6\]) may be represented by $$\left[\begin{array}{ccc}
    \mathrm A & \mathrm D^\top & -\mathrm N^\top \\
    -\mathrm D & -\kappa^2\mathrm M-\imath\kappa\mathrm T_{11} & \imath\kappa\mathrm T_{12} \\
    \mathrm N & -\imath\kappa\mathrm T_{12}^\top & \imath\kappa\mathrm T_{22} -\imath\kappa \mathrm C
\end{array}\right]
\left[\begin{array}{l} 
    \widetilde{\underline\sigma} \\ \underline u\\ \underline{\widehat u}
\end{array}\right]
=
\left[\begin{array}{c}
    \underline 0 \\ \widetilde{\underline f} \\ \widetilde{\underline g}_R
\end{array}\right],$$ where $\mathrm C$ is a boundary mass matrix related to the space ${{\mathbf M_h}}$ restricted to $\Gamma$. We are going to plot relative errors $$\|\mathbf u-\mathbf u_h\|_\Omega/\|\mathbf u\|_\Omega,
\qquad
\|\widetilde{{\boldsymbol}\sigma}-
\widetilde{{\boldsymbol}\sigma}_h\|_\Omega/\|\widetilde{{\boldsymbol}\sigma}\|_\Omega,$$ for a given exact plane pressure wave solution. We will fix the value $h\kappa$ and use finer grids on the cube, thus automatically increasing the frequency. Note that for plane waves $\| \mathbf u\|_\Omega=\mathcal O(1)$ but $\|\widetilde{{\boldsymbol}\sigma}\|_\Omega=\mathcal O(\kappa)=\mathcal O(h^{-1})$. We are expecting errors to stay bounded as $h\to 0$ (notice that we do not have a theory supporting this and we do not go very far with the frequency). Some experiments are reported in Figure \[fig:hk-constant-a\], using the lowest order method ($\mathcal P_1$ approximation for the stress, $\mathcal P_2$ for the displacement) with different values of $h\kappa$. That the error for the displacement seems to be declining as $h\to 0$, which shows the numerical method doing somewhat better than expected for the displacement.

\
\
\

#### Some conclusions.

We have introduced and analyzed a family of Hybridizable Discontinuous Galerkin methods for the equations of time-harmonic linear elastodynamics in three dimensions. Convergence orders are shown to be optimal in both variables (displacement and stress), even if different polynomial degrees are used for their approximation, which means that the lower approximation order of the stress does not pollute the approximation for the displacement. We have transferred the methods to an HDG semidiscretization in space of the transient model equations and shown that one of the choices leads to a conservative method. The experiments show that the methods behave well for fixed moderate frequences and that for increasing frequencies the scheme has good asymptotic properties.

Some additional proofs {#sec:A}
======================

The dual adjoint problem is equivalent to

\[eq:Y.1\] $$\begin{aligned}
{6}
\label{eq:Y.1a}
\nabla\cdot\mathcal C{\boldsymbol}\varepsilon({\boldsymbol}\phi) + \kappa^2\rho\,{\boldsymbol}\phi
     &= \imath\kappa {{\mathbf e_h^u}}&\qquad &\mbox{in $\Omega$},\\
{\boldsymbol}\phi &={\boldsymbol}0 & & \mbox{on $\Gamma_D$},\\
\mathcal C{\boldsymbol}\varepsilon({\boldsymbol}\phi)\,\mathbf n &={\boldsymbol}0 & & \mbox{on $\Gamma_N$},\end{aligned}$$

and therefore $\|{\boldsymbol}\phi\|_{1,\Omega} \le \kappa C_\kappa \| {{\mathbf e_h^u}}\|_\Omega$. From a variational formulation of we can prove the identity $$\| \kappa {\boldsymbol}\phi\|_\rho^2=\imath(\kappa {{\mathbf e_h^u}}, \overline{{\boldsymbol}\phi})_\Omega+
(\mathcal C{\boldsymbol}\varepsilon({\boldsymbol}\phi),{\boldsymbol}\varepsilon(\overline{{\boldsymbol}\phi}))_\Omega,$$ and therefore, using Young’s inequality and hiding constants only related to the equation’s coefficients, we can prove that $$\|\kappa {\boldsymbol}\phi\|_\Omega \lesssim (1+C_\kappa\kappa) \| {{\mathbf e_h^u}}\|_\Omega.$$ We now rewrite as $$\nabla\cdot\mathcal C{\boldsymbol}\varepsilon({\boldsymbol}\phi)-\rho{\boldsymbol}\phi=\imath\kappa{{\mathbf e_h^u}}-(\kappa^2+1)\rho{\boldsymbol}\phi$$ and use the regularity estimate to bound $\|{\boldsymbol}\phi\|_{2,\Omega}$.

The first order in space, second order in frequency system is

$$\begin{aligned}
{6}
\mathcal A\widetilde{{\boldsymbol}\sigma} -{\boldsymbol}\varepsilon({\mathbf}u) &=\mathbf 0
    &\qquad &\mbox{in $\Omega$},\\
\nabla\cdot\widetilde{{\boldsymbol}\sigma} +\kappa^2\,\rho\,{\mathbf}u &=\widetilde{{\mathbf}f}
    &\qquad &\mbox{in $\Omega$},\\
{\mathbf}u &=\mathbf g_D
    &\qquad &\mbox{on $\Gamma_D$},\\
\widetilde{{\boldsymbol}\sigma}\mathbf n &=\widetilde{\mathbf g}_N
    &\qquad &\mbox{on $\Gamma_N$}.\end{aligned}$$

From this moment on, we will drop all tildes in the formulas. It has to be understood though that the stress that we are computing with this method is the physical stress and not the one scaled by $\imath/\kappa$. The error equations are

\[eq:X.3\] $$\begin{aligned}
{6}
 (\mathcal A {{e_h^\sigma}},{\boldsymbol}\xi)_{{\mathcal T_h}}+({{\mathbf e_h^u}},\nabla\cdot{\boldsymbol}\xi)_{{\mathcal T_h}}-\langle {{\widehat{\mathbf e}_h^u}},{\boldsymbol}\xi\mathbf n\rangle_{\partial{{\mathcal T_h}}}
    &= (\mathcal A {{\varepsilon_h^\sigma}},{\boldsymbol}\xi)_{{\mathcal T_h}},\\
\nonumber
-(\nabla\cdot{{e_h^\sigma}},\mathbf w)_{{\mathcal T_h}}-\kappa^2 (\rho{{\mathbf e_h^u}},\mathbf w)_{{\mathcal T_h}}\\
    +\langle{\boldsymbol}\tau({{\mathbf P_M}}{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}}),{{\mathbf P_M}}\mathbf w\rangle_{\partial{{\mathcal T_h}}} 
    &=-\kappa^2(\rho{{{\boldsymbol}\varepsilon_h^u}},\mathbf w)_{{\mathcal T_h}}\\
\nonumber
    &\phantom{=}-\langle {{\varepsilon_h^\sigma}}\mathbf n,\mathbf w\rangle_{\partial{{\mathcal T_h}}}
        +\langle{\boldsymbol}\tau{{{\boldsymbol}\varepsilon_h^u}},{{\mathbf P_M}}\mathbf w\rangle_{\partial{{\mathcal T_h}}},\\
\langle {{e_h^\sigma}}\mathbf n-{\boldsymbol}\tau({{\mathbf P_M}}{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}}),{\boldsymbol}\mu\rangle_{\partial{{\mathcal T_h}}\setminus\Gamma_D}
    &=\langle{{\varepsilon_h^\sigma}}\mathbf n,{\boldsymbol}\mu\rangle_{\partial{{\mathcal T_h}}\setminus\Gamma_D}
        -\langle{\boldsymbol}\tau {{\mathbf P_M}}{{{\boldsymbol}\varepsilon_h^u}},{\boldsymbol}\mu\rangle_{\partial{{\mathcal T_h}}\setminus\Gamma_D}\\
\langle{{\widehat{\mathbf e}_h^u}},{\boldsymbol}\mu\rangle_{\Gamma_D} &=0,\end{aligned}$$

and a simple argument shows the new energy identity $$\begin{aligned}
{6}
\label{eq:X.0}
& \|{{e_h^\sigma}}\|_{\mathcal A}^2-\kappa^2\|{{\mathbf e_h^u}}\|_\rho^2+\|{{\mathbf P_M}}{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}}\|_\tau^2 \\
& \hspace{5pt} = (\mathcal A{{\varepsilon_h^\sigma}},\overline{{e_h^\sigma}})_{{\mathcal T_h}}-\kappa^2(\rho \overline{{{\boldsymbol}\varepsilon_h^u}},{{\mathbf e_h^u}})_{{\mathcal T_h}}-\langle\overline{{\varepsilon_h^\sigma}}\mathbf n,{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}}\rangle_{\partial{{\mathcal T_h}}}
    +\langle{\boldsymbol}\tau \, \overline{{{\boldsymbol}\varepsilon_h^u}}, {{\mathbf P_M}}{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}}\rangle_{\partial{{\mathcal T_h}}}.
    \nonumber\end{aligned}$$ The adjoint problem has to be written as

\[eq:X.1\] $$\begin{aligned}
{6}
\mathcal A{\boldsymbol}\psi +{\boldsymbol}\varepsilon({\boldsymbol}\phi) &=\mathbf 0
    &\qquad &\mbox{in $\Omega$},\\
-\nabla\cdot{\boldsymbol}\psi +\kappa^2\,\rho\,{\boldsymbol}\phi &= -{{\mathbf e_h^u}}&\qquad &\mbox{in $\Omega$},\\
{\boldsymbol}\phi &=\mathbf 0
    &\qquad &\mbox{on $\Gamma_D$},\\
{\boldsymbol}\psi\mathbf n &=\mathbf 0
    &\qquad &\mbox{on $\Gamma_N$},\end{aligned}$$

where the negative sign in the right hand side is added for convenience. With the usual regularity hypotheses, the scaled regularity inequalities for the solution of this problem are: $$\kappa(\| {\boldsymbol}\phi\|_{1,\Omega} + \| \rho{\boldsymbol}\phi\|_{1,{{\mathcal T_h}}})+\| \mathcal A {\boldsymbol}\psi\|_{1,\Omega}+ \|{\boldsymbol}\psi\|_{1,{{\mathcal T_h}}}
+ \|{\boldsymbol}\phi\|_{2,\Omega} \le E_\kappa\|{{\mathbf e_h^u}}\|_\Omega,$$ with $E_\kappa$ bounded as in the statement of the theorem. The proof of this inequality is very similar to the proof of above. After integration by parts and introduction of projections it can be shown that the solution of satisfies the following identities

\[eq:X.4\] $$\begin{aligned}
{6}
(\mathcal A{{e_h^\sigma}},\overline{{\boldsymbol}\psi})_{{\mathcal T_h}}-(\nabla\cdot{{e_h^\sigma}},{{\Pi_W}}\overline{{\boldsymbol}\phi})_{{\mathcal T_h}}+\langle{{e_h^\sigma}}\mathbf n,\overline{{\boldsymbol}\phi}\rangle_{\partial{{\mathcal T_h}}} &=0,\\
({{\mathbf e_h^u}},\nabla\cdot{{\Pi_V}}\overline{{\boldsymbol}\psi})_{{\mathcal T_h}}+\langle {{\mathbf e_h^u}},({{\Pi_V}}\overline{{\boldsymbol}\psi}-\overline{{\boldsymbol}\psi})\mathbf n\rangle_{\partial{{\mathcal T_h}}}
    -\kappa^2 (\rho\,{{\mathbf e_h^u}},\overline{{\boldsymbol}\phi})_{{\mathcal T_h}}& = \| {{\mathbf e_h^u}}\|_\Omega^2,\\
\langle {{\widehat{\mathbf e}_h^u}},\overline{{\boldsymbol}\psi}\mathbf n\rangle_{\partial{{\mathcal T_h}}} &=0.\end{aligned}$$

Testing the error equations with the conjugates of the projections of the adjoint problem and rearranging terms, we prove

\[eq:X.5\] $$\begin{aligned}
{6}
 (\mathcal A {{e_h^\sigma}},\overline{{\boldsymbol}\psi})_{{\mathcal T_h}}+({{\mathbf e_h^u}},\nabla\cdot {{\Pi_V}}\overline{{\boldsymbol}\psi})_{{\mathcal T_h}}-\langle {{\widehat{\mathbf e}_h^u}},\overline{{\boldsymbol}\psi}\mathbf n\rangle_{\partial{{\mathcal T_h}}}
    &=\ell_1({\boldsymbol}\psi),\\
-(\nabla\cdot{{e_h^\sigma}},{{\Pi_W}}\overline{{\boldsymbol}\phi})_{{\mathcal T_h}}-\kappa^2 (\rho{{\mathbf e_h^u}},\overline{{\boldsymbol}\phi})_{{\mathcal T_h}}&=\ell_2({\boldsymbol}\phi),\\
\langle {{e_h^\sigma}}\mathbf n,\overline{{\boldsymbol}\phi}\rangle_{\partial{{\mathcal T_h}}}
    &=\ell_3({\boldsymbol}\phi),\end{aligned}$$

where $$\begin{aligned}
\ell_1({\boldsymbol}\psi)&:=&
    (\mathcal A{{e_h^\sigma}},\overline{{\boldsymbol}\psi}-{{\Pi_V}}\overline{{\boldsymbol}\psi})_{{\mathcal T_h}}+(\mathcal A {{\varepsilon_h^\sigma}},{{\Pi_V}}\overline{{\boldsymbol}\psi})_{{\mathcal T_h}}-\langle{{\widehat{\mathbf e}_h^u}},(\overline{{\boldsymbol}\psi}-{{\Pi_V}}\overline{{\boldsymbol}\psi})\mathbf n\rangle_{\partial{{\mathcal T_h}}},\\
\ell_2({\boldsymbol}\phi) &:=& -\kappa^2\left(
    (\rho{{\mathbf e_h^u}},\overline{{\boldsymbol}\phi}-{{\Pi_W}}\overline{{\boldsymbol}\phi})_{{\mathcal T_h}}+(\rho{{{\boldsymbol}\varepsilon_h^u}},{{\Pi_W}}\overline{{\boldsymbol}\phi})_{{\mathcal T_h}}\right)\\
     & & -\langle {{\varepsilon_h^\sigma}}\mathbf n,{{\Pi_W}}\overline{{\boldsymbol}\phi}\rangle_{\partial{{\mathcal T_h}}}
        +\langle{\boldsymbol}\tau{{\mathbf P_M}}{{{\boldsymbol}\varepsilon_h^u}},{{\Pi_W}}\overline{{\boldsymbol}\phi}\rangle_{\partial{{\mathcal T_h}}}
         -\langle{\boldsymbol}\tau({{\mathbf P_M}}{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}}),{{\Pi_W}}\overline{{\boldsymbol}\phi}\rangle_{\partial{{\mathcal T_h}}} \\
\ell_3({\boldsymbol}\phi) &:=& \langle{{\varepsilon_h^\sigma}}\mathbf n,{{\mathbf P_M}}\overline{{\boldsymbol}\phi}\rangle_{\partial{{\mathcal T_h}}}
        -\langle{\boldsymbol}\tau {{\mathbf P_M}}{{{\boldsymbol}\varepsilon_h^u}},\overline{{\boldsymbol}\phi}\rangle_{\partial{{\mathcal T_h}}}
        +\langle {\boldsymbol}\tau({{\mathbf P_M}}{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}}),\overline{{\boldsymbol}\phi}\rangle_{\partial{{\mathcal T_h}}}.\end{aligned}$$ The sum of equations can then be compared with the sum of equations to prove the duality identity: $$\begin{aligned}
\label{eq:X.6}
\| {{\mathbf e_h^u}}\|_\Omega^2
    &=& ( (\mathcal A{{e_h^\sigma}},\overline{{\boldsymbol}\psi}-{{\Pi_V}}\overline{{\boldsymbol}\psi})_{{\mathcal T_h}}+(\mathcal A{{\varepsilon_h^\sigma}},{{\Pi_V}}\overline{{\boldsymbol}\psi})_{{\mathcal T_h}}\\
    \nonumber
    & & -\kappa^2 \left((\rho\,{{\mathbf e_h^u}},\overline{{\boldsymbol}\phi}-{{\Pi_W}}\overline{{\boldsymbol}\phi})_{{\mathcal T_h}}+(\rho\,{{{\boldsymbol}\varepsilon_h^u}},{{\Pi_W}}\overline{{\boldsymbol}\phi})_{{\mathcal T_h}}\right)\\
    \nonumber
    & & +\langle {{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}},(\overline{{\boldsymbol}\psi}-{{\Pi_V}}\overline{{\boldsymbol}\psi})\mathbf n\rangle_{\partial{{\mathcal T_h}}}\\
    \nonumber
    & & +\langle{\boldsymbol}\tau({{\mathbf P_M}}{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}})-{\boldsymbol}\tau{{\mathbf P_M}}{{{\boldsymbol}\varepsilon_h^u}},\overline{{\boldsymbol}\phi}-{{\Pi_W}}\overline{{\boldsymbol}\phi}\rangle_{\partial{{\mathcal T_h}}}
           - \langle {{\varepsilon_h^\sigma}}\mathbf n,{{\Pi_W}}\overline{{\boldsymbol}\phi}-{{\mathbf P_M}}\overline{{\boldsymbol}\phi}\rangle_{\partial{{\mathcal T_h}}}.\end{aligned}$$ What is left now is the proof of bounds for the right-hand sides of and . This process requires just going carefully over the proofs of Proposition and . Nothing essential is changed. We can write the results with our shorthand notation for errors $\Sigma:=\| {{e_h^\sigma}}\|_{\mathcal A}$, $\mathrm T:=\| {{\mathbf P_M}}{{\mathbf e_h^u}}-{{\widehat{\mathbf e}_h^u}}\|_\tau$, $\mathrm U:=\| {{\mathbf e_h^u}}\|_\Omega$, and approximation terms $\Sigma_h:= h^t |{\boldsymbol}\sigma|_{t,\Omega}$, $\mathrm U_h:=h^{s-1} |\mathbf u|_{s,\Omega}$. The bounds we obtain are: $$\begin{aligned}
\Sigma^2+ \mathrm T^2
    & \lesssim \Sigma\,\Sigma_h + \kappa^2 \mathrm U\,\mathrm U_h 
    + (\Sigma_h+\mathrm U_h) \mathrm T +\Sigma_h^2+\kappa^2\mathrm U^2, \\
\mathrm U & \lesssim \alpha^2 (\Sigma+\Sigma_h+\mathrm U_h +\mathrm T),\end{aligned}$$ where $\alpha:=E_\kappa h (1+\kappa)$. The condition that allows us to bootstrap is $C(\alpha\kappa)^2 \le 1/4$, where $C$ is a constant related to the constants hidden in the symbols $\lesssim$ above. After simplification, we prove $$\Sigma+\mathrm T  \lesssim \Sigma_h+(1+\kappa) \mathrm U_h \qquad 
\mathrm U  \lesssim \alpha (\Sigma_h+(1+\kappa)\mathrm U_h),$$ which is the statement of the theorem.

[^1]: Part of this work was developed while AH and FJS partially funded by NSF grant DMS 1216356.
