# Numerical investigation of stabilization in the Hybridizable Discontinuous Galerkin method for linear anisotropic elastic equation

Ha Pham, Florian Faucher, Hélène Barucq

► To cite this version:

Ha Pham, Florian Faucher, Hélène Barucq. Numerical investigation of stabilization in the Hybridizable Discontinuous Galerkin method for linear anisotropic elastic equation. Computer Methods in Applied Mechanics and Engineering, 2024, 428, pp.117080. 10.1016/j.cma.2024.117080. hal-04503407

HAL Id: hal-04503407

https://hal.science/hal-04503407v1

Submitted on 6 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Numerical investigation of stabilization in the Hybridizable Discontinuous Galerkin method for linear anisotropic elastic equation

Ha Pham*

Florian Faucher*

Hélène Barucq*

April 29, 2024

Abstract

This work is concerned with implementing the hybridizable discontinuous Galerkin (HDG) method to solve the linear anisotropic elastic equation in the frequency domain. First-order formulation with the compliance tensor and Voigt notation are employed to provide a compact description of the discretized problem and flexibility with highly heterogeneous media. We further focus on the question of optimal choices of stabilization in the definition of HDG numerical traces. For this purpose, we construct a hybridized Godunov-upwind flux for anisotropic elastic media possessing three distinct wavespeeds. This stabilization removes the need to choose a scaling factor, contrary to the identity and Kelvin–Christoffel based stabilizations which are popular choices in the literature. We carry out comparisons among these families for isotropic and anisotropic material, with constant background and highly heterogeneous ones, in two and three dimensions. These experiments establish the optimality of the Godunov stabilization which can be used as a reference choice for a generic material in which different types of waves propagate.

# 1 Introduction

In this work, we implement the Hybridizable Discontinuous Galerkin (HDG) method to solve the linear anisotropic elastic equation in the frequency domain and construct robust stabilization in the definition of the numerical traces required by the method. The investigation has in mind applications which are geared towards wave simulation and full waveform inversion (FWI) in heterogeneous Earth subsurfaces. To discretize with the HDG method, the elastic equation with its constitutive law is reformulated as a first-order system with a pair of unknowns consisting of displacement **u** and stress **σ**. This system is then written as a union of local problems defined on each mesh element with prescribed data given by the HDG global variable **λ**_u representing the trace of **u**. The local problems are linked by transmission conditions and boundary conditions defined on the interfaces and exterior boundary, respectively. Together with a definition of the numerical trace of the traction denoted by $\hat{\sigma}$ on the boundary of each element, the transmission conditions combined with the inversion of local problems result in a global problem in **λ**_u defined on the skeleton of the mesh. A choice of stabilization operator denoted by $\tau_u$ is required to define $\hat{\sigma}$, and the judicial choice of which is important to ensure the robustness and accuracy of the discretization method.

<footnote>*Project-Team Makutu, Inria, University of Pau and Pays de l'Adour, TotalEnergies, CNRS UMR 5142</footnote>
Works implementing HDG method in linear elasticity are e.g., [4] for time-harmonic elasticity, [32, 23] for static isotropic elasticity, [35, 22, 15, 17] for static elasticity including anisotropy, and [36, 30, 21] for time-domain isotropic elasticity. These works, including ours, follow the approach of the original HDG method devised by Cockburn and co-authors for second-order elliptic equations [13, 14]. The HDG method retains the attractive features of the Discontinuous Galerkin family, such as being locally conservative, flexible with h-p adaptivity strategy and complex geometry [21]. Additionally, it provides optimal convergence with $L^2$-based approximation with equal degree of polynomials for the primal variable and its derivatives. Finally, it is highly parallelizable and amenable to hybrid computing architecture, cf. [18]. For more detailed discussions of the method, we refer to [1, 14, 10, 35, 32, 23, 11], and for its implementation to various equations, [21, 24]. Other methods related to HDG are Weak Galerkin or Hybrid High Order, cf. [42] and [16, 12] for elasticity.

In view of applications in seismology, in particular with FWI, strong features of the HDG method are its flexibility with material heterogeneity and unstructured meshes, and that it can provide accurate and stable performance at tenable computational costs. The Finite Difference method can outperform Finite Element methods on media with low-variation properties; however, its need for structured grids creates problems for backgrounds with strong variation and topography (roughness of boundary), and with implementation of boundary conditions, cf. [40, Section 4.2]. In our work, the above advantages of the HDG method are further exploited by employing a first-order formulation in $(\mathbf{u}, \boldsymbol{\sigma})$ with the constitutive law written in terms of the compliance tensor. This allows for parameter variation within a mesh cell and, as shown in our numerical investigation, leads to a better representation of complex media and a more numerically stable resolution for these media. Another feature of our implementation is that the discretized problem is compactly described using Voigt representation, which maximizes matrix operations and provides an efficient book-keeping of physical parameters for any generic anisotropic material. These features are important not only in code building for the forward but also for the inverse solver.

In this work, we will focus on the question of choosing a robust choice for the stabilization matrix $\boldsymbol{\tau}_{\mathbf{u}}$. Most of the above-cited references which employ the HDG method for elasticity work with an identity-based stabilization, i.e., choosing $\boldsymbol{\tau}_{\mathbf{u}} = \tau \mathbb{I}$ with $\mathbb{I}$ the identity matrix and $\tau$ a scalar scaling factor. This is with the exception of [35, 22, 15] which investigate stabilization of the form $\boldsymbol{\tau}_{\mathbf{u}} = \tau \boldsymbol{\Gamma}$ where $\boldsymbol{\Gamma}$ is the Kelvin-Christoffel (KC) matrix for elastostatics, and [36] which constructs the hybridized Godunov-upwind flux for time-domain isotropic elasticity. For static isotropic elasticity, numerical investigations on choices for the scaling factor $\tau$ are carried out for identity-based stabilizations in [32], and in [35, 22] for KC-based stabilizations. We also refer to [26] for a study of scaling factors for identity-based stabilizations for the Helmholtz equation considering complex-valued scaling factors. On the other hand, there are no such investigations for elastodynamics, with isotropy or anisotropy.

We start our investigation by extending the hybridization technique of upwind fluxes for isotropy in [36] to anisotropy under the assumption of three distinct wave speeds. The ingredients for its construction in anisotropy, which remain the same as in isotropy, consist of the exact solutions of the Riemann problem and their Rankine-Hugoniot jump conditions. They are also employed in [43, 39, 38, 44] to construct upwind flux for DG implementation in the time domain, which depend however on information from both sides of an interface, in contrast with the hybridized HDG fluxes in our work and [36], which are defined from only one-sided data. Hybridization of the Godunov flux is also discussed in a series of works with different formalisms [5, 6, 7] and in which the Rankine-Hugoniot jump conditions are recognized in sequel works
([6, 7]) to provide a more natural and direct way to devise the HDG flux. Also employing the notion of Riemann problems to devise the HDG scheme for compressible fluid flow are, e.g., [41, 29] and the references therein, in which hybridization is carried out with Lax-Friedrichs, Roe and Harten-Lax-Van Leer approximate Riemann solvers.

In the second part of this work, we compare three families of stabilization: identity-based, KC-based, and hybridized Godunov (cf. Subsection 3.4) for various settings and types of waves. We show that the Godunov stabilization as derived is optimal, while the other two families require correct scaling factors which might not be universal in order to match the performance of the Godunov stabilization. The comparison is also carried out for setting supporting mixed types of wavefields in highly varying backgrounds, to conclude on the versatility of the Godunov stabilization and the fact that it provides a reference choice for generic anisotropic material and general simulations.

The organization of this article is as follows. We gather in Section 2 the notations needed in later sections, in particular for the discretization and in the construction of the hybridized flux. In Section 3, a first-order formulation for the elastic equation, its HDG formulation, and corresponding discretized problem are given. The construction for elastic anisotropy is carried out in Section 4, and numerical experiments are displayed in Section 5.

## 2 Notations

We introduce here the quantities that appear in the rest of the paper. In particular, the quantities with Voigt notation (Subsections 2.2 and 2.4) are employed to describe the discretized problem in Subsection 3.3, and to construct the hybridized Godunov-upwind flux in Subsection 3.3.

### 2.1 General notations

We denote by $S_2$ the set of 3 × 3 symmetric matrices and $S_4$ the set of fourth order tensors bearing the symmetry of the elastic tensor:

$$ 
\begin{aligned}
S_2 &:= \{ \mathbf{w} = (w_{ij})_{i,j=x,y,z} \in \mathbb{C}^{3 \times 3} \mid w_{ij} = w_{ji} \}, \\
S_4 &:= \{ \mathbf{T} = (t_{ijkl})_{i,j,k,l=x,y,z} \in \mathbb{C}^{3 \times 3 \times 3} \mid t_{ijkl} = t_{jikl} = t_{ijkl} = t_{klij} \}.
\end{aligned} 
\quad (2.1) 
$$

We list here the operations employed on the following objects: for vectors $\mathbf{v} = (v_i)_{i=1}^3$, $\mathbf{w} = (w_i)_{i=1}^3$, matrices $\mathbf{\chi} = (\chi_{kl})_{k,l=1}^3$, $\tilde{\mathbf{\chi}} = (\tilde{\chi}_{kl})_{k,l=1}^3$, and fourth-order tensor $\mathbf{T} = (t_{ijkl})_{i,j,k,l=1}^3$,

$$ \mathbf{v} \cdot \mathbf{w} = \sum_{j=1}^{3} v_j w_j, \quad (2.2a) \qquad \mathbf{\chi} \mathbf{v} = \sum_{j=1}^{3} \chi_{ij} v_j, \quad (2.2e) $$

$$ \mathbf{\chi} : \tilde{\mathbf{\chi}} = \sum_{i,j=1}^{3} \chi_{ij} \tilde{\chi}_{ij}, \quad (2.2b) \qquad (\mathbf{T}\mathbf{\chi})_{ij} = \sum_{k,l=1}^{3} t_{ijkl} \chi_{kl}, \quad (2.2f) $$

$$ (\mathbf{v} \otimes \mathbf{w})_{ij} = v_i w_j, \quad (2.2c) \qquad \mathbf{v} \odot \mathbf{w} = \frac{\mathbf{v} \otimes \mathbf{w} + \mathbf{w} \otimes \mathbf{v}}{2}, \quad (2.2g) $$

$$ \mathbf{v} \cdot \mathbf{\chi} \cdot \mathbf{w} = \sum_{i,j=1}^{3} v_i \chi_{ij} w_j, \quad (2.2d) \qquad \mathbf{\chi} : \mathbf{S} : \tilde{\mathbf{\chi}} = \sum_{i,j,k,l=1}^{3} \chi_{ij} S_{ijkl} \tilde{\chi}_{kl}. \quad (2.2h) $$
We also employ the following differential operators:

$$ (\nabla \mathbf{w})_{ij} = \partial_j w_i, \qquad \nabla^5 \mathbf{w} = \frac{\nabla \mathbf{w} + (\nabla \mathbf{w})^t}{2}, \qquad (\nabla \cdot \mathbf{\chi})_i = \sum_{j=1}^{3} \partial_j \chi_{ij}. \qquad (2.3) $$

## 2.2 Voigt identification

We define the following mapping $i_{sm}$ which identifies the indices of a symmetric matrix with integer values,

$$ i_{sm}: \begin{array}{l} \text{‘xx’} \mapsto 1, \qquad \text{‘yy’} \mapsto 2, \qquad \text{‘zz’} \mapsto 3, \\ \text{‘zy’, ‘yz’} \mapsto 4, \qquad \text{‘xz’, ‘zx’} \mapsto 5, \qquad \text{‘xy’, ‘yx’} \mapsto 6. \end{array} \qquad (2.4) $$

We denote by $\mathcal{I}_{sm}$ the ordered set consisting of the diagonal and lower-diagonal (equivalently upper-diagonal) indices,

$$ \begin{aligned} \mathcal{I}_{sm} &:= (\text{‘xx’}, \text{‘yy’}, \text{‘zz’}, \text{‘yz’}, \text{‘xz’}, \text{‘xy’}) && \text{identified} \\ &= (\text{‘xx’}, \text{‘yy’}, \text{‘zz’}, \text{‘zy’}, \text{‘zx’}, \text{‘yx’}) && \text{with} \quad (1, \dots, 6). \end{aligned} \qquad (2.5) $$

In Voigt notation, a symmetric 3 × 3 matrix $\mathbf{w}$ is identified with a vector of length 6, denoted by $\vec{\mathbf{w}}$, by employing identification (2.4),

$$ \begin{aligned} S_2 \ni \mathbf{w} = \begin{pmatrix} w_{xx} & w_{xy} & w_{xz} \\ w_{yx} & w_{yy} & w_{yz} \\ w_{zx} & w_{zy} & w_{zz} \end{pmatrix} & \leftarrow\leftarrow_{\text{Voigt identification}} & \vec{\mathbf{w}} = (w_{xx}, w_{yy}, w_{zz}, w_{yz}, w_{xz}, w_{xy})^t \\ & & = (w_{i_{sm}(\mathfrak{J})})_{\mathfrak{J} \in \mathcal{I}_{sm}}. \end{aligned} \qquad (2.6) $$

A fourth order symmetric tensor in $S_4$ is identified with a symmetric matrix of size 6 × 6 as

$$ S_4 \ni \mathbf{T} = (t_{ijkl})^3_{i,j,k,l=1} \quad \underset{\text{Voigt identification}}{\overleftarrow{\mathbf{T}}} = (t_{\mathfrak{J}\mathfrak{J}'}_{\mathfrak{J},\mathfrak{J}' \in \mathcal{I}_{sm}}) = (t_{i_{sm}(\mathfrak{J})i_{sm}(\mathfrak{J}')}_{\mathfrak{J},\mathfrak{J}' \in \mathcal{I}_{sm}}). \qquad (2.7) $$

We will also work with the following identifications,

$$ \vec{\mathbf{w}}^{\dagger} := D^{\dagger} \vec{\mathbf{w}}, \qquad \overline{\mathbf{T}}^{\dagger} = \overline{\mathbf{T}} D^{\dagger}, \qquad \overline{\mathbf{T}}^{\dagger} := D^{\dagger} \overline{\mathbf{T}} D^{\dagger}, \qquad (2.8) $$

$$ \text{with diagonal matrix } D^{\dagger} := \begin{pmatrix} \mathbb{1}_{3 \times 3} & \mathbf{0}_{3 \times 3} \\ \mathbf{0}_{3 \times 3} & 2 \mathbb{1}_{3 \times 3} \end{pmatrix}. \qquad (2.9) $$

The definition in (2.8) and matrix $D^\dagger$ arise when tensor operations are rewritten with Voigt identification. For example, the contraction (2.2b) between two symmetric matrices $\boldsymbol{\tau}_1, \boldsymbol{\tau}_2 \in S_2$ can be written as,

$$ \boldsymbol{\tau}_1 : \boldsymbol{\tau}_2 = \vec{\boldsymbol{\tau}}_1 \cdot \vec{\boldsymbol{\tau}}_2^{\dagger} = \vec{\boldsymbol{\tau}}_1^{\dagger} \cdot \vec{\boldsymbol{\tau}}_2 = \vec{\boldsymbol{\tau}}_1 \cdot D^{\dagger} \cdot \vec{\boldsymbol{\tau}}_2, \qquad (2.10) $$

where the operation in the second and third expression are vector dot products (2.2a), while the fourth expression is the bilinear form defined in (2.2d). The product (2.2f) between a symmetric tensor $\mathbf{C} \in S_4$ and a matrix $\mathbf{\chi} \in S_2$ can be written as,

$$ \overline{\mathbf{C}} \mathbf{\chi} = \overline{\mathbf{C}} D^{\dagger} \overline{\mathbf{\chi}} = \overline{\mathbf{C}}^{\dagger} \overline{\mathbf{\chi}} = \overline{\mathbf{C}} \overline{\mathbf{\chi}}^{\dagger}; \qquad (2.11a) $$

$$ \overline{\mathbf{C}} \overline{\mathbf{\chi}}^{\dagger} = D^{\dagger} \overline{\mathbf{C}} D^{\dagger} \overline{\mathbf{\chi}} \stackrel{(2.8)}{=} \overline{\mathbf{C}}^{\dagger} \overline{\mathbf{\chi}}^{\dagger} = D^{\dagger} \overline{\mathbf{C}} \overline{\mathbf{\chi}}^{\dagger}. \qquad (2.11b) $$

For other identities with Voigt notation, we refer to [31, Section 3].
## 2.3 Useful facts in linear elasticity

We note the following facts which will be employed in the discretization of the elastic equation, cf. (3.1). The stiffness tensor $C = (C_{ijkl})^3_{i,j,k,l=1}$, compliance tensor $S = (S_{ijkl})^3_{i,j,k,l=1}$, the strain matrix $\epsilon = (\epsilon_{IJ})_{I,J=x,y,z}$, and stress matrix $\sigma = (\sigma_{IJ})_{I,J=x,y,z}$, are such that,

$$ C, S \in S_4, \quad C = S^{-1}, \quad \epsilon, \sigma \in S_2. \qquad (2.12) $$

Their relations are given by the constitutive law of linear elasticity,

$$ \sigma = C\epsilon, \quad \text{or equivalently} \quad \epsilon = S\sigma. \qquad (2.13) $$

**Fact 1.** We write the unit vectors in Cartesian basis of $\mathbb{R}^3$ and $\mathbb{R}^6$, respectively as,

$$ e_I \in \mathbb{R}^3, \text{ with } I = x, y, z, \quad \text{and} \quad \hat{e}_J \in \mathbb{R}^6, \text{ with } J \in \mathcal{I}_{\text{sm}}. \qquad (2.14) $$

The latter basis is indexed by $\mathcal{I}_{\text{sm}}$ (2.5). We define the following symmetric matrix using (2.2g),

$$ \bar{\bar{e}}_J := \begin{cases} e_I \odot e_J, & I = J, \\ 2 e_I \odot e_J, & I \neq J. \end{cases} \qquad (2.15) $$

A basis for symmetric matrices is given by $\{\bar{\bar{e}}_J\}$, and in the Voigt identification (2.6), by $\{\hat{e}_J\}$. In particular for the stress tensor, $\sigma = (\sigma_{IJ})^3_{I,J=1} \in S_2$, we have

$$ \sigma = \sum_{J \in \mathcal{I}_{\text{sm}}} \sigma_J \bar{\bar{e}}_J, \quad \vec{\sigma} = \sum_{J \in \mathcal{I}_{\text{sm}}} \sigma_J \hat{e}_J. \qquad (2.16) $$

**Fact 2.** Employing identities in (2.11), the constitutive laws (2.13) can be written with Voigt notations (2.6) and (2.7) as,

$$ \sigma = C\epsilon \Leftrightarrow \vec{\sigma} = \bar{\bar{C}}D^{\dagger}\vec{\epsilon}, \quad \epsilon = S\sigma \Leftrightarrow \vec{\epsilon} = \bar{\bar{S}}D^{\dagger}\vec{\sigma}. \qquad (2.17) $$

The inverse relation between the stiffness tensor $C$ and compliance tensor $S$ can be written as,

$$ \bar{\bar{S}} = (\bar{\bar{C}}^{\dagger})^{-1}, \quad \bar{\bar{S}}^{\dagger} = (\bar{\bar{C}})^{-1}. \qquad (2.18) $$

This is useful for computing the compliance tensor of a generic anisotropic elastic material.

## 2.4 Kelvin–Christoffel (KC) matrix

The 3 × 3 Kevin-Christoffel matrix in direction $\nu$ is defined as, cf. [8],

$$ \Gamma(\nu) := \nu \cdot C \cdot \nu, \quad \text{with} \quad (\nu \cdot C \cdot \nu)_{jk} = \sum_{i,l=1}^{3} \nu_i c_{ijkl} \nu_l. \qquad (2.19) $$

From the symmetry of the stiffness tensor $C \in S_4$, the matrix $\Gamma(\nu)$ is symmetric and positive definite, and thus it is diagonalizable with positive eigenvalues, cf. [31]. The definition of the KC matrix can be rewritten in terms of the Voigt notation of $C$ as follows,

$$ \Gamma(\nu) = \mathcal{A}(\nu)^{\dagger} \bar{\bar{C}}^{\dagger} \mathcal{A}(\nu), \qquad (2.20) $$
where, for a vector $\xi \in \mathbb{C}^3$, we define

$$ \mathcal{A}(\xi) := \sum_{I=x,y,z} \xi_I \mathcal{A}_I, \quad \mathcal{A}^{\dagger}(\xi) := \sum_{I=x,y,z} \xi_I \mathcal{A}_I^{\dagger}, \qquad (2.21) $$

with the elementary matrices

$$ \mathcal{A}_x = \begin{pmatrix} 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}, \quad \mathcal{A}_y = \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}, \quad \mathcal{A}_z = \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \end{pmatrix}. \qquad (2.22) $$

With matrix $D^{\dagger}$ of (2.9), we also define,

$$ \mathcal{A}_I^{\dagger} = \mathcal{A}_I D^{\dagger}, \quad \text{where } I = x, y, z. \qquad (2.23) $$

## 2.5 Examples

**Isotropic elasticity** The stiffness tensor for an isotropic material has the following form in Voigt notation (2.7),

$$ \bar{\mathcal{C}}_{\text{iso}} = \begin{pmatrix} \lambda+2\mu & \lambda & \lambda & 0 & 0 & 0 \\ \lambda & \lambda+2\mu & \lambda & 0 & 0 & 0 \\ \lambda & \lambda & \lambda+2\mu & 0 & 0 & 0 \\ 0 & 0 & 0 & \mu & 0 & 0 \\ 0 & 0 & 0 & 0 & \mu & 0 \\ 0 & 0 & 0 & 0 & 0 & \mu \end{pmatrix}, \quad \bar{\mathcal{S}}_{\text{iso}} = \begin{pmatrix} \frac{1}{E} & -\frac{\nu}{E} & -\frac{\nu}{E} & 0 & 0 & 0 \\ -\frac{\nu}{E} & \frac{1}{E} & -\frac{\nu}{E} & 0 & 0 & 0 \\ -\frac{\nu}{E} & -\frac{\nu}{E} & \frac{1}{E} & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{4\mu} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1}{4\mu} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1}{4\mu} \end{pmatrix}. \qquad (2.24) $$

Here the two degrees of freedom are given by the Lamé parameters $\lambda$ and $\mu$, or by the Young modulus $E$ and Poisson coefficient (ratio),

$$ E := \frac{(3\lambda + 2\mu)\mu}{\lambda + \mu}, \quad \nu := \frac{\lambda}{2(\lambda + \mu)}, \quad \frac{\nu}{E} = \frac{\lambda}{2\mu(3\lambda + 2\mu)}. \qquad (2.25) $$

With $\rho$ denoting the material density, the KC matrix for isotropic elasticity is

$$ \Gamma_{\text{iso}}(\boldsymbol{\nu}) = \mu \mathbf{1} \mathbf{d} + (\lambda + \mu) \boldsymbol{\nu} \otimes \boldsymbol{\nu} = \rho (c_S^2 \mathbf{1} \mathbf{d} + (c_P^2 - c_S^2) \boldsymbol{\nu} \otimes \boldsymbol{\nu}). \qquad (2.26) $$

The second expression is written in terms of the P- and S-wave speed. The isotropic eigenvalues are independent of $\boldsymbol{\nu}$,

$$ \mu = \rho c_S^2 \text{ with multiplicity } 2, \quad \text{and} \quad \lambda + \mu = \rho c_P^2 \text{ with multiplicity } 1. \qquad (2.27) $$

**Vertical transverse isotropy (VTI)** Elastic materials in this family are rotational symmetric around $\mathbf{e}_z$. Their stiffness tensor in Voigt representation (2.7) is given with five independent coefficients, cf. [8, Equation (1.39)],

$$ \bar{\mathcal{C}}_{\text{VTI}} = \begin{pmatrix} C_{11} & C_{11}-2C_{66} & C_{13} & 0 & 0 & 0 \\ C_{11}-2C_{66} & C_{11} & C_{13} & 0 & 0 & 0 \\ C_{13} & C_{13} & C_{33} & 0 & 0 & 0 \\ 0 & 0 & 0 & C_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & C_{44} & 0 \\ 0 & 0 & 0 & 0 & 0 & C_{66} \end{pmatrix}. \qquad (2.28) $$

In this case, the KC matrix takes the following form,

$$ \Gamma_{\text{VTI}}(\boldsymbol{\nu}) = \begin{pmatrix} C_{11} \nu_x^2 + C_{66} \nu_y^2 + C_{44} \nu_z^2 & (C_{13}+C_{44})\nu_y\nu_z & (C_{13}+C_{44})\nu_x\nu_z \\ C_{66}\nu_x^2 + C_{11} \nu_y^2 + C_{44} \nu_z^2 & (C_{11}-C_{66})\nu_x\nu_y & C_{44}(\nu_x^2+\nu_y^2)+C_{33}\nu_z^2 \end{pmatrix}. \qquad (2.29) $$
### 3 HDG formulations for time-harmonic linear elasticity

We consider the propagation of time-harmonic waves in linear elasticity in terms of the displacement $\mathbf{u}$. For a domain $\Omega$ and an interior source $\mathbf{f}$, the wave equations are given by (e.g., [8]),

$$ \text{Formulation } (\mathbf{u}, \boldsymbol{\sigma})_S \quad \begin{cases} -\omega^2 \rho(\mathbf{x}) \mathbf{u}(\mathbf{x}, \omega) - \nabla \cdot \boldsymbol{\sigma}(\mathbf{x}, \omega) = \mathbf{f}(\mathbf{x}, \omega), & (3.1\text{a}) \\ \mathbf{S}(\mathbf{x})\boldsymbol{\sigma}(\mathbf{x}, \omega) = \nabla^s \mathbf{u}(\mathbf{x}, \omega). & (3.1\text{b}) \end{cases} $$

Here, $\boldsymbol{\sigma}$ is the stress tensor, $\rho$ is the density, $\mathbf{S}$ is the compliance tensor which is the inverse of the stiffness tensor $\mathbf{S} := \mathbf{C}^{-1}$.

**Remark 1** (Alternative formulations). Several variants to (3.1) have been introduced and studied in the literature, containing $\omega$ or $\omega^2$, and with possibly different choices of unknowns such as the strain tensor $\boldsymbol{\epsilon} := \nabla^s \mathbf{u}$. For instance, in our report [31], we compare the formulation (3.1) with one in terms of $(\mathbf{v}, \boldsymbol{\sigma})$ with the velocity $\mathbf{v} := -i\omega \mathbf{u}$, resulting in a system containing $\omega$ instead of $\omega^2$. We note that (3.1) gives the formulation for elastostatics as $\omega \to 0$, and we have not observed any difference between these two formulations, cf. [31, Sections 6 and 7].

**Remark 2.** We write the equation in terms of the compliance tensor $\mathbf{S}$ instead of the stiffness tensor $\mathbf{C} = \mathbf{S}^{-1}$. This allows us to easily encode highly-varying backgrounds on meshes containing large cells, see Subsection 5.2.2. In this case, high-order polynomial discretization is to be employed, which is critical for the efficiency of the HDG method compared to CG one, [18]. Secondly, it allows us to have a flexible representation of the physical properties on the discretized domain as seen in our experiments Section 5, and which is a key ingredient for inversion, [20].

#### 3.1 Discretization domain

For numerical resolution, we employ a conforming triangulation of the domain $\Omega$, denoted by $\mathcal{T}_h$, which consists of non-overlapping elements $K^e$, $1 \le e \le |\mathcal{T}_h|$,

$$ \Omega = \bigcup_{e=1}^{|\mathcal{T}_h|} K^e := \Omega_h, \quad \mathcal{T}_h = \{K^e | 1 \le e \le |\mathcal{T}_h|\}. \qquad (3.2) $$

The set of faces in the mesh $\mathcal{T}_h$, denoted by $\Sigma_h$, consists of faces $\mathbf{F}$ which also form the boundary of all elements $K$ in $\mathcal{T}_h$,

$$ \Sigma_h = \bigcup_{e=1}^{|\mathcal{T}_h|} \partial K^e = \{\mathbf{F}^k | 1 \le k \le |\Sigma_h|\}. \qquad (3.3) $$

We distinguish between the interior faces, also called interfaces, $\Sigma_{\text{int}}$ and boundary ones $\Sigma_{\partial}$ which form the boundary of $\Omega_h$:

$$ \Sigma_h = \Sigma_{\text{int}} \cup \Sigma_{\partial}, \quad \text{with } \Sigma_{\partial} = \partial \Omega_h = \Sigma_{\text{N}} \cup \Sigma_{\text{D}} \cup \Sigma_{\infty}. \qquad (3.4) $$

The latter set is partitioned into non-overlapping regions on which Neumann, Dirichlet, and Robin boundary conditions can be imposed. In our experiments, we only consider tetrahedral meshes in 3D and triangle meshes in 2D.
**Local and global indexes of an edge** An edge $F \in \Sigma_h$ can be referred to in two ways,

1. as $F^k$, with $1 \le k \le |\Sigma_h|$, here $k$ is the index of the face $F$ in the ordered set $\Sigma_h$,

2. as $F^{(e,\ell)}$ with $1 \le e \le |\mathcal{T}_h|$, $1 \le \ell \le n_{\text{face}}^e$, when it is considered as part of the boundary of an element $K^e$, and $\ell$ is its index among the set of faces of $K^e$. Note that $n_{\text{face}}^e$ is the number of edges in 2D, and the number of faces in 3D. Working with simplexes, we have,

$$ n_{\text{face}}^e = 4, \text{ in } 3\text{D}, \qquad n_{\text{face}}^e = 3, \text{ in } 2\text{D}. \qquad (3.5) $$

**Jump operator** At an interface $F$, the jumps $\llbracket \cdot \rrbracket$ of a vector $v$ and of the traction of a matrix $\chi$ defined on $\bar{K}^+ \cup \bar{K}^-$ with $\circ$ denoting the interior, are respectively defined as

$$ \llbracket v \rrbracket := v|_{K^+} - v|_{K^-}, \qquad \llbracket \chi \nu \rrbracket := \chi|_{K^+} \nu^+ + \chi|_{K^-} \nu^-, \qquad (3.6) $$

with $\nu^\pm$ the outward-pointing normal vector defined on $\partial K^\pm$.

## 3.2 Statement of HDG problem for formulation $(u, \sigma)_s$

In the HDG method, the original global problem is statically condensed into a problem in terms of a hybrid variable defined on the skeleton of the mesh [28]. Specifically, it is first written as a union of local boundary-valued problems defined for each mesh cell, having as boundary data the hybrid variable. Local problems are linked by transmission conditions (or jump conditions) which constrain weakly the continuity of solutions. In this way, a global problem is obtained and solved in terms of the hybrid variable; after this, the primal and mixed unknowns are found by solving the local problems, which can be realized in a parallel manner. For a more in-depth discussion, we refer the readers to [24].

**Strong form** The local problems are boundary value problems defined on each cell $K \in \mathcal{T}_h$ with prescribed Dirichlet condition $\lambda_u \in L^2(\partial K)$. From (3.1), we have,

$$ - \omega^2 \rho \mathbf{u} - \nabla \cdot \mathbf{\sigma} = \mathbf{f}, \qquad \text{on } K, \qquad (3.7a) $$

$$ \mathbf{S}\mathbf{\sigma} = \nabla^5 \mathbf{u}, \qquad \text{on } K, \qquad (3.7b) $$

$$ \mathbf{u} = \lambda_u, \qquad \text{on } \partial K. \qquad (3.7c) $$

The trace $\lambda_u$ is determined by the transmission condition,

$$ \llbracket \sigma \nu \rrbracket = 0 \quad \text{on } F \in \Sigma_{\text{int}}, \qquad (3.8) $$

and the boundary conditions on the boundary of the domain $\Sigma_\partial$ are given by,

$$ \begin{array}{lll} \lambda_u = \mathbf{g}^D & \text{on } \Sigma_D, & \text{Dirichlet boundary,} \\ \sigma \nu = 0 & \text{on } \Sigma_N, & \text{Neumann boundary,} \\ \sigma \nu = -i\omega\mathcal{Z}\lambda_u & \text{on } \Sigma_\infty, & \text{Robin boundary.} \end{array} \qquad (3.9) $$

The Robin condition employs impedance-like matrix $\mathcal{Z} = (\mathcal{Z}_{IJ})_{I,J=x,y,z}$.
**Finite element spaces** Below, tensor-valued functions whose components taking value in a function space $\mathcal{V}$ are written as,

$$ \text{vector-valued } \mathcal{V}^n, \quad \text{matrix-valued } \mathcal{V}^{n \times n}, \quad \text{symmetric} \quad \text{matrix-valued } \mathcal{V}_{\text{sym}}^{n \times n}. \qquad (3.10) $$

We introduce the global finite element spaces

$$ \begin{align} U_h &= \{ \mathbf{u} = (u_I)_{I=x,y,z} \in L^2(\Omega_h)^3 & : & u_I|_K \in U_h(K), \forall K \in \mathcal{T}_h \}; \nonumber \\ V_h &= \{ \mathbf{v} = (v_{IJ})_{I,J=x,y,z} \in L^2(\Omega_h)^{3 \times 3}_{\text{sym}} & : & v_{IJ}|_K \in V_h(K), \forall K \in \mathcal{T}_h \}; \tag{3.11} \end{align} $$

$$ M_h = \{ \mathbf{w} = (w_I)_{I=x,y,z} \in L^2(\Sigma_h)^3 : \mu_I|_{\mathcal{F}} \in M_h(\mathcal{F}), \forall \mathcal{F} \in \Sigma_h \}. $$

We employ polynomials of equal degree in (3.11), as was done in [35, 4, 21]; specifically, with $\mathcal{P}^k(D)$ denoting the space of polynomials of degree k defined on domain $D \subset \mathbb{R}^3$,

$$ U_h(K) = V_h(K) = \mathcal{P}^k(K), \quad M_h(\mathcal{F}) = \mathcal{P}^k(\mathcal{F}). \qquad (3.12) $$

**Approximate problem** To obtain the discrete problem, we integrate the local problem (3.7) against test functions $(\phi, \Psi) \in U_h \times V_h$, and the problems on the interfaces against test functions $\xi \in M_h$, and carry out integration by parts. In the weak form associated with (3.7a), we also employ the following definition for the numerical trace of the traction,

$$ \widehat{\sigma \nu} = \sigma_h \nu - \tau_u (u_h - \lambda_{uh}), \quad \text{with stabilization matrix } \tau_u = (\tau_{IJ})_{I,J=x,y,z}, \qquad (3.13) $$

and carry out an inverse integration by parts. For more details in obtaining (3.14), see [31].

The approximate problem reads as follows: Find $(u_h, \sigma_h, \lambda_{uh}) \in U_h \times V_h \times M_h$ that solves, 1. Local problems on element $K^e \in \mathcal{T}_h$, for all test functions $(\phi, \Psi) \in U_h \times V_h$,

$$ -\omega^2 \int_K \rho \mathbf{u}_h \cdot \bar{\phi} \, d\mathbf{x} - \int_K (\nabla \cdot \sigma_h) \cdot \bar{\phi} \, d\mathbf{x} + \int_{\partial K} \bar{\phi} \cdot \tau_u (\mathbf{u}_h - \lambda_{uh}) \, ds_{\mathbf{x}} = \int_K \mathbf{f} \cdot \bar{\phi} \, d\mathbf{x}, \quad (3.14a) $$

$$ \int_K \sigma_h : \mathbf{S} : \bar{\Psi} \, d\mathbf{x} = - \int_K \mathbf{u}_h \cdot \nabla \cdot \bar{\Psi} \, d\mathbf{x} + \int_{\partial K} \nu \cdot \bar{\Psi} \cdot \lambda_{uh} \, ds_{\mathbf{x}}. \qquad (3.14b) $$

2. Interface and boundary problems: the transmission and Neumann conditions give,

$$ \sum_{e=1}^{|\mathcal{T}_h|} \int_{\partial K^e \cap (\Sigma_{\text{int}} \cup \Sigma_{\text{N}})} (\nu \cdot \sigma_h \cdot \bar{\xi} - \bar{\xi} \cdot \tau_u \cdot (\mathbf{u}_h - \lambda_{uh})) \, ds_{\mathbf{x}} = 0, \quad \forall \xi \in M_h(\Sigma_{\text{int}} \cup \Sigma_{\text{N}}). \qquad (3.15) $$

The Robin boundary conditions (3.9) give, $\forall \xi \in M_h(\Sigma_\infty)$,

$$ \sum_{e=1}^{|\mathcal{T}_h|} \int_{\partial K^e \cap \Sigma_\infty} (\nu \cdot \sigma_h \cdot \bar{\xi} - \bar{\xi} \cdot \tau_u \cdot (\mathbf{u}_h - \lambda_{uh}) + i\omega \bar{\xi} \cdot \mathcal{Z} \cdot \lambda_{uh}) \, ds_{\mathbf{x}} = 0. \qquad (3.16) $$

The Dirichlet condition is imposed weakly, $\forall \xi \in M_h(\Sigma_D)$,

$$ \sum_{e=1}^{|\mathcal{T}_h|} \int_{\partial K^e \cap \Sigma_D} (\lambda_{uh} - \mathbf{g}^D) \cdot \bar{\xi} \, ds_{\mathbf{x}} = 0, \quad \text{on } \Sigma_D. \qquad (3.17) $$
## 3.3 Discretization of HDG problem

The discussion in this section employs the Voigt notation introduced in Subsection 2.2.

### 3.3.1 Discrete unknowns

**Local Basis functions** We have three groups of basis functions for the local finite element spaces introduced in (3.11),

$$ 
\begin{array}{l|l|l}
\text{Local finite element space} & \text{Basis functions} & \text{Dimension} \\
\hline
U_h(K^e), \quad 1 \le e \le |\mathcal{T}_h| & \phi_j^e, \quad 1 \le j \le n_e & \mathbf{n}_e \\
V_h(K^e), \quad 1 \le e \le |\mathcal{T}_h| & \psi_j^e, \quad 1 \le j \le m_e & \mathbf{m}_e \\
M_h(\mathbf{F}^k) = M_h(\mathbf{F}^{(e,\ell)}), \quad 1 \le k \le |\Sigma_h| & \xi_j^k = \xi_j^{(e,\ell)}, \quad 1 \le j \le \hat{n}_k & \hat{n}_k = \hat{n}_{(e,\ell)}.
\end{array}
\quad (3.18) $$

We also denote the total number of face degrees of freedom in each direction $I$ by,

$$ \hat{n}_I = \hat{n} = \sum_{k=1}^{|\Sigma_h|} \hat{n}_k, \quad I = x, y, z. \quad (3.19) $$

**Volume discrete unknowns** The displacement vector field $\mathbf{u}$ and strain tensor $\boldsymbol{\sigma}$ are approximated on each cell $K^e$, $1 \le e \le |\mathcal{T}_h|$ as follows. For the displacement we have,

$$ \mathbf{u}^h = \sum_{I=x,y,z} u_I^h \hat{\mathbf{e}}_I, \quad \text{with} \quad u_I^h|_{K^e} = \sum_{j=1}^{n_e} u_{Ij}^e \phi_j^e, \quad \text{for } I = x, y, z. \quad (3.20) $$

For the strain tensor, we work with its Voigt representation as a vector of length 6, cf. (2.16),

$$ \boldsymbol{\sigma}_h = \sum_{\mathcal{J} \in \mathcal{I}_{\text{sm}}} \sigma_{\mathcal{J}}^h \bar{\mathfrak{e}}_{\mathcal{J}} \quad \longleftrightarrow \quad \bar{\boldsymbol{\sigma}}_h = \sum_{\mathcal{J} \in \mathcal{I}_{\text{sm}}} \sigma_{\mathcal{J}}^h \hat{\mathfrak{e}}_{\mathcal{J}}, \quad \text{with } \sigma_{\mathcal{J}}^h|_{K^e} = \sum_{j=1}^{m_e} \mathbf{v}_{\mathcal{J}j}^e \psi_j^e. \quad (3.21) $$

We gather all volume unknowns into vector $\mathbf{W}$ such that,

$$ 
\begin{array}{l}
\mathbf{W} = (\mathbf{W}^e)_{1 \le e \le |\mathcal{T}_h|}, \quad \text{with} \quad \mathbf{W}^e := \begin{pmatrix} \mathbf{U}^e \\ \mathbf{V}^e \end{pmatrix} \quad \text{and} \quad \mathbf{U}^e = (\mathbf{U}_I^e)_{I=x,y,z}, \quad \mathbf{V}^e = (\mathbf{V}_\mathcal{J}^e)_{\mathcal{J} \in \mathcal{I}_{\text{sm}}}, \\
\text{with sub-blocks} \quad \mathbf{U}_I^e = (\mathbf{u}_{Ij}^e)_{j=1,...,n_e}, \quad \mathbf{V}_\mathcal{J}^e = (\mathbf{v}_{\mathcal{J}j}^e)_{j=1,...,m_e}.
\end{array}
\quad (3.22) $$

**Edge discrete unknowns** The unknowns defined on the edges $\Sigma_h$ are numerical approximations of the trace of the displacement $\mathbf{u}$ on the edges $\Sigma_h$. Specifically, on each edge $\mathbf{F}^k$, $k = 1, \dots, \hat{n}_k$,

$$ \boldsymbol{\lambda}_h = \sum_{I \in \{x,y,z\}} \lambda_I^h \hat{\mathbf{e}}_I, \quad \text{with} \quad \lambda_I^h|_{\mathcal{F}^k} = \sum_{j=1}^{\hat{n}_k} \lambda_{Ij}^k \xi_j^k. \quad (3.23) $$

We gather these coefficients into a global vector $\boldsymbol{\Lambda}$ which has the following substructure,

$$ \boldsymbol{\Lambda} = (\boldsymbol{\Lambda}^k)_{k=1,...,|\Sigma_h|}, \quad \text{with} \quad \boldsymbol{\Lambda}^k = (\boldsymbol{\Lambda}_I^k)_{I=x,y,z} \quad \text{and} \quad \boldsymbol{\Lambda}_I^k = (\lambda_{Ij}^k)_{1 \le j \le \hat{n}_k}. \quad (3.24) $$

**Remark 3.** The above expressions are written in terms of the global edge indices. They can also be written in terms of the local edge indices, specifically for $(e, \ell)$ identified with $k$, then

$$ \mathbf{F}^k = \mathbf{F}^{(e,\ell)}, \quad \lambda_{Ij}^{(e,\ell)} = \lambda_{Ij}^k, \quad \xi_j^{(e,\ell)} = \xi_j^k. \quad (3.25) $$
### 3.3.2 Discretized local problem

The discretized problem defined on each element $K^e$ is written as,

$$ \mathbb{A}^e \mathbf{W}^e + \mathbb{D}^e \mathbf{\Lambda}^e = \mathbf{S}^e, \quad \text{for } 1 \le e \le |\mathcal{T}_h|, \qquad (3.26) $$

with coefficient matrices and sources given by,

$$ \mathbb{A}^e = \begin{pmatrix} -\omega^2 \mathbb{M}_{\mathbf{u}}^e + \mathbb{M}_{\partial}^e & -\mathbb{K}_{\sigma}^e \\ -\mathbb{K}_{\mathbf{u}}^e & -\mathbb{M}_{\sigma}^e \end{pmatrix}, \quad \mathbb{D}^e = \begin{pmatrix} -\mathbb{D}_{\mathbb{m}}^{(e,1)} & \dots & -\mathbb{D}_{\mathbb{m}}^{(e,n_{\text{face}}^e)} \\ \mathbb{D}_{\mathbb{c}}^{(e,1)} & \dots & \mathbb{D}_{\mathbb{c}}^{(e,n_{\text{face}}^e)} \end{pmatrix}, \quad \mathbf{S}^e = \begin{pmatrix} \mathbf{S}_{\mathbb{m}}^e \\ \mathbf{0}_{(6\mathbb{m}_e)\times 1} \end{pmatrix}. \qquad (3.27) $$

We list below the components of the sub-blocks.

**Description of components** The non-zero entries of the source term come from integrating the volume source $\mathbf{f}$ in the equation of motion with test functions:

$$ \mathbf{S}_{\mathbb{m}}^{e} = (\mathbb{m} \mathbf{S}_{I}^{e})_{I=x,y,z}, \quad \text{with} \quad [\mathbb{m} \mathbf{S}_{I}^{e}]_{i} = \langle \mathbf{f} | K^{e}, \phi_{i}^{e} \hat{\mathbf{e}}_{I} \rangle, \quad 1 \le i \le n_{e}. \qquad (3.28) $$

The upper row blocks of $\mathbb{A}^e$ and $\mathbb{D}^e$ come from integrating the equation of motion (3.14a),

$$ \phi_i^e \hat{\mathbf{e}}_I, \quad I = x, y, z \quad 1 \le i \le n_e, \qquad (3.29) $$

while those of the second row-block come from<fnref n="1" /> integrating the constitutive equation (3.14b) against test functions,

$$ \Psi = \psi_i^e \bar{\mathfrak{e}}_\mathcal{J} \quad \Leftrightarrow \quad \overrightarrow{\Psi} = \psi_i^e \hat{\mathbf{e}}_\mathcal{J}, \quad \text{for } \mathcal{J} \in \mathcal{I}_{\text{sm}}, \quad 1 \le i \le m_e. \qquad (3.31) $$

For $1 \le i, j \le n_e$,

$$ \mathbb{M}_{\mathbf{u}}^{e} = (\mathbb{M}_{IJ}^{\mathbf{u}e})_{I,J=x,y,z}, \quad \text{with} \quad [\mathbb{M}_{II}^{\mathbf{u}e}]_{ij} = \begin{cases} \displaystyle \int_{K^e} \rho \overline{\phi_i^e} \phi_j^e \mathrm{d}\mathbf{x}, & I = J, \\ 0, & I \neq J \end{cases}, \qquad (3.32a) $$

$$ \mathbb{M}_{\partial}^{e} = (\mathbb{M}_{IJ}^{\partial e})_{I,J=x,y,z}, \quad \text{with} \quad [\mathbb{M}_{IJ}^{\partial e}]_{ij} = \sum_{\ell=1}^{n_{\text{face}}} \int_{\mathcal{F}^{(e,\ell)}} \tau_{IJ} \overline{\phi_i^e} \phi_j^e \mathrm{d}\mathbf{s}_x. \qquad (3.32b) $$

For $1 \le i, j \le m_e$,

$$ \mathbb{M}_{\sigma}^{e} = (\mathbb{M}_{\mathcal{J}\mathcal{J}'}^{\sigma e})_{\mathcal{J},\mathcal{J}' \in \mathcal{I}_{\text{sm}}} \quad \text{with} \quad [\mathbb{M}_{\mathcal{J}\mathcal{J}'}^{\sigma e}]_{ij} = \hat{\mathbf{e}}_{\mathcal{J}} \cdot \left( \int_{K^e} \overline{\psi_i^e} \psi_j^e \overline{\hat{\mathbf{S}}^\dagger} \mathrm{d}\mathbf{x} \right) \cdot \hat{\mathbf{e}}_{\mathcal{J}'} \qquad (3.33) $$

For $1 \le i \le n_e, 1 \le j \le m_e$,

$$ \mathbb{K}_{\sigma}^{e} = (\mathbb{K}_{I\mathcal{J}}^{\sigma e})_{I=x,y,z}, \quad \text{with} \quad [\mathbb{K}_{I\mathcal{J}}^{\sigma e}]_{ij} = \hat{\mathbf{e}}_I \cdot \left( \int_{K^e} \overline{\phi_i^e} \mathbb{A}^{\dagger}(\nabla \psi_j^e) \mathrm{d}\mathbf{x} \right) \cdot \hat{\mathbf{e}}_{\mathcal{J}}, \qquad (3.34a) $$

$$ \mathbb{K}_{\mathbf{u}}^{e} = (\mathbb{K}_{\mathcal{J}J}^{\mathbf{u}e})_{\mathcal{J} \in \mathcal{I}_{\text{sm}}, J=x,y,z} \quad \text{with} \quad [\mathbb{K}_{\mathcal{J}J}^{\mathbf{u}e}]_{ij} = \hat{\mathbf{e}}_J \cdot \left( \int_{K^e} \phi_j^e \overline{\mathbb{A}^{\dagger}(\nabla \psi_i^e)} \mathrm{d}\mathbf{x} \right) \cdot \hat{\mathbf{e}}_{\mathcal{J}}. \qquad (3.34b) $$

$$ (\nabla \cdot \boldsymbol{\sigma}) \cdot \mathbf{e}_I = \sum_{\mathcal{J} \in \mathcal{I}_{\text{sm}}} \mathbf{e}_I \cdot \mathbb{A}^{\dagger}(\nabla \sigma_{\mathcal{J}}) \cdot \hat{\mathbf{e}}_{\mathcal{J}}, \quad (\mathbf{S}\boldsymbol{\sigma}) : \mathbf{w} = \boldsymbol{\sigma} \cdot \mathbf{S} \cdot \mathbf{w} = \overrightarrow{\mathbf{w}} \cdot \hat{\mathbf{S}}^{\dagger} \cdot \overline{\hat{\mathbf{S}}^{\dagger}} \cdot \overline{\boldsymbol{\sigma}}. \qquad (3.30) $$

<footnote n="1">¹To describe the components of the matrices in Voigt quantities, the following identities are employed: for $\sigma \in S_2$ and vector $\mathbf{w}$,</footnote>
For $1 \le i \le m_e, 1 \le j \le \hat{n}_{(e,l)}$

$$ \mathbb{D}_{\mathrm{m}}^{(e,l)} = \left(_{\mathrm{m}}\mathbb{D}_{IJ}^{(e,l)}\right)_{I,J=x,y,z} \quad \text{with} \quad \left[_{\mathrm{m}}\mathbb{D}_{IJ}^{(e,l)}\right]_{ij} = \int_{\mathbb{F}^{(e,l)}} \tau_{IJ} \bar{\phi}_i^e \xi_j^{(e,l)} \mathrm{d}\mathbf{s}_x. \qquad (3.35) $$

For $1 \le i \le n_e, 1 \le j \le \hat{n}_{(e,l)}$,

$$ \mathbb{D}_{\mathrm{c}}^{(e,l)} = \left(_{\mathrm{c}}\mathbb{D}_{\mathfrak{J}J}^{(e,l)}\right)_{\mathfrak{J} \in \mathcal{I}_{\mathrm{sm}}, J=x,y,z} \quad \text{with} \quad \left[_{\mathrm{c}}\mathbb{D}_{\mathfrak{J}J}^{(e,l)}\right]_{ij} = \int_{\mathbb{F}^{(e,l)}} \overline{\psi_i^e} \xi_j^{(e,l)} \mathbb{A}^{\dagger}(\boldsymbol{\nu}^{(e,l)})_{\mathfrak{J}\mathfrak{J}'} \mathrm{d}\mathbf{s}_x. \qquad (3.36) $$

### 3.3.3 Discretized problem defined on edges

The discretization of (3.15)-(3.17) defined on $\Sigma$ takes the following form,

$$ \sum_{e=1}^{|\mathcal{T}|} \mathcal{R}_e^t (\mathbb{B}^e \mathbf{W}^e + \mathbb{L}^e \mathcal{R}_e \mathbf{\Lambda}) = \sum_{e=1}^{|\mathcal{T}|} \mathcal{R}_e^t \mathbf{s}^e, \qquad (3.37) $$

with coefficient matrices and sources,

$$ \mathbb{B}^e = \begin{pmatrix} -\mathbb{B}_{\mathbf{u}}^{(e,1)} & \mathbb{B}_{\boldsymbol{\sigma}}^{(e,1)} \\ \vdots & \vdots \\ -\mathbb{B}_{\mathbf{u}}^{(e,n_{\mathrm{face}}^e)} & \mathbb{B}_{\boldsymbol{\sigma}}^{(e,n_{\mathrm{face}}^e)} \end{pmatrix}, \quad \mathbb{L}^e = \begin{pmatrix} \mathbb{L}^{(e,1)} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbb{L}^{(e,2)} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbb{L}^{(e,n_{\mathrm{face}}^e)} \end{pmatrix}, \quad \mathbf{s}^e = \begin{pmatrix} \mathfrak{s}^{(e,1)} \\ \vdots \\ \mathfrak{s}^{(e,n_{\mathrm{face}}^e)} \end{pmatrix}, \qquad (3.38) $$

and the local-to-global map $\mathcal{R}_e$ which gives the restriction of a global vector to an element $K^e$,

$$ \mathcal{R}_e = \begin{pmatrix} [\mathcal{R}_e]_{(e,1),1} & [\mathcal{R}_e]_{(e,1),2} & \cdots & [\mathcal{R}_e]_{(e,1),|\Sigma_h|} \\ [\mathcal{R}_e]_{(e,2),1} & [\mathcal{R}_e]_{(e,2),2} & \cdots & [\mathcal{R}_e]_{(e,2),|\Sigma_h|} \\ \vdots & \vdots & \ddots & \vdots \\ [\mathcal{R}_e]_{(e,n_{\mathrm{face}}^e),1} & [\mathcal{R}_e]_{(e,n_{\mathrm{face}}^e),2} & \cdots & [\mathcal{R}_e]_{(e,n_{\mathrm{face}}^e),|\Sigma_h|} \end{pmatrix}. \qquad (3.39) $$

The components of their sub-blocks are listed below.

**Descriptions of components** The blocks of the local-to-global map operator are given by, for $1 \le \ell \le n_{\mathrm{face}}^e$, and $k$ with $1 \le k \le |\Sigma|$,

$$ [\mathcal{R}_e]_{(e,\ell),k} = \begin{cases} \mathbb{L}_{3\hat{n}_k \times 3\hat{n}_k}, & \text{if } \mathbb{F}^k \cap \partial K^e \neq \emptyset \text{ with } \mathbb{F}^k = \mathbb{F}^{(e,\ell)}, \\ \mathbf{0}_{3\hat{n}_{(e,\ell)} \times 3\hat{n}_k}, & \text{if } \mathbb{F}^k \cap \partial K^e = \emptyset. \end{cases} \qquad (3.40) $$

The vector $\mathbf{s}^e$, and matrices $\mathbb{B}^e$ and $\mathbb{L}^e$ have $n_{\mathrm{face}}^e$ row-blocks, labeled by $(e, \ell)$ and corresponding to faces of $K^e$. They have the following block structures:

$$ \mathbb{B}_{\mathbf{u}}^{(e,\ell)} = \left(_{\mathbf{u}}\mathbb{B}_{IJ}^{(e,\ell)}\right)_{I,J=x,y,z}, \quad \mathbb{B}_{\boldsymbol{\sigma}}^{(e,\ell)} = \left(_{\boldsymbol{\sigma}}\mathbb{B}_{\mathfrak{J}\mathfrak{J}'}^{(e,\ell)}\right)_{I=x,y,z, \mathfrak{J}' \in \mathcal{I}_{\mathrm{sm}}}, \\ \mathbb{L}^{(e,\ell)} = \left(\mathbb{L}_{IJ}^{(e,\ell)}\right)_{I,J=x,y,z}, \quad \mathfrak{s}^{(e,\ell)} = (\mathfrak{s}_I^{(e,\ell)})_{I=x,y,z}. $$
The j-th row of block-row (e, l) is obtained<fnref n="2" /> by integrating corresponding boundary conditions (3.15) and (3.16) with respect to test functions (3.42),

$$ \xi|_{F^{(e,l)}} = \xi_j^{(e,l)} \hat{e}_I, \quad \text{with } I = x, y, z. \qquad (3.42) $$

For $1 \le i \le \hat{n}_k$, $1 \le j \le n_e$,

$$ \left[ \mathbf{u} \mathbb{B}_{IJ}^{(e,l)} \right]_{ij} = \begin{cases} \int_{\mathbb{F}^{(e,l)}} \tau_{IJ} \overline{\xi_i^{(e,l)}} \phi_j^e \mathrm{d}\mathbf{s}_x, & \mathbb{F}^{(e,l)} \in \Sigma_{\text{int}} \cup \Sigma_{\text{N}} \cup \Sigma_{\infty}, \\ 0, & \mathbb{F}^{(e,l)} \in \Sigma_{\text{D}}. \end{cases} \qquad (3.43) $$

For $1 \le i, j \le \hat{n}_k$,

$$ \left[ \mathbb{L}_{IJ}^{(e,l)} \right]_{ij} = \begin{cases} \int_{\mathbb{F}^{(e,l)}} \tau_{IJ} \overline{\xi_i^{(e,l)}} \xi_j^{(e,l)} \mathrm{d}\mathbf{s}_x, & \mathbb{F}^{(e,l)} \in \Sigma_{\text{int}} \cup \Sigma_{\text{N}}, \\ \int_{\mathbb{F}^{(e,l)}} \overline{\xi_i^{(e,l)}} \xi_j^{(e,l)} (\tau_{IJ} + \mathcal{Z}_{IJ}^{\text{abc}}) \mathrm{d}\mathbf{s}_x, & \mathbb{F}^{(e,l)} \in \Sigma_{\infty}, \\ \int_{\mathbb{F}^{(e,l)}} \overline{\xi_i^{(e,l)}} \xi_j^{(e,l)} \mathrm{d}\mathbf{s}_x, & \mathbb{F}^{(e,l)} \in \Sigma_{\text{D}}. \end{cases} \qquad (3.44) $$

For $1 \le i \le \hat{n}_k$, $1 \le j \le m_e$,

$$ \left[ \boldsymbol{\sigma} \mathbb{B}_{I\bar{J}}^{(e,l)} \right]_{ij} = \begin{cases} \int_{\mathbb{F}^{(e,l)}} \overline{\xi_i^{(e,l)}} \psi_j^e \mathcal{A}^{\dagger}(\boldsymbol{\nu}^{(e,l)})_{I\bar{J}} \mathrm{d}\mathbf{s}_x, & \mathbb{F}^{(e,l)} \in \Sigma_{\text{int}} \cup \Sigma_{\text{N}} \cup \Sigma_{\infty}, \\ 0, & \mathbb{F}^{(e,l)} \in \Sigma_{\text{D}}. \end{cases} \qquad (3.45) $$

We recall the Dirichlet boundary vector source $\mathbf{g}^D = (\mathbf{g}_I^D)_{I=x,y,z}$. The components of the row block for the source $\mathbf{s}^e$ are, for $1 \le i \le \hat{n}_k$,

$$ [\mathbf{s}_I^{(e,l)}]_i = \begin{cases} 0, & \mathbb{F}^{(e,l)} \in \Sigma_{\text{int}} \cup \Sigma_{\text{N}} \cup \Sigma_{\infty}, \\ \int_{\mathbb{F}^{(e,l)}} \overline{\mathbf{g}_I^D} \xi_j^{(e,l)} \mathrm{d}\mathbf{s}_x, & \mathbb{F}^{(e,l)} \in \Sigma_{\text{D}}. \end{cases} \qquad (3.46) $$

### 3.3.4 Summary of discrete problem in HDG method

The discretization of (3.14)–(3.17) takes the following form with discrete unknowns $(\mathbf{W}, \mathbf{\Lambda})$,

$$ \mathbb{A}^e \mathbf{W}^e + \mathbb{D}^e \mathcal{R}_e \mathbf{\Lambda} = \mathbf{S}^e, \quad \forall e = 1, \dots, |\mathcal{T}_h|, \qquad (3.47a) $$

$$ \sum_{e=1}^{|\mathcal{T}|} \mathcal{R}_e^t (\mathbb{B}^e \mathbf{W}^e + \mathbb{L}^e \mathcal{R}_e \mathbf{\Lambda}) = \sum_{e=1}^{|\mathcal{T}|} \mathcal{R}_e^t \mathbf{s}^e. \qquad (3.47b) $$

Problem (3.47) can be reduced to one in terms of $\mathbf{\Lambda}$ only, called the *global problem*,

$$ \mathbb{K} \mathbf{\Lambda} = \mathcal{S}, \quad \text{with} \quad \mathcal{S} := \sum_{e=1}^{|\mathcal{T}|} \mathcal{R}_e^t (\mathbf{s}^e - \mathbb{B}_e \mathbb{A}_e^{-1} \mathbf{S}^e), \qquad (3.48a) $$

$$ \text{and } \mathbb{K} := \sum_{e=1}^{|\mathcal{T}_h|} \mathcal{R}_e^t \mathbb{K}^e \mathcal{R}_e, \quad \text{where } \mathbb{K}^e := -\mathbb{B}_e (\mathbb{A}^e)^{-1} \mathbb{D}^e + \mathbb{L}_e. \qquad (3.48b) $$

$$ \boldsymbol{\sigma}\boldsymbol{\nu} = \mathbb{A}^{\dagger}(\boldsymbol{\nu})\boldsymbol{\sigma}; \quad \boldsymbol{\sigma}\boldsymbol{\nu} \cdot \mathbf{w} = \mathbf{w} \cdot \boldsymbol{\sigma} \cdot \boldsymbol{\nu} = \mathbf{w} \cdot \mathbb{A}^{\dagger}(\boldsymbol{\nu}) \cdot \boldsymbol{\sigma}. \qquad (3.41) $$

<footnote n="2">²We have employed the identity: for two vectors $\boldsymbol{\nu}$, $\mathbf{w}$ and a symmetric matrix $\boldsymbol{\sigma} \in S_2$,</footnote>
This means that in the HDG method the problem is solved in two stages. Firstly, one solves the global problem (3.49a) in terms of the trace $\Lambda$; secondly, the values of the volume unknowns $\mathbf{W}^e$ are retrieved element-by-element with the right-hand side of (3.49b),

$$ (3.47) \quad \Leftrightarrow \quad \left\{ \begin{array}{l} \mathbb{K} \Lambda = S, \hfill (3.49\text{a}) \\ \mathbf{W}^e = (\mathbb{A}^e)^{-1} (-\mathbb{D}^e \mathcal{R}_e \Lambda + \mathbf{S}^e), \quad \forall e = 1, \dots, |\mathcal{T}_h|. \hfill (3.49\text{b}) \end{array} \right. $$

## 3.4 Stabilization matrix $\tau_u$

The stabilization matrix $\tau_u$ appears in the definition of the numerical trace of the traction (3.13) and its subsequent discretization, e.g., (3.43) and (3.44). A judicial choice of $\tau_u$ is necessary to provide accurate results. In its most general form, $\tau_u$ is a symmetric positive definite matrix. The most common choice is an identity-based stabilization, in which $\tau_u$ is a scalar multiple of the identity matrix. Also proposed in the literature for elastostatics, cf., [35, 22, 15] is the Kelvin–Christoffel (KC) stabilization which corresponds to a scalar scaling multiple of the 3 × 3 Kelvin–Christoffel matrix $\Gamma$ (2.19).

In addition to these two families of stabilization, we construct in Section 4 the hybridized Godunov stabilization defined in terms of the Godunov matrix $M_{\text{Godunov}}$. Its definition for a mesh cell $K$ with outward-pointing unit normal vector $\nu$ is,

$$ M_{\text{Godunov}}(\nu) = \begin{cases} M_{\text{Giso}}(\nu), & \text{isotropy;} \\ M_{\text{Gani}}(\nu), & \text{anisotropy with 3 distinct speeds } \rho c_{\text{qP}}^2 > \rho c_{\text{qS1}}^2 > \rho c_{\text{qS2}}^2, \end{cases} \quad (3.50) $$

$$ \text{with} \quad M_{\text{Giso}}(\nu) = \frac{\rho}{c_P + c_S} \left( c_P c_S \mathbb{I}_d + \frac{\Gamma_{\text{iso}}(\nu)}{\rho} \right) = \rho (c_S \mathbb{I}_d + (c_P - c_S) \nu \otimes \nu), \quad (3.51) $$

$$ \text{and } M_{\text{Gani}}(\nu) := \rho (c_{\text{qS1}} + c_{\text{qS2}} + c_{\text{qP}}) \left( \mathbb{I}_d + \gamma \left( \frac{\Gamma(\nu)}{\rho} + \mathfrak{p}_2 \mathbb{I}_d \right)^{-1} \right), \quad (3.52) $$

$$ \text{where } \mathfrak{p}_2 := c_{\text{qS1}} c_{\text{qS2}} + c_{\text{qS1}} c_{\text{qP}} + c_{\text{qS2}} c_{\text{qP}}, \quad \text{and} \quad \gamma := \frac{c_{\text{qS1}} c_{\text{qS2}} c_{\text{qP}}}{c_{\text{qS1}} + c_{\text{qS2}} + c_{\text{qP}}} - \mathfrak{p}_2. $$

Here, the density $\rho$ and wavespeeds are associated with mesh cell $K$. Hybridization for isotropy was carried out in [36], we also refer to [31] which shows the derivation of $M_{\text{Giso}}$ in terms of $\Gamma$. The derivation for $M_{\text{Gani}}$ is given in Section 4. In Section 5, we will investigate the three families (identity-based, KC-based and Godunov-based) and show the optimality of the Godunov stabilization.

## 4 Construction of hybridized Godunov stabilization operators

In this section, we extend the hybridization method in [36] for isotropic elasticity to anisotropy, and construct the hybridized Godunov stabilization operator $M_{\text{Godunov}}(\nu)$ (3.50). We recall that a stabilization matrix is employed to define the numerical trace of traction $\hat{\sigma}\nu$ on each mesh element. We consider an elastic material whose KC matrix $\Gamma(\nu)$ (2.19) has three distinct eigenvalues,

$$ \rho c_{\text{qP}}^2 > \rho c_{\text{qS1}}^2 > \rho c_{\text{qS2}}^2. \quad (4.1) $$

The discussion in this section employs the Voigt notation introduced in Subsection 2.2– Subsection 2.4.
In Subsection 4.1, we introduce a first-order formulation of the time-dependent elastic equation (4.3) with unknown $\mathbf{q} = (\epsilon, \rho\mathbf{v})$, and make appear the flux term $\mathfrak{B}(\nu)\mathbf{q}$ (4.8) whose last three components give $\sigma\nu$. This means that the numerical trace of traction $\hat{\sigma}\nu$ at a face F can be obtained from the numerical trace $\hat{\mathfrak{B}}(\nu)\mathbf{q}$. Additionally, we consider $\mathbf{q}$ as a collection of values on each mesh element, i.e, $\mathbf{q} = \cup_{K \in \mathcal{T}_h} \mathbf{q}^K$. At an interface F shared by two elements $K^+$ and $K^-$, we denote by $\mathbf{q}^{\pm}$ the trace of $\mathbf{q}$ from each side of F associated with the neighbor cells $K^{\pm}$, i.e.

$$F = \partial K^- \cap \partial K^+, \quad F^{\pm} := \partial K^{\pm} \cap F. \qquad (4.2)$$

For the construction of $\hat{\mathfrak{B}}(\nu)\mathbf{q}$, we work with the exact solutions of the Riemann problem in a neighborhood of F, having as initial data $\mathbf{q}^{\pm}$, denoted by $\mathrm{RP}(\mathbf{q}^{-}, \mathbf{q}^{+})$. The solution of $\mathrm{RP}(\mathbf{q}^{-}, \mathbf{q}^{+})$ to the immediate left and right of F is called intermediate states $\mathbf{q}^{*}$.

The main construction is given in Subsection 4.2 and consists of 3 main steps.

- The numerical flux $\hat{\mathfrak{B}}(\nu)\mathbf{q}$ is given by the Godunov flux, which by definition is the value of $\mathfrak{B}(\nu)$ at the intermediate state $\mathbf{q}^*$, i.e. $\mathfrak{B}(\nu)\mathbf{q}^*$.

- The relations between the flux of the intermediate states $\mathbf{q}^*$ and the left and right data $\mathbf{q}^{\pm}$ are given by a system of Rankine-Hugoniot (RH) jump conditions, cf. (4.18). From this, we derive the usual transmission conditions for elasticity (4.19), and an expression of the Godunov flux, cf. (4.37), and thus of the traction, cf. (4.40), which is in terms of one-sided data (i.e., either $\mathbf{q}^{+}$ or $\mathbf{q}^{-}$).

- Finally, in a step called 'hybridization', these above ingredients are combined to obtain the hybridized HDG numerical flux (4.42).

For the hybridization in isotropy, we refer to [36] and [31]; in the latter, the derivation is written with Voigt notation which brings out the connection between the stabilization operator and the Kelvin-Christoffel matrix. Details for the identities employed in this section with Voigt notation as well as a review of the Riemann problem can be found in [31].

**Remark 4.** We adopt the name 'Godunov-upwind' from [5, 6, 7], in which the hybridization was carried out following different approaches; however, employing Rankine-Hugoniot jump condition is recognized in sequel work [6, 7] (where the hybridized flux is also called 'Rankine-Hugoniot' flux) to be more natural and direct.

## 4.1 First order system and flux operator

We rewrite the time-dependent elastic equation as a first order system with unknowns $(\epsilon, \rho\mathbf{v})$, in recalling $\sigma = C\epsilon$, $\epsilon = \nabla^s\mathbf{u}$, and $\mathbf{v} = \partial_t\mathbf{u}$:

$$\left\{ \begin{array}{l} \partial_t \epsilon - \nabla^s \mathbf{v} = 0, \\ \partial_t \rho \mathbf{v} - \nabla \cdot (C \epsilon) = 0 \end{array} \right. \Leftrightarrow \partial_t \begin{pmatrix} \epsilon \\ \rho \mathbf{v} \end{pmatrix} = \begin{pmatrix} 0 & \nabla^s \rho^{-1} \\ \nabla \cdot C & 0 \end{pmatrix} \begin{pmatrix} \epsilon \\ \rho \mathbf{v} \end{pmatrix}. \qquad (4.3)$$

With matrix $A_I$ (2.22), define matrix-valued differential operator $A(\partial_x)$ and its transpose,

$$A(\partial_x) := \sum_{I=x,y,z} \partial_I A_I, \quad A(\partial_x)^t = \sum_I \partial_I A_I^t. \qquad (4.4)$$
Employing the following identity, cf. [31],

$$ \overrightarrow{\nabla^g \mathbf{v}} = \mathfrak{A}(\partial_x)^t \mathbf{v}, \quad \nabla \cdot (\mathbf{C} \epsilon) = \mathfrak{A}(\partial_x)^{\dagger} \overline{\mathbf{C}^{\dagger} \epsilon}, \qquad (4.5) $$

equation (4.3) can be written as

$$ \partial_t \mathbf{q} + \mathfrak{B}(\partial_x) \mathbf{q} = \begin{pmatrix} 0 \\ \mathbf{f} \end{pmatrix}, \qquad (4.6) $$

$$ \text{with } \mathfrak{B}(\partial_x) := - \begin{pmatrix} 0 & \mathfrak{A}^t(\partial_x) \rho^{-1} \\ \mathfrak{A}(\partial_x)^{\dagger} \overline{\mathbf{C}^{\dagger}} & 0 \end{pmatrix} \quad \text{and} \quad \mathbf{q} := \begin{pmatrix} \overrightarrow{\epsilon} \\ \rho \mathbf{v} \end{pmatrix}. \qquad (4.7) $$

**Flux term** Along $\partial K$ which has outward pointing normal vector $\nu$, we define the matrix

$$ \mathfrak{B}(\nu) := - \begin{pmatrix} 0 & \mathfrak{A}^t(\nu)\rho^{-1} \\ \mathfrak{A}(\nu)^{\dagger} \overline{\mathbf{C}^{\dagger}} & 0 \end{pmatrix}. \qquad (4.8) $$

This is called the flux matrix at $\partial K$ and arises from integrating both sides of (4.3) against test function $(\chi/\mathbf{w})$ with $\chi \in S_2$ and vector $\mathbf{w} \in \mathbb{C}^3$, and carrying out an integration by parts,

$$ \int_K \partial_t \mathbf{q} \cdot \left( \begin{pmatrix} \overrightarrow{\chi} \\ \mathbf{w} \end{pmatrix} \right) d\mathbf{x} + \int_K \mathfrak{B}(\partial_x) \left( \begin{pmatrix} \overrightarrow{\epsilon} \\ \rho \mathbf{v} \end{pmatrix} \right) \cdot \left( \begin{pmatrix} \overrightarrow{\chi} \\ \mathbf{w} \end{pmatrix} \right) d\mathbf{x} \\ = \int_K \partial_t \mathbf{q} \cdot \left( \begin{pmatrix} \overrightarrow{\chi} \\ \mathbf{w} \end{pmatrix} \right) d\mathbf{x} + \int_K \left( \begin{pmatrix} \overrightarrow{\epsilon} \\ \rho \mathbf{v} \end{pmatrix} \right) \cdot \mathfrak{B}^{\text{ad}}(\partial_x) \left( \begin{pmatrix} \overrightarrow{\chi} \\ \mathbf{w} \end{pmatrix} \right) d\mathbf{x} + \int_{\partial K} \mathfrak{B}(\nu) \left( \begin{pmatrix} \overrightarrow{\epsilon} \\ \rho \mathbf{v} \end{pmatrix} \right) \cdot \left( \begin{pmatrix} \overrightarrow{\chi} \\ \mathbf{w} \end{pmatrix} \right) d\mathbf{s}_x. \qquad (4.9) $$

Here 'adjoint' operator associated with $\mathfrak{B}$ is

$$ \mathfrak{B}^{\text{ad}}(\partial_x) := \begin{pmatrix} 0 & \overline{\mathfrak{C}^{\dagger}} \mathfrak{A}(\partial_x)^t \\ \rho^{-1} \mathfrak{A}(\partial_x) & 0 \end{pmatrix}. \qquad (4.10) $$

**Eigenvalues of $\mathfrak{B}(\nu)$** Consider the eigenvalue problem of the flux operator $\mathfrak{B}(\nu)$,

$$ \begin{pmatrix} 0 & \mathfrak{A}^t(\nu)\rho^{-1} \\ \mathfrak{A}(\nu)^{\dagger} \overline{\mathbf{C}^{\dagger}} & 0 \end{pmatrix} \begin{pmatrix} \overrightarrow{\tau} \\ \mathbf{w} \end{pmatrix} = \alpha \begin{pmatrix} \overrightarrow{\tau} \\ \mathbf{w} \end{pmatrix} \quad \Leftrightarrow \quad \begin{cases} \mathfrak{A}^t(\nu)\rho^{-1}\mathbf{w} = \alpha \overrightarrow{\tau} \\ \mathfrak{A}(\nu)^{\dagger} \overline{\mathbf{C}^{\dagger} \tau} = \alpha \mathbf{w} \end{cases}. \qquad (4.11) $$

Applying $\mathfrak{A}(\nu)^{\dagger} \overline{\mathbf{C}^{\dagger}}$ to both sides of the first equation, we obtain the eigenproblem for the Kelvin-Christoffel matrix $\Gamma(\nu)$, which is also called the Christoffel equation,

$$ \mathfrak{A}(\nu)^{\dagger} \overline{\mathbf{C}^{\dagger}} \mathfrak{A}^t(\nu) \mathbf{w} = \rho \alpha^2 \mathbf{w} \quad \Leftrightarrow \quad \Gamma(\nu) \mathbf{w} = \rho \alpha^2 \mathbf{w}. \qquad (4.12) $$

We work under the assumption that $\Gamma(\nu)$ has the three distinct eigenvalues,

$$ \rho c_{qP}^2 < \rho c_{qS1}^2 < \rho c_{qS2}^2. \qquad (4.13) $$

This implies that the eigenvalues of $\mathfrak{B}(\nu)$ are

$$ -c_{qP}, -c_{qS1}, -c_{qS2}, 0, c_{qS2}, c_{qS1}, c_{qP}. \qquad (4.14) $$
multiplicity 3

**Remark 5.** For elastic isotropy, (2.24), the eigenvalues of $\mathfrak{B}(\nu)$ are

$$ -c_P, -c_S, 0, c_S, c_P. \qquad (4.15) $$
multiplicity 2
multiplicity 3
multiplicity 2
## 4.2 Derivation for anisotropic elasticity with distinct waves speeds

Along an interface F shared by mesh cells $K^+$ and $K^-$, cf. (4.2), we denote by $\nu^\pm$ the outward-pointing normal vectors of $K^\pm$, and by $\nu$ the normal vector pointing from $K^-$ to $K^+$, i.e. $\nu = \nu^- = -\nu^+$. To distinguish the background parameters on each side of F, we write

$$ \rho^{\pm}, C^{\pm}, \Gamma(\nu)^{\pm}, \mathfrak{B}(\nu)^{\pm} \text{ with eigenvalues } c_{\alpha}^{\pm}. \qquad (4.16) $$

Here, $\mathfrak{B}(\nu)^-$ is defined in (4.8) with the physical parameters of $K^-$, i.e., $(\rho^-, \lambda^-, \mu^-)$ and its normal vector $\nu^-$; similarly for $\mathfrak{B}(\nu)^+$. The KC matrix associated with each mesh cell, $\Gamma(\nu)^\pm$, is defined in the same manner with definition (2.19).

<figure><img src="image_2.png" /><figcaption>Figure 1: The eigenvalues, spectral structure and the Godunov states appearing in the Rankine-Hugoniot jump condition (4.18) for anisotropic elasticity with 3 distinct speeds.</figcaption></figure>

Under the assumption of spectral structure (4.14) of $\mathfrak{B}(\nu)^\pm$, we have 7 possible discontinuities propagating at the speed of corresponding eigenvalues indicated in (4.14). Denote the states in between the discontinuities by, cf. Figure 1,

$$ \mathbf{q}^{-}, \mathbf{q}^{a-}, \mathbf{q}^{b-}, \mathbf{q}^{*-}, \mathbf{q}^{*+}, \mathbf{q}^{b+}, \mathbf{q}^{a+}, \mathbf{q}^{+}. \qquad (4.17) $$

They satisfy jump condition, cf. [38, Equation (19)],

$$ \mathfrak{B}(\nu)^{-} (\mathbf{q}^{a-} - \mathbf{q}^{-}) = -c_{qP}^{-} (\mathbf{q}^{a-} - \mathbf{q}^{-}); \qquad (4.18a) $$

$$ \mathfrak{B}(\nu)^{-} (\mathbf{q}^{b-} - \mathbf{q}^{a-}) = -c_{qS1}^{-} (\mathbf{q}^{b-} - \mathbf{q}^{a-}); \qquad (4.18b) $$

$$ \mathfrak{B}(\nu)^{-} (\mathbf{q}^{*-} - \mathbf{q}^{b-}) = -c_{qS2}^{-} (\mathbf{q}^{*-} - \mathbf{q}^{b-}); \qquad (4.18c) $$

$$ \mathfrak{B}(\nu)^{-} \mathbf{q}^{*-} = \mathfrak{B}(\nu)^{+} \mathbf{q}^{*+}; \qquad (4.18d) $$

$$ \mathfrak{B}(\nu)^{+} (\mathbf{q}^{b+} - \mathbf{q}^{*+}) = c_{qS1}^{+} (\mathbf{q}^{b+} - \mathbf{q}^{*+}); \qquad (4.18e) $$

$$ \mathfrak{B}(\nu)^{+} (\mathbf{q}^{b+} - \mathbf{q}^{a+}) = c_{qS1}^{+} (\mathbf{q}^{b+} - \mathbf{q}^{a+}); \qquad (4.18f) $$

$$ \mathfrak{B}(\nu)^{+} (\mathbf{q}^{a+} - \mathbf{q}^{+}) = c_{qS}^{+} (\mathbf{q}^{a+} - \mathbf{q}^{+}). \qquad (4.18g) $$

**Transmission conditions for the interface states** We obtain the transmission conditions from the first 3 and last 3 components of equation (4.18d) associated with the non-propagative state $\mathbf{q}^{*\pm}$,

$$ \left\{ \begin{array}{l} \nu_I^- v_I^{*-} = \nu_I^+ v_I^{*+}, \quad I = x, y, z \\ (\sigma \nu)^{*} = (\sigma \nu)^{*+} \end{array} \right. \quad \Rightarrow \quad \left\{ \begin{array}{l} \mathbf{v}^{*-} = \mathbf{v}^{*+}; \\ (\sigma \mathbf{v})^{*-} = (\sigma \mathbf{v})^{*+}. \end{array} \right. \qquad (4.19) $$
Here $\nu_I^{\pm}$ are the components of the normal vectors $\nu^{\pm}$, i.e., $\nu^{\pm} = (\nu_x^{\pm}, \nu_y^{\pm}, \nu_z^{\pm})^t$.

**Relation between the intermediate states with the interface and boundary states**

Working on each side of the interface F, we will express the state $\mathbf{q}^{a\pm}$ and $\mathbf{q}^{b\pm}$ in terms of the interface states $\mathbf{q}^{*\pm}$ and the boundary states $\mathbf{q}^{\pm}$. To alleviate the notation in the derivation, we write

$$ \mathfrak{B}_{\pm} := \mathfrak{B}(\boldsymbol{\nu})^{\pm}. \qquad (4.20) $$

**Step 1a** We start by working with the first three equations of (4.18). Summing them yields,

$$ \mathfrak{B}_{-} (\mathbf{q}^{*-} - \mathbf{q}^{-}) = (-c_{qP}^{-} + c_{qS1}^{-}) \mathbf{q}^{a-} + (-c_{qS1}^{-} + c_{qS2}^{-}) \mathbf{q}^{b-} - c_{qS2}^{-} \mathbf{q}^{*-} + c_{qP}^{-} \mathbf{q}^{-}. \qquad (4.21) $$

After rearrangement, this leads to,

$$ (c_{qP}^{-} - c_{qS1}^{-}) \mathbf{q}^{a-} + (c_{qS1}^{-} - c_{qS2}^{-}) \mathbf{q}^{b-} = -\mathfrak{B}_{-} (\mathbf{q}^{*-} - \mathbf{q}^{-}) - c_{qS2}^{-} \mathbf{q}^{*-} + c_{qP}^{-} \mathbf{q}^{-}. \qquad (4.22) $$

Next, taking $c_{qP}^{-} \times (4.18a) + c_{qS1}^{-} \times (4.18b) + c_{qS2}^{-} \times (4.18c)$, we obtain,

$$ \mathfrak{B}_{-} \left[ \begin{aligned} & (c_{qP}^{-} - c_{qS1}^{-}) \mathbf{q}^{a-} + (c_{qS1}^{-} - c_{qS2}^{-}) \mathbf{q}^{b-} \\ & + [(c_{qS1}^{-})^2 - (c_{qP}^{-})^2] \mathbf{q}^{a-} + [(c_{qS2}^{-})^2 - (c_{qS1}^{-})^2] \mathbf{q}^{b-} - (c_{qS2}^{-})^2 \mathbf{q}^{*-} + (c_{qP}^{-})^2 \mathbf{q}^{-} \end{aligned} \right] = c_{qP}^{-} \mathfrak{B}_{-} \mathbf{q}^{-} - c_{qS2}^{-} \mathfrak{B}_{-} \mathbf{q}^{*-} \qquad (4.23) $$

**Step 1b** Using (4.22) to rewrite the left-hand side of equation (4.23), we obtain,

$$ \mathfrak{B}_{-} \left( \begin{aligned} & -\mathfrak{B}_{-} (\mathbf{q}^{*-} - \mathbf{q}^{-}) - c_{qS2}^{-} \mathbf{q}^{*-} + c_{qP}^{-} \mathbf{q}^{-} \\ & + \left( -(c_{qP}^{-})^2 + (c_{qS1}^{-})^2 \right) \mathbf{q}^{a-} + \left( -(c_{qS1}^{-})^2 + (c_{qS2}^{-})^2 \right) \mathbf{q}^{b-} - (c_{qS2}^{-})^2 \mathbf{q}^{*-} + (c_{qP}^{-})^2 \mathbf{q}^{-} \end{aligned} \right) = c_{qP}^{-} \mathfrak{B}_{-} \mathbf{q}^{-} - c_{qS2}^{-} \mathfrak{B}_{-} \mathbf{q}^{*-} \qquad (4.24) $$

After simplification, we arrive at,

$$ -\mathfrak{B}_{-}^{2} (\mathbf{q}^{*-} - \mathbf{q}^{-}) = \left( -(c_{qP}^{-})^2 + (c_{qS1}^{-})^2 \right) \mathbf{q}^{a-} + \left( -(c_{qS1}^{-})^2 + (c_{qS2}^{-})^2 \right) \mathbf{q}^{b-} \\ - (c_{qS2}^{-})^2 \mathbf{q}^{*-} + (c_{qP}^{-})^2 \mathbf{q}^{-}. \qquad (4.25) $$

Equations (4.22) and (4.25) give a linear system which determines uniquely $(\mathbf{q}^{a-}, \mathbf{q}^{b-})$ in terms of $\mathbf{q}^{*-}$ and $\mathbf{q}^{-}$,

$$ \begin{pmatrix} \alpha & \beta \\ \alpha\tilde{\alpha} & \beta\tilde{\beta} \end{pmatrix} \begin{pmatrix} \mathbf{q}^{a-} \\ \mathbf{q}^{b-} \end{pmatrix} = \begin{pmatrix} -\mathfrak{B}_{-} (\mathbf{q}^{*-} - \mathbf{q}^{-}) - c_{qS2}^{-} \mathbf{q}^{*-} + c_{qP}^{-} \mathbf{q}^{-} \\ \mathfrak{B}_{-}^{2} (\mathbf{q}^{*-} - \mathbf{q}^{-}) - (c_{qS2}^{-})^2 \mathbf{q}^{*-} + (c_{qP}^{-})^2 \mathbf{q}^{-} \end{pmatrix}, \qquad (4.26) $$

$$ \text{where } \alpha = c_{qP}^{-} - c_{qS1}^{-}, \quad \tilde{\alpha} = c_{qS1}^{-} + c_{qP}^{-}, \quad \beta = c_{qS1}^{-} - c_{qS2}^{-}, \quad \tilde{\beta} = c_{qS1}^{-} + c_{qS2}^{-}. \qquad (4.27) $$

Note that the determinant of the coefficient matrix in (4.26) is

$$ \alpha\beta(\tilde{\beta} - \tilde{\alpha}) = (c_{qP}^{-} - c_{qS1}^{-}) (c_{qS1}^{-} - c_{qS2}^{-}) (c_{qS2}^{-} - c_{qP}^{-}). \qquad (4.28) $$

This matrix is thus invertible under the assumption of distinct wave speeds (4.13). In this case, we obtain the expression of $\mathbf{q}^{a-}$ and $\mathbf{q}^{b-}$ in terms of $\mathbf{q}^{*-}$ and $\mathbf{q}^{-}$,

$$ \begin{pmatrix} \mathbf{q}^{a-} \\ \mathbf{q}^{b-} \end{pmatrix} = \frac{1}{\tilde{\beta} - \tilde{\alpha}} \begin{pmatrix} \frac{\tilde{\beta}}{\alpha} & -\frac{1}{\alpha} \\ -\frac{\tilde{\alpha}}{\beta} & \frac{1}{\beta} \end{pmatrix} \begin{pmatrix} -\mathfrak{B}_{-} (\mathbf{q}^{*-} - \mathbf{q}^{-}) - c_{qS2}^{-} \mathbf{q}^{*-} + c_{qP}^{-} \mathbf{q}^{-} \\ \mathfrak{B}_{-}^{2} (\mathbf{q}^{*-} - \mathbf{q}^{-}) - (c_{qS2}^{-})^2 \mathbf{q}^{*-} + (c_{qP}^{-})^2 \mathbf{q}^{-} \end{pmatrix}. \qquad (4.29) $$
From here, we can proceed by working with either expression of $\mathbf{q}^{a-}$ and (4.18a), or with $\mathbf{q}^{b-}$ and (4.18c). The first option is chosen in the following step,

$$
\begin{align}
\alpha(\tilde{\beta} - \tilde{\alpha})\mathbf{q}^{a-} &= \tilde{\beta} \left( -\mathfrak{B}_{-}(\mathbf{q}^{*-} - \mathbf{q}^{-}) - c_{qS2}^{-} \mathbf{q}^{*-} + c_{qP}^{-} \mathbf{q}^{-} \right) \nonumber \\
&\quad - \left( \mathfrak{B}_{-}^{2}(\mathbf{q}^{*-} - \mathbf{q}^{-}) - (c_{qS2}^{-})^{2} \mathbf{q}^{*-} + (c_{qP}^{-})^{2} \mathbf{q}^{-} \right) \nonumber \\
&= -\left( \tilde{\beta}\mathfrak{B}_{-} + \mathfrak{B}_{-}^{2} \right) (\mathbf{q}^{*-} - \mathbf{q}^{-}) + c_{qS2}^{-} \left( -\tilde{\beta} + c_{qS2}^{-} \right) \mathbf{q}^{*-} + c_{qP}^{-} \left( \tilde{\beta} - c_{qP}^{-} \right) \mathbf{q}^{-}. \tag{4.30}
\end{align}
$$

With some algebraic manipulations<fnref n="3" />, we obtain the following identity,

$$ c_{qP}^{-} (\tilde{\beta} - c_{qP}^{-}) - \alpha(\tilde{\beta} - \tilde{\alpha}) = c_{qS1}^{-} c_{qS2}^{-}. \quad (4.32) $$

We employ this to further rewrite (4.30) as,

$$ \alpha(\tilde{\beta} - \tilde{\alpha}) (\mathbf{q}^{a-} - \mathbf{q}^{-}) = - (\tilde{\beta}\mathfrak{B}_{-} + \mathfrak{B}_{-}^{2} + c_{qS1}^{-} c_{qS2}^{-}) (\mathbf{q}^{*-} - \mathbf{q}^{-}). \quad (4.33) $$

**Step 2a** Using (4.18a), we have the following equalities,

$$ (\mathfrak{B}_{-} + c_{qP}^{-}) (\mathbf{q}^{a-} - \mathbf{q}^{-}) = 0 \quad \Rightarrow \quad (\mathfrak{B}_{-} + c_{qP}^{-}) \alpha(\tilde{\beta} - \tilde{\alpha}) (\mathbf{q}^{a-} - \mathbf{q}^{-}) = 0. \quad (4.34) $$

Substitute the difference given by (4.33) into the above expression to obtain,

$$ (\mathfrak{B}_{-} + c_{qP}^{-}) (\tilde{\beta}\mathfrak{B}_{-} + \mathfrak{B}_{-}^{2} + c_{qS1}^{-} c_{qS2}^{-}) (\mathbf{q}^{*-} - \mathbf{q}^{-}) = 0. \quad (4.35) $$

In rewriting left-hand side and introducing $\mathbf{p}_2^-$ and $\mathbf{p}_3^-$,

$$ \mathbf{p}_3^- := c_{qS1}^- c_{qS2}^- c_{qP}^-, \quad \mathbf{p}_2^- := c_{qS1}^- c_{qS2}^- + c_{qS1}^- c_{qP}^- + c_{qS2}^- c_{qP}^- = c_{qS1}^- c_{qS2}^- + c_{qP}^- \tilde{\beta}, \quad (4.36) $$

we rewrite the above equation as,

$$ (\mathfrak{B}_{-}^{3} + \mathbf{p}_{2}^{-} \mathfrak{B}_{-} + (c_{qS1}^{-} + c_{qS2}^{-} + c_{qP}^{-}) \mathfrak{B}_{-}^{2} + \mathbf{p}_{3}^{-} \mathfrak{B}_{-}) (\mathbf{q}^{*-} - \mathbf{q}^{-}) = 0. \quad (4.37) $$

**Step 2b** By algebraic computation, cf. [31], we have the following identities for the action of powers of $\mathfrak{B}$ on vector $\mathbf{q} = \begin{pmatrix} \vec{\epsilon} \\ \rho \mathbf{v} \end{pmatrix}$,

$$ \mathfrak{B}(\nu)\mathbf{q} = - \begin{pmatrix} \mathfrak{A}^t(\nu)\mathbf{v} \\ \sigma\nu \end{pmatrix}, \quad \mathfrak{B}(\nu)^2\mathbf{q} = \begin{pmatrix} \rho^{-1}\mathfrak{A}^t(\nu)\sigma\nu \\ \mathbf{\Gamma}(\nu)\mathbf{v} \end{pmatrix}, \quad \mathfrak{B}(\nu)^3\mathbf{q} = -\frac{1}{\rho} \begin{pmatrix} \mathfrak{A}^t(\nu)\mathbf{\Gamma}(\nu)\mathbf{v} \\ \mathbf{\Gamma}(\nu)\sigma\nu \end{pmatrix}. \quad (4.38) $$

$$ \begin{align} (c_{qS1}^{-} - c_{qP}^{-})(c_{qS2}^{-} - c_{qP}^{-}) &= c_{qS1}^{-}c_{qS2}^{-} - c_{qP}^{-}(c_{qS1}^{-} + c_{qS2}^{-} - c_{qP}^{-}) \tag{4.31} \\ \Rightarrow \quad c_{qP}^{-}(\tilde{\beta} - c_{qP}^{-}) - \alpha(\tilde{\beta} - \tilde{\alpha}) &= c_{qP}^{-}(c_{qS1}^{-} + c_{qS2}^{-} - c_{qP}^{-}) - (c_{qP}^{-} - c_{qS1}^{-})(c_{qS2}^{-} - c_{qP}^{-}) = c_{qS1}^{-}c_{qS2}^{-}. \notag \end{align} $$

<footnote n="3">³This is seen as,</footnote>
We employ these identities with $\mathfrak{B}_-$ in (4.37), whose last three components give,

$$ -\left(\frac{\Gamma(\nu)^-}{\rho^-} + \mathfrak{p}_2^-\right) ((\sigma\nu)^{*-} - (\sigma\nu)^-) + \rho^- \left(\left(c_{qS1}^- + c_{qS2}^- + c_{qP}^-\right) \frac{\Gamma(\nu)^-}{\rho^-} + \mathfrak{p}_3^-\right) (\mathbf{v}^{*-} - \mathbf{v}^-) = 0 \\ \Rightarrow (\sigma\nu)^{*-} - (\sigma\nu)^- = \rho^- \left(\frac{\Gamma(\nu)^-}{\rho^-} + \mathfrak{p}_2^-\right)^{-1} \left(\left(c_{qS1}^- + c_{qS2}^- + c_{qP}^-\right) \frac{\Gamma(\nu)^-}{\rho^-} + \mathfrak{p}_3^-\right) (\mathbf{v}^{*-} - \mathbf{v}^-). \quad (4.39) $$

With some algebraic manipulations, we arrive at

$$ (\sigma\nu)^{*-} - (\sigma\nu)^- = \rho^- \left(c_{qS1}^- + c_{qS2}^- + c_{qP}^-\right) \left(1 + \gamma^- \left(\frac{\Gamma(\nu)^-}{\rho^-} + \mathfrak{p}_2^-\right)^{-1}\right) (\mathbf{v}^{*-} - \mathbf{v}^-), \quad (4.40) \\ \text{where } \gamma^- := \frac{\mathfrak{p}_3^-}{c_{qS1}^- + c_{qS2}^- + c_{qP}^-} - \mathfrak{p}_2^- $$

A similar relation is obtained for mesh cell $K^+$, with ‘-’ replaced by ‘+’ in the above expression.

**Derivation of HDG traces** From this point on, to arrive at the hybridized traces along $F^\pm$ (4.2), we follow the hybridization procedure employed in [36] for isotropy. First, we let the numerical trace along $F^\pm$ be given by the intermediate states,

$$ \hat{\mathbf{v}}^{\pm} := \mathbf{v}^{*\pm}, \quad \hat{\sigma\nu}^{\pm} := (\sigma\nu)^{*\pm}. \quad (4.41) $$

Next, employing the transmission condition (4.19) and relation (4.40), and introducing the quantity $\lambda_\nu$ defined on $F$ to represent $\mathbf{v}^{*\pm}$, the numerical traces (4.41) along $F^-$ are rewritten as,

$$ \hat{\mathbf{v}}^{-} = \lambda_{\mathbf{v}}, \quad \text{and} \quad \hat{\sigma\nu}^{-} = (\sigma\nu)^{-} + M_{\text{Gani}}(\nu)^{-} (\lambda_{\mathbf{v}} - \mathbf{v}^{-}), \quad (4.42) $$

$$ \text{with } M_{\text{Gani}}(\nu)^- := \rho^- \left( c_{qS1}^- + c_{qS2}^- + c_{qP}^- \right) \left( 1 + \gamma^- \left( \frac{\Gamma(\nu)^-}{\rho^-} + \mathfrak{p}_2^- \right)^{-1} \right). \quad (4.43) $$

The numerical traces $\hat{\mathbf{v}}^+$, $\hat{\sigma\nu}^+$ along $F^+$, are defined similarly with ‘-’ replaced by ‘+’ above.

# 5 Numerical experiments

The HDG method is implemented in the open-source parallel software hawen<fnref n="4" />, [19]. We carry out numerical experiments to evaluate the accuracy of numerical solutions depending on the choice of stabilization. Three families of stabilization are investigated: based on the identity matrix, the Kelvin-Christoffel (KC) matrix $\Gamma$ ((2.26) and (2.29)), and the Godunov matrix $M_{\text{Godunov}}$ ((3.51) and (4.43)). Each of them is first considered with a scaling factor $\tau$:

$$ \text{Choices of } \tau_u \text{ for} \quad \hat{\sigma\nu} = \sigma_h - \tau_u(\mathbf{u}_h - \lambda_{uh}), \quad \tau_u = \begin{cases} \text{Identity-based: } \tau_{u\mathbb{R}\pm}(\tau) = \pm\omega \tau \mathbb{1}_d, \\ \tau_{u\mathbb{1}\pm}(\tau) = \pm i\omega \tau \mathbb{1}_d; \\ \text{KC based: } \tau_{u\Gamma}(\tau) = -i\omega \tau \Gamma; \\ \text{Godunov based: } \tau_{uG}(\tau) = -i\omega \tau M_{\text{Godunov}}. \end{cases} \quad (5.1) $$

**Remark 6.** The scaling factor $i\omega$ is motivated by the relation between the velocity and the displacement, and we refer to our extended report [31] for more details. This result is also illustrated in the comparison within the identity-based family, with $\tau$ purely real or imaginary.

<footnote n="4">⁴https://ffaucher.gitlab.io/hawen-website/</footnote>
We consider isotropic and anisotropic materials. Our investigation consists of two steps:

1. In Subsection 5.1, we investigate the optimal choice for the scaling factor $\tau$ within each of the three families of stabilization in (5.1). The investigation is carried out with planewaves, and thus concerns single-typed waves propagating in a homogeneous material. Here, analytical solutions are used as references to evaluate the accuracy of the numerical simulations.

2. In Subsection 5.2, we compare the best candidates from each family in a highly heterogeneous medium with a Dirac point source. Here, the wavefield contains all types of waves and, as there are no analytical solutions, a reference solution is constructed by using a highly refined mesh and high-order polynomial basis functions. This experiment is carried out for multiple frequencies and orders of polynomial, in isotropic and TTI medium.

To evaluate the difference between a reference solution (either an analytic or a numerical one computed with a refined mesh) and simulations, we introduce relative errors $\epsilon$ and $\mathcal{E}$. With $w$ representing a component of the displacement field $\mathbf{u}$ or the stress tensor $\boldsymbol{\sigma}$, and $\mathbf{x}_k$ the $k^{\text{th}}$ position where the solutions are evaluated, we define,

$$ \epsilon[w] := \frac{1}{N} \sum_{k=1}^{N} \frac{|w^{\text{ref}}(\mathbf{x}_k) - w(\mathbf{x}_k)|}{\|w^{\text{ref}}\|}, \quad \text{with} \quad \|w^{\text{ref}}\| = \sqrt{\sum_{k=1}^{N} |w^{\text{ref}}(\mathbf{x}_k)|^2}. \quad (5.2) $$

The error for the total field is given by,

$$ \mathfrak{E}[\mathbf{u}] := \frac{1}{n_d} \sum_{j \in I} \mathfrak{E}[u_j]; \quad \mathfrak{E}[\boldsymbol{\sigma}] := \frac{1}{\mathfrak{N}_d} \sum_{j \in \mathcal{I}_{\text{sm}}} \mathfrak{E}[\boldsymbol{\sigma}_j], \quad (5.3) $$

with $n_d = 2$ or 3 for two and three dimensions respectively, and $\mathfrak{N}_d = 3$ or 6. In the experiments, we use a Cartesian grid for the positions $\mathbf{x}_k$ where the solutions are evaluated even though our computational mesh can be unstructured. We however exclude the positions near the boundaries, and near Dirac point source (i.e., for experiments of Subsection 5.2).

## 5.1 3D experiments with planewaves

We consider a three-dimensional domain $\Omega$ corresponding to the cube $(-1, 1)^3$, with boundary denoted by $\partial\Omega$. We consider equation (3.1) on $\Omega$ with boundary condition $\mathbf{u} = \mathbf{u}_{\text{pw}}$ on $\partial\Omega$. We consider different types of planewaves and refer to [31, Appendix A] for more details on the derivation of planewaves in linear elasticity.

### 5.1.1 Elastic isotropic medium

We consider the propagation of P- and S-planewaves, with corresponding Dirichlet condition $\mathbf{u}_{\text{pw}}^{\text{P}}$ and $\mathbf{u}_{\text{pw}}^{\text{S}}$ respectively. We select the planewave direction $\mathbf{d} = (1/\sqrt{2}, 0, 1/\sqrt{2})^t$, and further impose a strong contrast between the P- and S-wavespeeds, with the following configuration:

$$ \left\{ \begin{array}{l} c_P = 2.5 \times 10^{-3} \, \text{m s}^{-1}, \quad c_S = 10^{-4} \, \text{m s}^{-1}, \quad \rho = 1 \, \text{kg m}^{-3}, \\ \Leftrightarrow \lambda = 6.24 \times 10^{-6} \, \text{Pa} \quad \text{and} \quad \mu = 10^{-8} \, \text{Pa}; \\ \text{P-planewave propagation:} \quad \mathbf{u}_{\text{pw}}^{\text{P}}(\mathbf{x}) = (1, 0, 1)^t e^{i \frac{\omega}{c_P} (\mathbf{d} \cdot \mathbf{x})}, \\ \text{S-planewave propagation:} \quad \mathbf{u}_{\text{pw}}^{\text{S}}(\mathbf{x}) = (0, 1, 0)^t e^{i \frac{\omega}{c_S} (\mathbf{d} \cdot \mathbf{x})}. \end{array} \right. \quad (5.4) $$
The simulations use a mesh consisting of 48 000 tetrahedra and polynomials of order 4. Due to the strong contrast between the P- and S-wavelength, we use frequency $\omega/(2\pi) = 8$mHz for the P-planewave, and $\omega/(2\pi) = 0.4$mHz for the S-planewave, resulting in wavelength of size 0.3m and 0.25m, respectively.

**Remark 7.** Our experiment (5.4) on a unitary cube is equivalent to working with a cube ($-1,1$)km$^3$ with $c_P = 2500$m s$^{-1}$ and $c_S = 100$m s$^{-1}$, using P-planewave frequency 8Hz and S-planewave frequency 0.4Hz. This is obtained by maintaining the number of wavelengths propagating in a domain of size $L^3$, which is given by $(L \times \text{frequency})/\text{wavespeed}$.

In Figure 2, we investigate the accuracy of the solution as the scaling factor $\tau$ in (5.1) varies. For the identity-based stabilization, we investigate the scaling factor on the purely real and complex lines and determine the optimal sign. For the sake of conciseness, the relative error $\epsilon$ is shown only for some of the non-zero wave fields: $u_x$ and $\sigma_{zz}$ for the P-planewave, and $u_y$ for the S-planewave.

<figure><img src="image_3.png" /><figcaption>Figure 2: Investigation of the optimal scaling $\tau$ within the three families of stabilization in (5.1) for planewave propagation in isotropic medium using parameters (5.4). Our criterion of accuracy is the relative error $\epsilon$ (5.2).</figcaption></figure>

We can draw the following observations.

- From the plots using the identity-based stabilization, Figure 2a, it is clear that the real-valued stabilization $\tau_{u\mathbb{R}^+}$ is inaccurate compared to the imaginary one $\tau_{u\mathbb{R}^-}$. On the other hand, the choice of sign ($\pm i$) barely affects the accuracy.

- The optimal value of the coefficient $\tau$ in the identity-based stabilization $\tau_{u\mathbb{R}\pm}(\tau)$ (5.1) depends on the type of waves propagating and appears to be given by the impedance.
Namely, for P-planewave, the optimal coefficient would be the P-impedance, while it is the S-impedance for the S-planewave, see Figure 2a.

- For the Kelvin-Christoffel stabilization, Figure 2c, the optimal scaling coefficient $\tau$ is proportional to the inverse of the velocity, that is, the slowness. It seems to depend only slightly on the type of waves. For P-planewave, the value of the optimal scaling factor is hard to identify and lies between $\tau = 1/c_P$ and $\tau = 1/c_S$. The S-planewave shows more flexibility, with the optimal value for the scaling given in range $\tau \in (1/c_P, 1/c_S)$.

- For the Godunov stabilization, Figure 2b, the optimal scaling coefficient does not depend on the type of waves, and is simply given by $\tau = 1$.

- Comparing the accuracy of the different families of stabilization, the Godunov stabilization (with $\tau = 1$) and the Kelvin-Christoffel one (with appropriate slowness scaling) give similar accuracy. The identity-based stabilization, although close, does not match this accuracy.

Therefore, the Godunov stabilization appears as the most versatile choice as it accurately treats both the P-planewave and the S-planewave, without any specific scaling coefficient, with $\tau = 1$. For the identity-based, the optimal scaling factor clearly depends on the types of waves. For the Kelvin-Christoffel stabilization, the same value of the scaling factor can be chosen for both types of waves. This can be expected since the Kelvin-Christoffel matrix contains the information of the material.

### 5.1.2 Elastic vertical transverse isotropic medium

We modify the previous experiment to include anisotropy in the medium. We work with Thom- sen's parameters, [37], with notation of [31, Section 3.5.2] for VTI, and select<fnref n="5" />

$$ \lambda^{\text{TI}} = 3.6107 \times 10^{-6} \text{Pa}, \quad \mu^{\text{TI}} = 3.0490 \times 10^{-9} \text{Pa}, \quad \rho = 1 \text{kg m}^{-3}, \\ \epsilon = 1.12, \quad \delta = -0.235 \quad \gamma = 2.28. \qquad (5.5) $$

This amounts to the following values of the stiffness tensor coefficients:

$$ \begin{gathered} C_{11} = C_{22} = 1.1719 \times 10^{-5}, \quad C_{33} = 3.6168 \times 10^{-6}, \quad C_{44} = C_{55} = 3.0490 \times 10^{-9}, \\ C_{66} = 1.6953 \times 10^{-8}, \quad C_{12} = C_{11} - 2C_{66} = 1.1685 \times 10^{-5}, \quad C_{23} = C_{23} = 2.6268 \times 10^{-6}. \end{gathered} $$

In direction $\mathbf{d} = (d_x, 0, d_z)$, the qP- and sH-planewaves are given by, [31, Appendix A],

$$ \mathbf{u}_{\text{pw}}^{\text{qP}}(\mathbf{x}) = \frac{e^{i \frac{\omega}{c_{\text{qP}}} (\mathbf{d} \cdot \mathbf{x})}}{\sqrt{(C_{11} + C_{55})d_x^2 + (C_{55} + C_{33})d_z^2 - 2\rho c_{\text{qP}}}} \begin{pmatrix} \sqrt{C_{55}d_x^2 + C_{33}d_z^2 - \rho c_{\text{qP}}^2} \\ 0 \\ \sqrt{C_{11}d_x^2 + C_{55}d_z^2 - \rho c_{\text{qP}}^2} \end{pmatrix}, \quad (5.7a) \\ \mathbf{u}_{\text{pw}}^{\text{sH}}(\mathbf{x}) = (0, 1, 0)^t e^{i \frac{\omega}{c_{\text{sH}}} (\mathbf{d} \cdot \mathbf{x})}, \quad (5.7b) $$

<footnote n="5">⁵The values of the anisotropic coefficients $\epsilon$, $\delta$ and $\gamma$ correspond to muscovite crystal in [37]. The values of $\lambda^{\text{TI}}$ and $\mu^{\text{TI}}$ are selected such that the qP- and sH-wavespeeds in (5.8) correspond to the P- and S-wavespeed of the isotropic experiment (5.4), i.e., $2.5 \times 10^{-3}$ and $10^{-4}$ m s$^{-1}$ respectively.</footnote>
with corresponding qP- and sH-wavespeeds, denoted by $c_{qP}$ and $c_{sH}$,

$$c_{qP}^2(d_x, 0, d_z) = \frac{1}{2\rho} \left( C_{11}d_x^2 + C_{33}d_z^2 + C_{55} + \sqrt{((C_{11} - C_{55})d_x^2 + (C_{55} - C_{33})d_z^2)^2 + 4(C_{13} + C_{55})^2 d_x^2 d_z^2} \right), \quad (5.8a)$$

$$c_{sH}^2(d_x, 0, d_z) = \frac{C_{66}d_x^2 + C_{55}d_z^2}{\rho}. \quad (5.8b)$$

In the experiment, we work with direction $d = (1/\sqrt{2}, 0, 1/\sqrt{2})^t$, and illustrate the corresponding analytic solutions in Figure 3.

<figure><img src="image_4.png" /><figcaption>Figure 3: Reference solution for the 3D planewave propagation in VTI medium (5.5).</figcaption></figure>

In Figure 4, we evaluate the accuracy of the solutions depending on the scaling of the coefficient $\tau$ for the three families of stabilization in (5.1). The qP-planewave is computed at frequency 8 mHz and the sH-planewave at frequency 0.4 mHz. The observations for the elastic VTI medium are close to those of the elastic isotropic case of Subsection 5.1.1:

- The optimal value of the scaling coefficient for the Godunov stabilization is $\tau = 1$, and it does not depend on the type of wave propagating.

- The optimal value of the scaling for the Kelvin-Christoffel stabilization is in between $1/c_{qP}$ and $1/c_{sH}$, with the value of the optimal scaling hard to find for the qP-planewave.

- For the identity-based stabilization, the optimal coefficient (which is imaginary) depends on the type of waves and is given by the impedances. Furthermore, the accuracy of the results, even with the optimal scaling, is less than the accuracy obtained using other stabilizations.

In short, the Godunov stabilization is the most robust to treat both types of waves, with its default scaling factor, $\tau = 1$, being optimal.
<figure><img src="image_5.png" /><figcaption>(a) Identity-based stabilization using real and imaginary scaling factors. Relative error for P-planewave $\epsilon[u_x]$ (left) and $\epsilon[\sigma_{zz}]$ (middle); for S-planewave $\epsilon[u_y]$ (right).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(b) Using Godunov stabilization.</figcaption></figure>

<caption>(c) Using Kelvin-Christoffel stabilization.</caption>

<caption>Figure 4: Investigation of the optimal scaling $\tau$ within the three families of stabilization in (5.1) for planewave propagation in VTI medium with parameters (5.5). Our criterion of accuracy is the relative error $\epsilon$ (5.2).</caption>

## 5.2 2D experiments with highly varying properties and point-source

To study further the efficiency of stabilizations, we consider a heterogeneous medium on a two-dimensional disk domain $\Omega$. The variation of the properties is too strong for a piecewise-constant representation to be reliable and we instead use piecewise-polynomial representations, whose implementation is straightforward with the formulation based on the compliance tensor, (3.1). In this experiment, the piecewise polynomial representation is independent on each cell, and allows discontinuities across the interfaces. We consider a problem with a Dirac point-source,

$$ \{-\omega^2 \rho \mathbf{u} - \nabla \cdot \mathbf{\sigma} = \delta_y, \quad \mathbf{S} \mathbf{\sigma} = \nabla^5 \mathbf{u}, \quad \text{in } \Omega, \quad (5.9a) $$

$$ \mathbf{\sigma} \nu = -i\omega\mathcal{Z}\mathbf{u}, \quad \text{on } \partial\Omega. \quad (5.9b) $$

Its solution is a superposition of all types of waves supported by the medium. In (5.9), $\mathbf{y}$ is the position of the source (the center of the disk in our experiments), and the impedance coefficient $\mathcal{Z}$ corresponds to a low-order absorbing boundary condition (e.g., (5.11)).
Following our previous experiments to select scaling, we compare the following stabilizations:

$$
\begin{align}
\text{Identity-based:} \quad & \tau_{\mathbf{u}I_s} := \tau_{\mathbf{u}\parallel-}(\rho c_{qS2}) &&= -i\omega \rho c_{qS2} \parallel d; \nonumber \\
& \tau_{\mathbf{u}I} := \tau_{\mathbf{u}\parallel-}(1) &&= -i\omega \parallel d; \nonumber \\
\text{KC based:} \quad & \tau_{\mathbf{u}\Gamma_p} := \tau_{\mathbf{u}\Gamma}(1/c_{qP}) &&= -\frac{i\omega}{c_{qP}} \Gamma; \tag{5.10} \\
& \tau_{\mathbf{u}\Gamma_1} := \tau_{\mathbf{u}\Gamma}(1) &&= -i\omega \Gamma; \nonumber \\
\text{Godunov based:} \quad & \tau_{\mathbf{u}G_1} := \tau_{\mathbf{u}G}(1) &&= -i\omega M_{\text{Godunov}}. \nonumber
\end{align}
$$

**Remark 8.** We also investigated the qP-wavespeed in the identity-based and the sH-wavespeed for the Kelvin-Christoffel stabilizations, denoted respectively as $\tau_{\mathbf{u}\parallel-}$ and $\tau_{\mathbf{u}\Gamma}$ in (5.10). They give (slightly) worse results, and are not shown for the sake of conciseness. For the identity-based stabilization, the choice of the S-wavespeed for $\tau$ is further discussed in Subsection 5.3.

### 5.2.1 Numerical representation of the models of parameters

One advantage in writing the elastic system in terms of the compliance matrix $S$ rather than the stiffness tensor $C$ is that it facilitates a non-constant description of background materials within a mesh cell. For this experiment, we construct a background model inspired from model $S$, a standard solar model in helioseismology, [9, 2]. Model $S$ provides us with radially symmetric density and wave speed, which is employed as a P-wave speed for our model. Here, we consider a two-dimensional disk geometry of radius 1, on which the radial profiles are defined.

This model is particularly challenging because of the high variation of amplitude, with an exponential decrease of the density close to the boundary, and wave speed variation of about two orders of magnitude, [2, Figure 3]. In order to capture correctly the variation of the models, using a piecewise-constant model representation would require extreme mesh refinement near the boundary. Exploiting the flexibility of the formulation with compliance tensor (3.1), we instead represent the physical parameters (wave speeds and density) using a basis of Lagrange polynomials on each cell. In Figure 5 we show the mesh of the 2D disk employed in our simulations. In the same figure, we compare the two options of representing the density model near the surface: piecewise constant and piecewise polynomial (here with order 2 Lagrange polynomials) representations. The Lagrange polynomial representation preserves the spherical nature of the model, which is not the case for the piecewise constant representation. Note that we still have a finer discretization near the surface, to ensure accuracy as the wavelength drastically decreases (as the wave speed drops), and to have an accurate circular geometry.

### 5.2.2 Elastic isotropic medium

We first consider an isotropic medium, with the P-wavespeed and density given by the solar background, and S-wavespeed selected such that $c_S = 0.70c_P$. The Dirac point-source is positioned at the center of the unit disk, and absorbing conditions are implemented on the boundary, [25, 27], such that

$$ \sigma \cdot \nu - i\omega \left( \sqrt{(\lambda + 2\mu)\rho} \nu \otimes \nu + \sqrt{\mu\rho} (\parallel d - \nu \otimes \nu) \right) \mathbf{u} = 0, \quad \text{on } \partial\Omega. \quad (5.11) $$

Since there are no analytic solutions in this case, to investigate accuracy, a reference solution is constructed with a refined mesh of 80 000 cells (while for the simulation we have 50 000 cells, Figure 5a), polynomials of order 7, and with stabilization $\tau_{\mathbf{u}G_1}$. The reference solutions are
<figure><img src="image_7.png" /><figcaption>Figure 5: Computational mesh and solar-like density represented on the mesh near surface. The piecewise-constant representation is unable to preserve the radial symmetry, as shown in the last (radial) layers.</figcaption></figure>

shown in Figure 6 where we also compare with the solution obtained with a piecewise-constant model representation. Due to the exponentially decreasing nature of the density, we scale the displacement fields by $\sqrt{\rho}$ to give a better visualization of the solution, [2, 3]. For a similar reason, the components of the stress tensor are scaled by $1/\sqrt{\rho}$.

<figure><img src="image_8.png" /><figcaption>Figure 6: Comparison of the reference solutions at 4 mHz of the solar-like isotropic elastic background models, obtained using piecewise-constant representation (bottom) or with piecewise polynomial representation (Lagrange basis of order 2) (top). The computations use HDG formulation $(\mathbf{u}, \boldsymbol{\sigma})_S$ with stabilization $\tau_{\mathbf{u}G_1}$ and polynomial order 7.</figcaption></figure>

The solutions corresponding to a piecewise-constant model representation (both displacement **u** and stress tensor **σ**) show ripples and artifacts, in addition, their radial nature is not preserved. On the contrary, results employing the Lagrange basis model representation capture well the
spherical pattern and provide ‘clean’ solutions. With such high variations in the background models, we see that it is mandatory to design an efficient representation, and that piecewise-constant model representation is not appropriate. In this way, the HDG formulation based on the compliance tensor $\mathbf{S}$ is useful as it allows us to vary easily the models within each cell.

In Figure 7, we plot the relative error $\mathcal{E}$ (5.3) in terms of the frequency and polynomial orders, for the stabilization coefficients of (5.10). Note that for some frequencies and orders, stabilization $\tau_{ul}$ results in an ill-conditioned matrix that cannot be factorized or leads to erroneous solutions. In these cases, the discretized wave system with $\tau_{ul}$ cannot be solved. However, when it works, using $\tau_{ul}$ still gives worst performance with high levels of error.

- We observe that the error increases with increasing frequency, Figure 7a, which is common as the wavelength is reduced.

- As expected, the error also decreases with increasing polynomial orders, cf. Figure 7b.

- The stabilization $\tau_{uG_1}$ gives accurate results, with an accuracy that can be matched by using the appropriate scaling in other stabilizations, with either $\tau_{ul_s}$ or $\tau_{u\Gamma_p}$.

Overall, this test confirms that $\tau_{uG_1}$ is the most versatile choice. While an equivalent accuracy can be obtained with other stabilizations, it requires a judicial choice of scaling.

<figure><img src="image_9.png" /><figcaption>Figure 7: Relative error $\mathcal{E}$ of 5.3 with frequency (left) and polynomial order (right) for heterogeneous elastic isotropic medium, with the stabilization coefficients of (5.10).</figcaption></figure>

### 5.2.3 Elastic tilted transverse isotropic medium

We now consider the propagation of waves in a tilted transverse isotropic (TTI) medium, which is a rotated version of a VTI one by an angle $\theta$, cf. [31, Section 3.5.3]. The anisotropy is represented with constant Thomsen's parameters $\epsilon = 0.25$, $\delta = 0.15$, and $\theta = 45^\circ$. Note that compared to [31, Section 3.5.3] written for 3D, the 2D case only has one angle, and no anisotropic coefficient $\gamma$. Here, the heterogeneous solar-like background velocity $c$ and density $\rho$ are used to define $\mu^{\text{TI}} = (0.70c)^2\rho$ and $\lambda^{\text{TI}} = c^2\rho - 2\mu^{\text{TI}}$. In Figure 8, we show the reference solution at frequency 4 mHz, computed on a mesh with 80 000 triangles (while following simulations use 50 000 elements), with polynomial approximation of order 7 and Godunov stabilization. Compared to the isotropic case of Figure 6, the solution is less smooth and contains more ripples. Note that for the TTI absorbing boundary conditions, we follow [4]. It is out of the scope to discuss the accuracy of the anisotropic boundary conditions here.

In Figure 9, we plot the relative error $\mathcal{E}$ depending on the frequency and polynomial orders. Similar to the isotropic case, stabilization $\tau_{ul}$ sometimes results in an ill-conditioned matrix
<figure><img src="image_10.png" /></figure>

Figure 8: Reference solutions for the elastic TTI wave equation in a solar-like background model at 4 mHz obtained with stabilization $\tau_{uG_1}$ and polynomial approximation of order 7. The background models are represented by piecewise polynomials on each cell with Lagrange basis of order 2.

leading to a linear system that cannot be solved. The results confirm the behaviour observed in the isotropic case: the Godunov stabilization gives the most accurate results, and the level of accuracy can be met with other stabilization involving a well-chosen scaling factor.

<figure><img src="image_11.png" /></figure>

Figure 9: Relative error $\mathcal{E}$ of (5.3) with frequency (left) and the polynomial order (right) for heterogeneous elastic TTI medium, with stabilization coefficients of (5.10).

## 5.3 Summary of numerical experiments

In our first experiments working with planewaves, we identify the optimal scaling factor for each family of stabilization in (5.1). We find that the Godunov stabilization is the most robust as it does not need scaling and $\tau = 1$ gives the best results for all types of waves. On the other hand, the Kelvin–Christoffel and identity-based families need a suitable scaling factor in order to reach the same level of accuracy obtained with the Godunov matrix. The optimal value for the scaling does not seem to change for the Kelvin–Christoffel stabilization and leans towards the P-slowness. However, the identity-based stabilization shows a different optimal scaling factor depending on the type of waves propagating. Our experiments also compare the performance of the stabilizations in both isotropy and anisotropy, for which the same conclusions are obtained.

In our second experiment, we consider waves generated by a point source in a heterogeneous medium, thus a superposition of different types of waves. In this case, the Godunov stabilization is the most accurate. The identity-based stabilization with the S-impedance scaling factor and
and the Kelvin-Christoffel stabilization with the P-slowness also yield comparable accuracy. In the identity-based family, the better performance of the S-impedance compared to the P-impedance can be explained by the equipartition phenomenon, [33, 34], which says that S-waves are more energetic than P-waves. In another word, although all types of waves are present, the magnitude of S-waves is stronger, hence stabilization with S-impedance provides higher accuracy.

As a common feature among the optimal form from each family, the components of the stabilization operators take the value of an impedance (i.e., density multiplied by wavespeed). This is evident for the identity-based and Godunov stabilization (cf. (3.51) and (4.43)). For the Kelvin–Christoffel stabilization, the components of the matrix Γ take the form of ρc² (2.26), and the optimal scaling τ is a slowness, hence τΓ has a magnitude of an impedance.

# 6 Conclusion

In this work, we employed Voigt notation in the HDG method to describe compactly the discrete problem for anisotropic elasticity. Additionally, a first-order formulation with the compliance tensor is used and allows for mesh-wise variation of parameters. This, together with Voigt’s notation which provides efficient bookkeeping of physical parameters, will be indispensable in quantitative reconstruction of elastic parameters in inverse problems. Secondly, to determine an optimal choice of stabilization, we constructed the hybridized Godunov-upwind flux for anisotropic elasticity, which offers a versatile choice and removes the need for scaling factor tuning. It is worth noting that the identity-based stabilization, which is very popular, lacks a universal scaling factor, making it less robust for wave simulation in a general setting. We have carried out numerical experiments in two and three dimensions, considering isotropic elasticity and anisotropy, with constant backgrounds as well as one containing high variation and contrast. They demonstrate the performance and versatility of the Godunov stabilization, which is well suited for generic anisotropic material and different types of waves.

# Acknowledgments

This project was provided with computer and storage resources by GENCI at CINES thanks to the grant gda2306 on the supercomputer Adastra’s GENOA partition. This work was partially supported by the EXAMA (Methods and Algorithms at Exascale) project under grant ANR-22-EXNU-0002. FF acknowledges funding by the European Union with ERC Project INCORWAVE - grant 101116288. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency (ERCEA). Neither the European Union nor the granting authority can be held responsible for them.

# References

[1] D. N. ARNOLD, F. BREZZI, B. COCKBURN, AND L. D. MARINI, *Unified analysis of discontinuous Galerkin methods for elliptic problems*, SIAM journal on numerical analysis, 39 (2002), pp. 1749–1779.

[2] H. BARUCQ, F. FAUCHER, D. FOURNIER, L. GIZON, AND H. PHAM, *Efficient and accurate algorithm for the full modal Green's kernel of the scalar wave equation in helioseismology*, SIAM Journal on Applied Mathematics, 80 (2020), pp. 2657–2683.
[3] H. BARUCQ, F. FAUCHER, AND H. PHAM, Outgoing solutions and radiation boundary conditions for the ideal atmospheric scalar wave equation in helioseismology, ESAIM: Mathematical Modelling and Numerical Analysis, 54 (2020), pp. 1111–1138.

[4] M. BONNASSE-GAHOT, H. CALANDRA, J. DIAZ, AND S. LANTERI, *Hybridizable discontinuous Galerkin method for the 2-D frequency-domain elastic wave equations*, Geophysical Journal International, 213 (2018), pp. 637–659.

[5] T. BUI-THANH, *From Godunov to a unified hybridized discontinuous Galerkin framework for partial differential equations*, Journal of Comp. Physics, 295 (2015), pp. 114–146.

[6] T. BUI-THANH, *From Rankine-Hugoniot condition to a constructive derivation of HDG methods*, in Spectral and High Order Methods for Partial Differential Equations ICOSAHOM 2014, Springer, 2015, pp. 483–491.

[7] T. BUI-THANH, *Construction and analysis of HDG methods for linearized shallow water equations*, SIAM Journal on Scientific Computing, 38 (2016), pp. A3696–A3719.

[8] J. M. CARCIONE, *Wave fields in real media: Wave propagation in anisotropic, anelastic, porous and electromagnetic media*, Elsevier, 2007.

[9] J. CHRISTENSEN-DALSGAARD, W. DÄPPEN, S. V. AJUKOV, E. R. ANDERSON, H. M. ANTIA, S. BASU, V. A. BATURIN, G. BERTHOMIEU, B. CHABOYER, S. M. CHITRE, A. N. COX, P. DEMARQUE, J. DONATOWICZ, W. A. DZIEMBOWSKI, M. GABRIEL, D. O. GOUGH, D. B. GUENTHER, J. A. GUZIK, J. W. HARVEY, F. HILL, G. HOUDEK, C. A. IGLESIAS, A. G. KOSOVICHEV, J. W. LEIBACHER, P. MOREL, C. R. PROFFITT, J. PROVOST, J. REITER, E. J. RHODES, F. J. ROGERS, I. W. ROXBURGH, M. J. THOMPSON, AND R. K. ULRICH, *The current state of solar modeling*, Science, 272 (1996), pp. 1286–1292. 10.1126/science.272.5266.1286.

[10] B. COCKBURN, *Static condensation, hybridization, and the devising of the HDG methods*, Building bridges: connections and challenges in modern approaches to numerical partial differential equations, (2016), pp. 129–177.

[11] B. COCKBURN, *Hybridizable discontinuous Galerkin methods for second-order elliptic problems: overview, a new result and open problems*, Japan Journal of Industrial and Applied Mathematics, (2023), pp. 1–40.

[12] B. COCKBURN, D. A. DI PIETRO, AND A. ERN, *Bridging the hybrid high-order and hybridizable discontinuous Galerkin methods*, ESAIM: Mathematical Modelling and Numerical Analysis, 50 (2016), pp. 635–650.

[13] B. COCKBURN, B. DONG, AND J. GUZMÁN, *A superconvergent LDG-hybridizable Galerkin method for second-order elliptic problems*, Math. of Comp., 77 (2008), pp. 1887–1916.

[14] B. COCKBURN, J. GOPALAKRISHNAN, AND R. LAZAROV, *Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems*, SIAM Journal on Numerical Analysis, 47 (2009), pp. 1319–1365.

[15] B. COCKBURN AND K. SHI, *Superconvergent HDG methods for linear elasticity with weakly symmetric stresses*, IMA Journal of Numerical Analysis, 33 (2013), pp. 747–770.
[16] D. A. DI PIETRO AND A. ERN, *A hybrid high-order locking-free method for linear elasticity on general meshes*, Comp. Meth. in App. Mechanics and Engineering, 283 (2015), pp. 1-21.

[17] S. DU AND F.-J. SAYAS, *New analytical tools for HDG in elasticity, with applications to elastodynamics*, Mathematics of Computation, 89 (2020), pp. 1745-1782.

[18] M. S. FABIEN, *A GPU-accelerated hybridizable discontinuous Galerkin method for linear elasticity*, Commun Comput Phys, 27 (2020), pp. 513-545.

[19] F. FAUCHER, *hawen: time-harmonic wave modeling and inversion using hybridizable discontinuous Galerkin discretization*, Journal of Open Source Software, 6 (2021).

[20] F. FAUCHER AND O. SCHERZER, *Adjoint-state method for Hybridizable Discontinuous Galerkin discretization, application to the inverse acoustic wave problem*, Computer Methods in Applied Mechanics and Engineering, 372 (2020), p. 113406.

[21] P. FERNANDEZ, A. CHRISTOPHE, S. TERRANA, N. C. NGUYEN, AND J. PERAIRE, *Hybridized discontinuous Galerkin methods for wave propagation*, Journal of Scientific Computing, 77 (2018), pp. 1566-1604.

[22] G. FU, B. COCKBURN, AND H. STOLARSKI, *Analysis of an HDG method for linear elasticity*, International Journal for Numerical Methods in Engineering, 102 (2015), pp. 551-575.

[23] M. GIACOMINI AND R. SEVILLA, *Discontinuous Galerkin approximations in computational mechanics: hybridization, exact geometry and degree adaptivity*, SN Applied Sciences, 1 (2019), p. 1047.

[24] M. GIACOMINI, R. SEVILLA, AND A. HUERTA, *HDGlab: An open-source implementation of the hybridisable discontinuous Galerkin method in matlab*, Archives of Computational Methods in Engineering, 28 (2021), pp. 1941-1986.

[25] D. GIVOLI AND J. B. KELLER, *Non-reflecting boundary conditions for elastic waves*, Wave motion, 12 (1990), pp. 261-279.

[26] J. GOPALAKRISHNAN, S. LANTERI, N. OLIVARES, AND R. PERRUSSEL, *Stabilization in relation to wavenumber in HDG methods*, Advanced Modeling and Simulation in Engineering Sciences, 2 (2015), pp. 1-24.

[27] R. L. HIGDON, *Absorbing boundary conditions for elastic waves*, Geophysics, 56 (1991), pp. 231-241. 10.1190/1.1443035.

[28] A. HUNGRIA, D. PRADA, AND F.-J. SAYAS, *HDG methods for elastodynamics*, Computers & Mathematics with Applications, 74 (2017), pp. 2671-2690.

[29] N. C. NGUYEN AND J. PERAIRE, *Hybridizable discontinuous Galerkin methods for partial differential equations in continuum mechanics*, Journal of Computational Physics, 231 (2012), pp. 5955-5988.

[30] N. C. NGUYEN, J. PERAIRE, AND B. COCKBURN, *High-order implicit hybridizable discontinuous Galerkin methods for acoustics and elastodynamics*, Journal of Computational Physics, 230 (2011), pp. 3695-3718.
[31] H. PHAM, F. FAUCHER, AND H. BARUCQ, *On the implementation of Hybridizable Discontinuous Galerkin discretization for linear anisotropic elastic wave equation: Voigt-notation and stabilization*, Research Report RR-9533, Inria Bordeaux Sud-Ouest; Project-Team Makutu, December 2023. https://hal.archives-ouvertes.fr/hal-04356602/file/RR-9533.pdf.

[32] R. SEVILLA, M. GIACOMINI, A. KARKOULIAS, AND A. HUERTA, *A superconvergent hybridisable discontinuous Galerkin method for linear elasticity*, International Journal for Numerical Methods in Engineering, 116 (2018), pp. 91–116.

[33] N. M. SHAPIRO, M. CAMPILLO, L. MARGERIN, S. K. SINGH, V. KOSTOGLODOV, AND J. PACHECO, *The Energy Partitioning and the Diffusive Character of the Seismic Coda*, Bulletin of the Seismological Society of America, 90 (2000), pp. 655–665.

[34] R. SNIEDER, *Coda wave interferometry and the equilibration of energy in elastic media*, Phys. Rev. E, 66 (2002), p. 046615.

[35] S.-C. SOON, B. COCKBURN, AND H. K. STOLARSKI, *A hybridizable discontinuous Galerkin method for linear elasticity*, International journal for numerical methods in engineering, 80 (2009), pp. 1058–1092.

[36] S. TERRANA, J.-P. VILOTTE, AND L. GUILLOT, *A spectral hybridizable discontinuous Galerkin method for elastic-acoustic wave propagation*, Geophysical Journal International, 213 (2018), pp. 574–602.

[37] L. THOMSEN, *Weak elastic anisotropy*, Geophysics, 51 (1986), pp. 1954–1966.

[38] B. TIE AND A.-S. MOURONVAL, *Systematic development of upwind numerical fluxes for the space discontinuous Galerkin method applied to elastic wave propagation in anisotropic and heterogeneous media with physical interfaces*, Computer Methods in Applied Mechanics and Engineering, 372 (2020), p. 113352.

[39] B. TIE, A.-S. MOURONVAL, V.-D. NGUYEN, L. SERIES, AND D. AUBRY, *A unified variational framework for the space discontinuous Galerkin method for elastic wave propagation in anisotropic and piecewise homogeneous media*, Computer Methods in Applied Mechanics and Engineering, 338 (2018), pp. 299–332.

[40] P.-H. Tournier, P. JOLIVET, V. DOLEAN, H. S. Aghamiry, S. OPERTO, AND S. RIFFO, *3d finite-difference and finite-element frequency-domain wave simulation with multilevel optimized additive schwarz domain-decomposition preconditioner: A tool for full-waveform inversion of sparse node data sets*, Geophysics, 87 (2022), pp. T381–T402.

[41] J. VILA-PÉREZ, M. GIACOMINI, R. SEVILLA, AND A. HUERTA, *Hybridisable discontinuous Galerkin formulation of compressible flows*, Archives of Computational Methods in Engineering, 28 (2021), pp. 753–784.

[42] R. WANG AND R. ZHANG, *A weak Galerkin finite element method for the linear elasticity problem in mixed form*, Journal of Computational Mathematics, 36 (2018), p. 469.

[43] L. C. WILCOX, G. STADLER, C. BURSTEDDE, AND O. GHATTAS, *A high-order discontinuous Galerkin method for wave propagation through coupled elastic-acoustic media*, Journal of Computational Physics, 229 (2010), pp. 9373–9396.
[44] Q. ZHAN, Q. REN, M. ZHUANG, Q. SUN, AND Q. H. LIU, *An exact Riemann solver for wave propagation in arbitrary anisotropic elastic media with fluid coupling*, Computer Methods in Applied Mechanics and Engineering, 329 (2018), pp. 24–39.