# Propagation and scattering of ultrasonic waves in macroscopically anisotropic polycrystalline materials with fiber texture

Juan Camilo Victoria Giraldo, Denis Solas, Jérôme Laurent, Alain Lhémery, Bing Tie

► To cite this version:

Juan Camilo Victoria Giraldo, Denis Solas, Jérôme Laurent, Alain Lhémery, Bing Tie. Propagation and scattering of ultrasonic waves in macroscopically anisotropic polycrystalline materials with fiber texture. Ultrasonics, In press. hal-05521418

HAL Id: hal-05521418

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

Submitted on 24 Feb 2026

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 CC BY 4.0 - Attribution - International License
Propagation and scattering of ultrasonic waves in macroscopically
anisotropic polycrystalline materials with fiber texture

Juan Camilo Victoria Giraldo$^a$, Denis Solas$^b$, Jérôme Laurent$^c$, Alain Lhémery$^c$, Bing Tie$^{a,*}$

$^a$Université Paris-Saclay, CentraleSupélec, ENS Paris-Saclay, CNRS, LMPS, Gif-Sur-Yvette, 91190, France

$^b$Université Paris-Saclay, CNRS, ICMMO, Orsay, 91400, France

$^c$Université Paris-Saclay, CEA List, Palaiseau, 91120, France

Abstract

This work investigates the scattering of ultrasonic elastic waves in polycrystalline materials with macroscopic transverse isotropy (TI). Using the Dyson equation with the First-order smoothing approximation (FOSA) for the mass operator, we obtain the complex wavenumber of the $qL$, $T_1$, and $qT_2$ waves in 3D and $qL$ and $qT_2$ in 2D. The equations are properly modified by including the Green solution for the reference TI homogeneous medium. Additionally, 2D numerical finite element (FE) models are developed, allowing for a direct comparison with the 2D theoretical amplitude attenuation and phase velocities estimations. The numerical simulations are carried out with a set of periodic boundary conditions (PBC), leading to a proper plane wave propagation even in the case of macroscopic anisotropy. The influence of the orientation of the axis of symmetry on the attenuation and phase velocities is discussed. The role of the symmetry-axis orientation is analyzed, and scattering mechanisms are contrasted between 2D and 3D. Numerical and theoretical results show excellent agreement, confirming the increased wave dispersion when propagation occurs along the cross-fiber (or cross-symmetry-axis) direction.

Keywords: ultrasonic scattering, attenuation, phase velocity, elastic wave propagation,

transverse isotropy, macroscopic fiber texture

1. Introduction

Polycrystalline materials are composed of numerous single-crystal grains, each having a unique
shape and anisotropic crystallographic properties. Both the morphology and crystallographic ori-
entation of these grains are generally determined by the specific fabrication procedure used to
produce the material. In such media, there are two possible sources leading to the anisotropy of
ultrasonic scattering phenomena. The first one is the morphological anisotropy induced by the
average grain shape, while the macroscopic elasticity of the polycrystal may remain isotropic, that
is, untextured. The second one is the elastic anisotropy induced by the macroscopic or average
elastic properties. In this work, the macroscopic elastic anisotropy will be assessed. Classical met-
allurgical fabrication processes typically aim to produce macroscopically isotropic polycrystalline
microstructures, in which the grains have random crystallographic orientations. However, in many
practically relevant situations, microstructures exhibit pronounced crystallographic texture, of-
ten associated with elongated grains. Such features may arise during solidification, as observed
in austenitic stainless-steel welds [19, 11], or during thermomechanical processing, such as in tita-
nium alloys, where macrozones sharing a common crystallographic orientation can develop [33, 32].

The presence of macrozones does not necessarily result in macroscopic elastic anisotropy; however,
such anisotropy may emerge when macrozones reach sufficiently large dimensions. Comparable

<footnote>*Corresponding author
Email address: bing.tie@centralesupelec.fr (Bing Tie)</footnote>
microstructural characteristics resulting in macroscopic elastic anisotropy have also been studied and modelled by Kong et al. [17] in the context of recently developed Wire and Laser Additive Manufacturing (WLAM) process. More specifically, in WLAM components, all grains may share a common aligned axis – forming a distribution of crystallographic orientations known as fiber texture – where the fiber direction corresponds to the favored aligned axis. In polycrystalline materials, elastic wave propagation is scattered at the grain boundaries due to the stiffness variations. This scattering leads to wave attenuation and induces phase shifts during propagation. Unlike in the macroscopically isotropic case, wave propagation in a fiber-textured microstructure strongly depends on the propagation direction relative to the fiber axis. Appropriate modeling and accurate prediction of elastic wave attenuation and phase velocity in such specific media are highly useful for Non-Destructive Testing (NDT) and material characterization [4, 3], with the use of inversion algorithms such as the ones proposed by Allaire et al. [2].

The problem of elastic wave scattering in macroscopically isotropic media has been widely examined. To date, the two main contributions are those of Stanke and Kino [24], building on the work of Karal and Keller [16], and Weaver [30]. Both approaches are comparable and lead to the same results under the Born approximation. Their derivations were initially presented for a simplified isotropic polycrystalline microstructure with equiaxed grains. Their development, however, is general and can be extended to more general hypotheses. Several works have extended these base models to consider the grains to have an elongated shape [31, 14, 6, 5, 15].

Nevertheless, the elastic wave scattering in textured polycrystalline microstructures has not been extensively studied. The first works [12, 13, 1] to treat textured polycrystalline microstructures with equiaxed grains considered the background reference medium to be isotropic. As would later be demonstrated by Turner [27], the textured microstructure has, on average, a macroscopic anisotropy that strictly requires the reference homogeneous medium to be also anisotropic. Turner [27] was the first to consider a transversely isotropic (TI) reference medium by introducing the anisotropic Green's dyadic and the analytical expressions for the single- and two-point elastic correlation tensors. More recently, Li and Rokhlin [18] proposed a general approach to treat any type of macroscopic anisotropy in addition to the geometric anisotropy of elongated grains. Their work is based on the expansion of the anisotropic Green's dyadic in terms of the polarization eigenmodes obtained from Christoffel's equation. As they clearly stated, the limitation lies in the increased complexity of the single- and two-point elastic correlation tensors. For the most general case, the solution is only possible numerically.

In this work, we follow up on Turner's [27] work to obtain 3D and 2D analytical solutions for the elastic wave attenuation and phase velocity in an Inconel fiber-textured microstructure with equiaxed grains. Furthermore, 2D numerical finite element (FE) simulations with a space-discontinuous Galerkin framework are performed to validate the 2D theoretical results. These numerical simulations are carried out with a set of periodic boundary conditions (PBC), which leads to proper plane wave propagation even in the case of elastic anisotropy. The numerical-theoretical comparison and the study of the dimensionality and scattering mechanisms will contribute to the understanding of elastic wave scattering in such complex and specific media as the TI one.

The outline of the paper is as follows. In Section 2, we present the first-order smooth approximation (FOSA), Turner's [27] proposed anisotropic Green's dyadic, and our 2D generalization. Then, Section 3 presents the FE models and the methodology followed to estimate the wave attenuation and phase velocity. Finally, in Section 4, the directional and frequency dependencies of the wave scattering is analyzed, and the theoretical and numerical results with special interest in the dimensionality analysis between the 2D and 3D elastic scattering behavior are presented and discussed.
## 2. Theoretical expressions for attenuation and phase velocity for the fiber texture case

### 2.1. The elastodynamic wave equation

Let us consider the elastic wave propagation in a heterogeneous medium with a spatially varying elastic stiffness tensor $C(x)$, decomposed into the following two parts,

$$C(x) = C^0 + \delta C(x), \qquad (1)$$

with $C^0$ the elastic stiffness tensor of the homogeneous reference medium and $\delta C$ the spatially varying fluctuation in $C$.

In the frequency domain, the wave equation reads as:

$$\mathcal{L}u = (\nabla_x \cdot (C : \nabla_x) + \rho\omega^2 I) u = f. \qquad (2)$$

Or with indicial notations:

$$\frac{\partial}{\partial x_j} \left( C_{ijkl} \frac{\partial u_k}{\partial x_l} \right) + \rho \omega^2 u_i = f_i,$$

with the Einstein summation convention.

By decomposing the linear differential operator $\mathcal{L}$ into to operators $\mathcal{L}^0$ and $\mathcal{L}^1(x)$ according to the decomposition Eq. (1),

$$\mathcal{L}^0 = \nabla_x \cdot (C^0 : \nabla_x) + \rho\omega^2 I, \quad \text{with} \quad \mathcal{L}_{ik}^0 = C_{ijkl}^0 \frac{\partial^2}{\partial x_j \partial x_l} + \rho\omega^2 \delta_{ik}, \qquad (3a)$$

$$\mathcal{L}^1(x) = \nabla_x \cdot (\delta C(x) : \nabla_x), \quad \text{with} \quad \mathcal{L}_{ik}^1(x) = \frac{\partial}{\partial x_j} \left( \delta C_{ijkl}(x) \frac{\partial}{\partial x_l} \right), \qquad (3b)$$

where $\delta_{ik}$ denotes the Kronecker delta representing an identity operator. Then, Eq. (2) can be re-written

$$(\mathcal{L}^0 + \mathcal{L}^1(x)) u(x, \omega) = f(x, \omega). \qquad (4)$$

In what follows, the dependency on $\omega$ will be omitted for simplicity.

### 2.2. The Dyadic Green's function

In the following, the dyadic Green's function is developed as previously presented by Weaver [30], Frisch [10], Roy and Kube [22]. Let denote by $g_m(x, x')$, the solution of the wave equation Eq. (2) with $f = \delta(x - x')e_m$, a unit 3D-impulse along the $e_m$ direction, that is,

$$\mathcal{L}g_m(x, x') = \delta(x - x')e_m, \qquad (5)$$

or with indicial notations:

$$\frac{\partial}{\partial x_j} \left( C_{ijkl} \frac{\partial (g_m)_k}{\partial x_l} \right) + \rho \omega^2 \delta_{ik} (g_m)_k = \delta(x - x') \delta_{im}.$$

Denoting by $G = g_m \otimes e_m$, the dyadic Green's function, then it is the solution of the following Green equation, a tensorial equation:

$$\mathcal{L}G(x, x') = \delta(x - x')I, \qquad (6)$$

or with indicial notations:

$$\mathcal{L}_{ik} G_{km}(x, x') = \delta(x - x') \delta_{im}.$$

Due to the linearity, the solution of Eq. (2) can be expressed as:

$$u(x) = u^H(x) + \int G(x, x') \cdot f(x') dx', \qquad (7a)$$
$$u_i = u_i^H + \int G_{im}(x, x')f_m(x')dx' \left( = \int (g_m)_i f_m dx' \right), \quad (7b)$$

where $u^H$ is the solution of the homogeneous equation $\mathcal{L}u^H = 0$.

For the homogeneous reference medium, the dyadic Green's function $G^0$ is the solution of the following governing equation:

$$\mathcal{L}^0 G^0(x, x') = \delta(x - x')I, \quad (8)$$

We note that assuming statistical homogeneity of the medium implies that the variation of the stiffness tensor is a zero-mean random function, with $\langle \delta C(x) \rangle = 0$. The reference medium is therefore homogeneous, and in this case, we have $G^0(x, x') = G^0(x - x')$.

Combining Eqs. (6) and (8), the following equation is straightforward:

$$\mathcal{L}^0 G(x, x') = \mathcal{L}^0 G^0(x, x') - \mathcal{L}^1(x) G(x, x'), \quad (9)$$

then, applying Eq. (7a) results in,

$$G(x, x') = G^0(x, x') + u^{0H}(x, x') - \int G^0(x, x_1) \mathcal{L}^1(x_1) G(x_1, x') dx_1, \quad (10)$$

where $u^{0H}(x, x')$ is the solution of the homogeneous equation $\mathcal{L}^0 u^{0H} = 0$, which is a plane wave propagating in the homogeneous medium. As the dyadic Green's function $G(x, x')$ is the solution of a 3D-impulse loading, so, it must verify $\lim_{|x-x'| \to +\infty} G(x, x') = 0$, this results in $u^{0H} = 0$.

Finally, we get the following Lippmann-Schwinger-type equation, which expresses the solution, the heterogeneous Green dyadic $G(x, x')$, in a recursive form:

$$G(x, x') = G^0(x, x') - \int G^0(x, x_1) \mathcal{L}^1(x_1) G(x_1, x') dx_1. \quad (11)$$

Iterating the right-hand side of equation (11) yields an infinite Neumann series, known as Born series,

$$G(x, x') = G^0(x - x') + G^{(1)}(x, x') + \dots + G^{(p)}(x, x') + \dots, \quad (12)$$

where the terms represent increasing orders of perturbation to the homogeneous solution and each term $G^{(p)}$ involves $\delta C(x_1) \otimes \dots \otimes \delta C(x_p)$, the p-order p-point correlation function of $\delta C$ (see Appendix A1: The Born series).

A mathematical representation of a polycrystalline microstructure is achievable in a statistical sense, which is why the solution Equation (11) is sought in terms of the average Green's function. The elastodynamic Dyson equation for the mean Green's dyadic is obtained using the diagrammatic method [8]:

$$\langle G(x, x') \rangle = G^0(x, x') + \iint G^0(x, x_1) \cdot m(x_1, x_\infty) \left( \langle G(x_\infty, x') \rangle \right) d^3x_1 d^3x_\infty, \quad (13)$$

where $m$ is known as the "mass" or "self-energy" operator. This requires ensemble-averaged statistics of the medium, which, together with statistical assumptions, can lead to feasible solutions [21, 22]. The critical assumption of statistical homogeneity of the medium allows the even-order terms of the resulting equation to be expressed solely in terms of two-point correlation averages [8, 21, 22]. Furthermore, when the variation of the stiffness tensor is assumed to follow a multivariate Gaussian distribution, all odd-order correlation averages vanish.

Weaver [30] originally considered the first term in the mass operator, after application of what is called the first-order smooth approximation (FOSA). The FOSA of the mass operator can be written as [30, 21, 22]:

$$m^1(x_1, x_2)(\bullet) = -\nabla_{x_1} \cdot \left( \langle \delta C(x_1) \otimes \delta C(x_2) \rangle : \left( \nabla^2_{x_1 x_2} G^0(x_1, x_2) \otimes \nabla_{x_2}(\bullet) \right) \right). \quad (14)$$
The FOSA accounts for multiple scattering up to a certain degree. As the series is infinite, the first terms are expected to provide the dominant contribution to multiple scattering. However, FOSA may be limited when applied to highly anisotropic materials. Readers are referred to the work of Roy [21], Roy and Kube [22] for further discussion on the contribution of higher-order terms in the mass operator.

## 2.3. Solution in the wavenumber domain

The Dyson equation (13) is more easily solved in the frequency-wavenumber domain. It should be noted that the transformation to the wavenumber domain contributes to simplifying the equation by expressing the Green's function in the wavenumber domain in a simpler and more elegant form than its spatial counterpart. The spatial derivatives are also treated by the transform, avoiding any increase in the necessary terms. Again, thanks to the assumption of statistical homogeneity, the Dyson equation – which can be interpreted as a double convolution with respect to $(x - x')$ – can be written in the frequency-wavenumber domain as [30, 21, 22]:

$$ \langle \mathbf{G}(\mathbf{p}) \rangle = \mathbf{G}^0(\mathbf{p}) - \mathbf{G}^0(\mathbf{p})\boldsymbol{\sigma}(\mathbf{p})\langle \mathbf{G}(\mathbf{p}) \rangle, \quad (15) $$

where $\boldsymbol{\sigma}$, the spatial transform of the self-energy, is defined as:

$$ -\boldsymbol{\sigma}(\mathbf{p})\delta^3(\mathbf{p} - \mathbf{q}) = \frac{1}{(2\pi)^3} \int \int \mathbf{m}(\mathbf{x}, \mathbf{x}') \exp(-i(\mathbf{p} \cdot \mathbf{x} - \mathbf{q} \cdot \mathbf{x}')) d^3\mathbf{x}d^3\mathbf{x}'. \quad (16) $$

Equation (15) is simplified by grouping the averaged Green tensor and dividing by the homogeneous Green tensor:

$$ \langle \mathbf{G}(\mathbf{p}) \rangle = \left( (\mathbf{G}^0(\mathbf{p}))^{-1} + \boldsymbol{\sigma}(\mathbf{p}) \right)^{-1}. \quad (17) $$

## 2.4. Dispersion equations for a macroscopically anisotropic medium

In the following, we follow up on Turner [27] work to obtain 3D and 2D analytical solutions for the elastic wave attenuation and phase velocity in a fiber-textured microstructure with equiaxed grains. Considering a polycrystalline medium with texture leads to a macroscopically anisotropic medium. Depending on the texture different levels of anisotropy can be obtained. Therefore, a general homogeneous anisotropic Green function is necessary to solve the Dyson equation (17).

Recalling the elastodynamic Green function equation in the frequency-wavenumber domain for a homogeneous medium is normally expressed as:

$$ (\rho\omega^2\mathbf{I} - \mathbf{\Gamma}) \mathbf{G}^0(\mathbf{p}) = \mathbf{I}, \quad (18) $$

where $\mathbf{\Gamma} = \mathbf{p}\mathbf{C}^0\mathbf{p}$ is the Christoffel's tensor for an anisotropic medium. Here we keep the subscript 0 to denote the homogeneous stiffness tensor, however it can be either isotropic or anisotropic. Solving for the Green tensor:

$$ \mathbf{G}^0(\mathbf{p}) = (\rho\omega^2\mathbf{I} - \mathbf{\Gamma})^{-1}. \quad (19) $$

As explained in Turner [27], Li and Rokhlin [18], the solution can be expanded in terms of the orthonormal eigenvectors of $\mathbf{\Gamma}$ as:

$$ \mathbf{G}^0(\mathbf{p}) = \sum_{M=1}^{3} \frac{\hat{\mathbf{u}}_M(\hat{\mathbf{p}}) \otimes \hat{\mathbf{u}}_M(\hat{\mathbf{p}})}{\rho(\omega^2 - p^2 V_{0M}^2(\hat{\mathbf{p}}))} = \sum_{M=1}^{3} \frac{\hat{\mathbf{u}}_M(\hat{\mathbf{p}}) \otimes \hat{\mathbf{u}}_M(\hat{\mathbf{p}})}{\rho V_{0M}^2(\hat{\mathbf{p}})(p_{0M}^2 - p^2)} = \sum_{M=1}^{3} g_{0M}(\mathbf{p}) \hat{\mathbf{u}}_M(\hat{\mathbf{p}}) \otimes \hat{\mathbf{u}}_M(\hat{\mathbf{p}}), \quad (20) $$

where $\mathbf{p} = p\hat{\mathbf{p}}$, with $\hat{\mathbf{p}}$ the unit wavevector, $V_{0M}(\hat{\mathbf{p}})$ and $\hat{\mathbf{u}}_M(\hat{\mathbf{p}})$ are respectively eigenvalues and eigenvectors of the Christoffel tensor – that is, the wave velocities and polarization vectors of the propagating modes, in physical terms. It is clear, in contrast with the isotropic case, that both the wave velocities and the polarization vectors now depend on the wave propagation direction. Indeed, the orthonormal basis formed by the polarization vector is not aligned with the orthonormal basis formed by the propagation direction. In general, any level of anisotropy could be considered by solving the eigenvalue problem numerically for any anisotropy of $\mathbf{C}^0$.
In the case of transverse isotropy, which remains common in practice and is also called “fiber texture”, the crystals share a common aligned axis that serves as a symmetry axis of the medium. In the following, this axis is defined by a unit vector $\hat{n}$. Turner [27] dealt with this case and defined an orthonormal basis $(\hat{u}_M)_{M=T_1,qL,qT_2}$, in terms of $\hat{u}_{T_1} \equiv \hat{u}_1$, the direction vector of the transverse horizontal mode ($T_1$), which is perpendicular to both $\hat{n}$ and $\hat{p}$; $\hat{u}_{qL} \equiv \hat{u}_2$, the direction vector of the quasi-longitudinal mode ($qL$); and $\hat{u}_{qT_2} \equiv \hat{u}_3$, the direction vector of the quasi-transverse mode ($qT_2$). The latter two lie in the plane defined by $\hat{n}$ and $\hat{p}$. By denoting $\psi$ as the angle between $\hat{u}_{qL}$ and $\hat{p}$, and $\theta$ as the angle between $\hat{p}$ and $\hat{n}$, we have:

$$ \hat{u}_{T_1} = \frac{\hat{p} \times \hat{n}}{|\hat{p} \times \hat{n}|}, \qquad (21) $$

$$ \hat{p} = \hat{u}_{qL} \cos \psi + \hat{u}_{qT_2} \sin \psi, \qquad (22) $$

$$ \hat{n} = \hat{u}_{qL} \cos(\psi + \theta) + \hat{u}_{qT_2} \sin(\psi + \theta), \qquad (23) $$

Turner [27] developed relations that allow the calculation of the angle $\psi$ and the wave velocities $V_{0T_1}$, $V_{0qL}$, and $V_{0qT_2}$ from the five independent constants of a transverse isotropic stiffness tensor.

As for the mass operator, it must follow Green's symmetry and therefore can be expanded in terms of the polarization vectors [27, 18]:

$$ \sigma(\mathbf{p}) = \sum_{M=T_1,qL,qT_2} \sigma^M(\mathbf{p}) \hat{u}_M \otimes \hat{u}_M. \qquad (24) $$

Setting the inverse of (17) to zero and knowing the form of $G^0$ in (20), the following dispersion relations can be obtained: for $M = T_1, qL$, and $qT_2$,

$$ p_{0M}^2(\hat{p}) - p_M^2(\hat{p}) + p_M^2 \bar{\sigma}^M(\mathbf{p}_M) = 0, \qquad (25) $$

with $p_M = p_M \hat{p}$ and $\bar{\sigma}^M = \sigma^M / (\rho p_M^2 V_{0M}^2)$. It is important to note that the dependence on $\hat{p}$, the wave propagation direction with respect to the symmetry axis $\hat{n}$, in terms of $\theta$, is implicitly expressed.

As far as $\langle \delta C(x_1) \otimes \delta C(x_2) \rangle$, the two-point correlation function of $\delta C(x)$ in Equation (14), is concerned, and under the assumption of statistical homogeneity in the monophase case, it – along with its spatial Fourier transform – can be decomposed into decoupled tensorial and spatial parts as follows:

$$ \langle \delta C(x_1) \otimes \delta C(x_2) \rangle = \Xi W(x_1 - x_2), \quad \Lambda(p) = \Xi \tilde{W}(p) \qquad (26) $$

where $\Xi = \langle \delta C \otimes \delta C \rangle$, with $\delta C$ the single crystal constants, $W(x_1 - x_2)$ and $\tilde{W}(p)$ are respectively the spatial two-point correlation (TPC) function and its Fourier transform.

As previously mentioned in 2.4, the macroscopic transverse isotropy is governed by an axis of symmetry $\hat{n}$, and in this case, $\Xi$, the eighth-rank covariance, is also a function of $\hat{n}$. The expressions for the Voigt average tensor and covariance tensor in the case of cubic symmetry are given by Turner [27].

Using Equation (20), the decomposition of the homogeneous Green's function into different wave components, the mass operators can be further decoupled as: for $M = T_1, qL, qT_2$,

$$ \bar{\sigma}^M(\mathbf{p}_M) = \sum_{Q=T_1,qL,qT_2} \bar{\sigma}^{M \to Q}(\mathbf{p}_M). \qquad (27) $$

The expressions for the decoupled mass operators can be found to be [21, 22]:

$$ \bar{\sigma}^{M \to Q}(\mathbf{p}_M) = -\frac{1}{\rho V_{0M}^2(\hat{p})} \int IP_{M \to Q}(\hat{p}, \hat{s}) g_{0Q}(s_Q) p_Q^2 \tilde{W}(\mathbf{p}_M - s_Q) d^3 s_Q, \qquad (28) $$

where $s_Q = p_Q \hat{s}$, with $\hat{s}$ the unit wave vector of the Q-type scattered wave, and $IP_{M \to Q}$ are the inner products of the covariance stiffness tensor $\Xi$ with the unit wave vectors $\hat{p}$ and $\hat{s}$ of the incident and scattered waves, as well as their respective polarization vectors (cf. Equation (32)). Equation (28) is typically formulated for the 3D case, but it remains valid in 2D when the appropriate dimension-specific inner products and spatial TPC function are used.
## 2.5. 3D and 2D solution of the dispersion equations: attenuation and phase velocity

The dispersion equations (25) are non-linear equations that could be solved iteratively. However, here we make use of the Born approximation, which leads to linear expressions of the attenuation and wave phase velocity. As has been discussed in other works [30, 8, 24], the Born approximation limits the validity of the frequency regimes up to the stochastic. Calvet and Margerin [6] proposes to solve the dispersion equations with a spectral method, which allows obtaining results up to the geometric region. Roy's work [21, 22] also contains different numerical strategies to compute the attenuation and phase velocity iteratively up to the geometric region. These approaches are not implemented here, but are an interesting perspective for future work.

The Born approximation allows to simplify the dispersion equations (25) as: for $M = qT_1, qL, qT_2$,

$$p_M(\hat{\boldsymbol{p}}) = p_{0M}(\hat{\boldsymbol{p}}) - \frac{p_{0M}(\hat{\boldsymbol{p}})}{2} \bar{\sigma}^M(\boldsymbol{p}_M). \qquad (29)$$

Since the attenuation is obtained as the imaginary part of the complex wavenumber $p_M$, and the phase velocity from the real part, then one can find that:

$$\alpha_M(\boldsymbol{p}_M) = -\frac{p_{0M}(\hat{\boldsymbol{p}})}{2} \operatorname{Im}(\bar{\sigma}^M(\boldsymbol{p}_M)), \qquad (30)$$

$$V_M(\boldsymbol{p}_M) = V_{0M}(\hat{\boldsymbol{p}}) (1 - \operatorname{Re}(\bar{\sigma}^M(\boldsymbol{p}_M)))^{-1}.$$

Equation (30) can be further decomposed as follows:

$$\alpha_M(\boldsymbol{p}_M) = \sum_{Q=T_1, qL, qT_2} \alpha_{M \to Q}(\boldsymbol{p}_M), \qquad (31)$$
$$V_M(\boldsymbol{p}_M) = \sum_{Q=T_1, qL, qT_2} V_{M \to Q}(\boldsymbol{p}_M).$$

To derive the expressions for attenuation and phase velocity, it is necessary to consider both the inner products and the spatial TPC function in the 3D and 2D cases, with the 2D case representing a novel contribution of the present work.

The vectors appearing in the inner products arise from the spatial derivatives of the functions involved in the mass operator, following the application of the spatial Fourier transform. Therefore, the inner product has the following general form:

$$\begin{aligned} IP_{M \to Q} &= (\hat{\boldsymbol{p}} \otimes \hat{\boldsymbol{u}}_M(\hat{\boldsymbol{p}}) \otimes \hat{\boldsymbol{u}}_Q(\hat{\boldsymbol{s}}) \otimes \hat{\boldsymbol{s}}) : \Xi : (\hat{\boldsymbol{s}} \otimes \hat{\boldsymbol{u}}_Q(\hat{\boldsymbol{s}}) \otimes \hat{\boldsymbol{u}}_M(\hat{\boldsymbol{p}}) \otimes \hat{\boldsymbol{p}}) \\ &= \sum_{ijkl} \Xi_{ijkl}^{\alpha\beta\gamma\delta} \hat{p}_{\alpha} \hat{u}_{M\beta} \hat{u}_{Q\gamma} \hat{s}_{\delta} \hat{s}_{\epsilon} \hat{u}_{Qj} \hat{u}_{Mk} \hat{p}_{l}, \end{aligned} \qquad (32)$$

In the case of general anisotropy, the parametrization of the polarization vectors is complex, whereas it becomes simpler for the transverse isotropic medium – for example, when using the orthonormal basis shown in Equations (21). The expression can be further simplified by choosing the symmetry axis to be aligned with the $e_3$ axis and constraining the propagating direction $\hat{\boldsymbol{p}}$ to lie within the $e_2-e_3$ plane. This simplification is valid in the case studied in the present work, which concerns polycrystals with equiaxed grains and macroscopic transverse isotropy. $\hat{n}$, $\hat{\boldsymbol{p}}$, and $\hat{s}$, the unit vectors of the symmetry axis, as well as those of the incident and scattered wave vectors (see Figure 1), can therefore be written as:

$$\begin{aligned} \hat{n} &= e_3, \\ \hat{\boldsymbol{p}} &= \sin \theta_p e_2 + \cos \theta_p e_3, \\ \hat{s} &= \sin \theta_s (\cos \varphi_s e_1 + \sin \varphi_s e_2) + \cos \theta_s e_3. \end{aligned} \qquad (33)$$

Then, the polarization vectors for the incident and scattered waves can be written as [27]:

$$\begin{aligned} \hat{u}_{T_1}(\hat{\boldsymbol{p}}) &= e_1, & \hat{u}_{T_1}(\hat{s}) &= \sin \varphi_s e_1 - \cos \varphi_s e_2, \\ \hat{u}_{qL}(\hat{\boldsymbol{p}}) &= \sin \gamma_p e_2 + \cos \gamma_p e_3, & \hat{u}_{qL}(\hat{s}) &= \sin \gamma_s (\cos \varphi_s e_1 + \sin \varphi_s e_2) + \cos \gamma_s e_3, \qquad (34) \\ \hat{u}_{qT_2}(\hat{\boldsymbol{p}}) &= -\cos \gamma_p e_2 + \sin \gamma_p e_3, & \hat{u}_{qT_2}(\hat{s}) &= -\cos \gamma_s (\cos \varphi_s e_1 + \sin \varphi_s e_2) + \sin \gamma_s e_3, \end{aligned}$$
where $\gamma = \theta + \psi(\theta)$. Therefore, the inner products depend on the angle $\theta$ between the symmetry axis and the propagation direction as well as the angle $\psi$ between the propagating direction and the polarization vector. It can be observed that the left group of equations is obtained due to the condition $\varphi_p = 90^{\circ}$.

<figure><img src="image_2.png" /><figcaption>Figure 1. Illustration of the incident wave $\hat{p}$ polarization basis dependent on the angle $\theta_p$ with respect to the fiber axis $\hat{n}$.</figcaption></figure>

Finally, the required inner products $IP_{M \to Q} = \nu^2 I_{M \to Q}$, with $\nu = c_{11} - c_{12} - 2c_{44}$ for single cubic crystal anisotropy, can be found to be [27]:

$$ 
\begin{align*}
I_{T_1 \to T_1} &= \frac{\sin^2 \theta_s \sin^2 \theta_p}{32}, \\
I_{T_1 \to qL} &= \frac{\sin^2 \theta_s \sin^2 \gamma_s \sin^2 \theta_p}{32}, & I_{T_1 \to qT_2} &= \frac{\sin^2 \theta_s \cos^2 \gamma_s \sin^2 \theta_p}{32}, \\
I_{qL \to T_1} &= \frac{\sin^2 \theta_s \sin^2 \gamma_p \sin^2 \theta_p}{32}, & I_{qT_2 \to T_1} &= \frac{\sin^2 \theta_s \cos^2 \gamma_p \sin^2 \theta_p}{32}, \tag{35} \\
I_{qL \to qL} &= \frac{\sin^2 \theta_s \sin^2 \gamma_s \sin^2 \gamma_p \sin^2 \theta_p}{32}, & I_{qL \to qT_2} &= \frac{\sin^2 \theta_s \cos^2 \gamma_s \sin^2 \gamma_p \sin^2 \theta_p}{32}, \\
I_{qT_2 \to qL} &= \frac{\sin^2 \theta_s \sin^2 \gamma_s \cos^2 \gamma_p \sin^2 \theta_p}{32}, & I_{qT_2 \to qT_2} &= \frac{\sin^2 \theta_s \cos^2 \gamma_s \cos^2 \gamma_p \sin^2 \theta_p}{32}.
\end{align*}
$$

Since, in the case of transverse isotropy, the parametrization is carried out in a frame defined by the plane spanned by the symmetry axis and the incident wave vector, along with its perpendicular direction, the problem decouples into two parts: a 2D problem within the plane, and a pure transverse mode polarized in the perpendicular direction. Therefore, in the 2D case, the polarization vector of the scattered waves lies within the plane spanned by the symmetry axis and the incident wave vector. It can be shown that the inner products in this 2D configuration are identical to those in the corresponding 3D case for the $qL$ and $qT_2$ modes (cf. the terms in the last two lines of Equation (35)).

As for the spatial TPC function, its general 3D form for a grain of ellipsoidal shape is commonly expressed as an exponential function. Nevertheless, it is known that in real microstructures, the grain shape and their size distribution does not follow a perfect exponential function and numerical/experimental methods are utilized to approximate more realistic functions [9, 23, 28]. In this section, however, we keep the idealized exponential function for the following theoretical development. The exponential function is written as follows [31]:

$$ W(\boldsymbol{r}) = \exp\left(-\sqrt{\frac{r_1^2}{a_1^2} + \frac{r_2^2}{a_2^2} + \frac{r_3^2}{a_3^2}}\right), \qquad (36) $$
242 and its spatial Fourier transform reads as:

$$ \tilde{W}(\mathbf{q}) = \frac{1}{(2\pi)^3} \int W(\mathbf{r}) \exp(-i\mathbf{q} \cdot \mathbf{r}) d^3\mathbf{r} \equiv \bar{W}^{3D}(\mathbf{q})/(2\pi)^3, \qquad (37) $$

243 where:

$$ \bar{W}^{3D}(\mathbf{q}) = \frac{8\pi a_1 a_2 a_3}{(1 + (\mathbf{q} \cdot \mathbf{q})(\mathbf{a} \cdot \mathbf{a}))^2}. \qquad (38) $$

244 It is noteworthy that, by giving the expression of the spatial TPC function like Equation (36), it is assumed that the principal axes of the elongated grain is aligned with ($e_1, e_2, e_3$).

246 For the 2D case, defined within the $e_2-e_3$ plane as previously parameterized, we have, by analogy, $\tilde{W}(\mathbf{q}) = \bar{W}^{2D}(\mathbf{q})/(2\pi)^2$, where:

247

$$ \bar{W}^{2D}(\mathbf{q}) = \frac{2\pi a_2 a_3}{(1 + (\mathbf{q} \cdot \mathbf{q})(\mathbf{a} \cdot \mathbf{a}))^{3/2}}. \qquad (39) $$

248 In the following, analytical formulas for attenuation and phase velocity are provided for the case of elongated grains. However, as it has been mentioned before, in this work, we will focus on the fiber texture effect on the elastic wave scattering, while keeping the assumption of grains with equiaxed shape. With this assumption, the case studied in the present work reduces to having $a_1 = a_2 = a_3$.

250

251 With this assumption, the case studied in the present work reduces to having $a_1 = a_2 = a_3$. Studying the combined effect of morphological and macroscopic elastic anisotropies on elastic wave scattering is left for future work.

253

254 We can now derive formulas for attenuation and phase velocity.

255 2.5.1. Attenuation

256 To calculate the attenuation, as shown in equation (30), one needs to compute the imaginary part of the mass operator. Using the Born approximation, which replaces the heterogeneous wavenumber with the homogeneous one, makes the homogeneous Green function (20) unstable. Huang et al. [14] utilizes a normalization that helps to avoid this instability, especially for the numerical integration routines. The normalization consists in defining a new variable $\xi = p_M/p_{0M}$, which allows to rewrite the imaginary part of the homogeneous Green function as:

258

259

260

261

$$ \mathrm{Im}(g_{0M}) = \frac{\pi}{2\rho V_{0M}^2(\hat{p})p_{0M}^2} \delta(1 - \xi) \qquad (40) $$

262 With the same normalization, the Fourier transform of the TPC spatial function can also be rewritten as [14]:

263

$$ \bar{W}^{3D}(\mathbf{q}) = \frac{8\pi a_1 a_2 a_3}{(1 + C(\xi))^2}, \quad \text{and,} \quad \bar{W}^{2D}(\mathbf{q}) = \frac{2\pi a_1 a_2}{(1 + C(\xi))^{3/2}}, \qquad (41) $$

264 where $C(\xi) = C_0 - \xi C_1 + \xi^2 C_2$, with $C_0, C_1, C_2$ the coefficients depending on the incident and scattered wave basis since we have that $\mathbf{q} = p_{0M}\hat{\mathbf{p}} - p_{0Q}\hat{\mathbf{s}}$ and can be written as follows, considering the form of the incident unit wave vector $\hat{\mathbf{p}}$ and scattered unit wave vector $\hat{\mathbf{s}}$ in Eq. (33): in 3D, we have,

265

266

267

$$ \begin{aligned} C_0^{3D} &= p_{0M}^2 (a_2^2 \sin^2 \theta_p + a_3^2 \cos^2 \theta_p), \\ C_1^{3D} &= 2p_{0M}p_{0Q} (a_2^2 \sin \varphi_s \sin \theta_p \sin \theta_s + a_3^2 \cos \theta_p \cos \theta_s), \\ C_2^{3D} &= p_{0Q}^2 (a_1^2 \cos^2 \varphi_s + a_2^2 \sin^2 \varphi_s) \sin^2 \theta_s + a_3^2 \cos^2 \theta_s), \end{aligned} \qquad (42) $$

268 and in 2D, we have:

$$ \begin{aligned} C_0^{2D} &= p_{0M}^2 (a_2^2 \sin^2 \theta_p + a_3^2 \cos^2 \theta_p), \\ C_1^{2D} &= 2p_{0M}p_{0Q} (a_2^2 \sin \varphi_s \sin \theta_p \sin \theta_s + a_3^2 \cos \theta_p \cos \theta_s), \\ C_2^{2D} &= p_{0Q}^2 (a_1^2 \cos^2 \varphi_s + a_2^2 \sin^2 \varphi_s) \sin^2 \theta_s + a_3^2 \cos^2 \theta_s). \end{aligned} \qquad (43) $$
269 Substituting in the mass operator (28), performing integration in the spherical domain, and being
270 consistent with the Born approximation by taking $p_M = p_{0M}$ and $p_Q = p_{0Q}$, one has:

$$
\mathrm{Im}(\bar{\sigma}^{M \to Q}(\mathbf{p}_M)) = - \frac{a_1 a_2 a_3 \nu^2 \omega^3}{2\pi \rho^2 V_{0M}^2(\theta_p)} \int_0^{2\pi} \int_0^{\pi} \frac{I_{M \to Q}(\theta_s)}{V_{0Q}^5(\theta_s) (1+C(\xi))^2} \sin \theta_s d\theta_s d\varphi_s \int_0^{\infty} \xi^4 \delta(1-\xi) d\xi. \quad (44)
$$

271 The delta function drops the integral over $\xi$, leading to:

$$
\mathrm{Im}(\bar{\sigma}^{M \to Q}(\mathbf{p}_M)) = - \frac{a_1 a_2 a_3 \nu^2 \omega^3}{2\pi \rho^2 V_{0M}^2(\theta_p)} \int_0^{2\pi} \int_0^{\pi} \frac{I_{M \to Q}(\theta_s)}{V_{0Q}^5(\theta_s) (1+C(1))^2} \sin \theta_s d\theta_s d\varphi_s. \quad (45)
$$

272 By analogy, the 2D mass operator reduces to:

$$
\mathrm{Im}(\bar{\sigma}^{M \to Q}(\mathbf{p}_M)) = - \frac{a_1 a_2 \nu^2 \omega^2}{4 \rho^2 V_{0M}^2(\theta_p)} \int_0^{2\pi} \frac{I_{M \to Q}^{2D}(\theta_s)}{V_{0Q}^4(\theta_s) (1 + C(1))^{3/2}} d\theta_s. \quad (46)
$$

Finally, the equations for the attenuation can be written as:

$$
\alpha_{M \to Q}^{3D} = \frac{a_1 a_2 a_3 \nu^2 \omega^4}{4\pi \rho^2 V_{0M}^3(\theta_p)} \int_0^{2\pi} \int_0^{\pi} \frac{I_{M \to Q}(\theta_s)}{V_{0Q}^5(\theta_s) (1+C(1))^2} \sin\theta_s d\theta_s d\phi_s, \quad (47)
$$

$$
\alpha_{M \to Q}^{2D} = \frac{a_1 a_2 \nu^2 \omega^3}{8 \rho^2 V_{0M}^3(\theta_p)} \int_0^{2\pi} \frac{I_{M \to Q}^{2D}(\theta_s)}{V_{0Q}^4(\theta_s) (1 + C(1))^{3/2}} d\theta_s. \quad (48)
$$

273 2.5.2. Phase Velocity

274 The calculation of the phase velocity is only a little trickier because of the Principal Value coming
275 from the homogeneous Green function. Using the same normalization mentioned above, that is
276 $\xi = p_M/p_{0M}$, the real part of the homogeneous Green function can be rewritten as:

$$
\mathrm{Re}(g_{0M}) = \frac{1}{\rho V_{0M}^2(\hat{\mathbf{p}}) p_{0M}^2} \mathrm{PV} \left( \frac{1}{1-\xi^2} \right). \qquad (49)
$$

277 Replacing in the mass operator (28), leads to:

$$
\mathrm{Re}(\bar{\sigma}^{M \to Q}(\mathbf{p}_M)) = - \frac{a_1 a_2 a_3 \nu^2 \omega^3}{\pi^2 \rho^2 V_{0M}^2(\theta_p)} \int_0^{2\pi} \int_0^{\pi} \frac{I_{M \to Q}(\theta_s)}{V_{0Q}^5(\theta_s)} \sin \theta_s d\theta_s d\varphi_s \mathrm{PV} \int_0^{\infty} \frac{\xi^4}{(1-\xi^2)(1+C(\xi))^2} d\xi. \quad (50)
$$

278 As done by Huang et al. [14], an analytical solution can be obtained for the Cauchy Principal Value
279 integral with contour integration and the residue theorem. However, as noted by the authors, it
280 appears that the solution provided in Huang et al. [14], namely Equation (13), yields incorrect
281 results for the phase velocity. Our final formulas for the real part of the mass operator can be
282 written as:

$$
\mathrm{Re}(\bar{\sigma}^{M \to Q}(\mathbf{p}_M)) = - \frac{a_1 a_2 a_3 \nu^2 \omega^3}{\pi^2 \rho^2 V_{0M}^2(\theta_p)} \int_0^{2\pi} \int_0^{\pi} \frac{I_{M \to Q}(\theta_s)}{V_{0Q}^5(\theta_s)} \chi(\theta_s, \varphi_s) \sin \theta_s \, d\theta_s \, d\varphi_s, \quad (51)
$$

283 with,

$$
\chi = - \frac{\pi \left( 2\tilde{C}_0^5 - 4\tilde{C}_0^3 C_1^2 + \tilde{C}_0^2 C_2^2 (14\tilde{C}_0 + 6C_2) + \tilde{C}_0 C_2 (10\tilde{C}_0^3 - 6\tilde{C}_0 C_1^2 - 6C_1^2 C_2) + C_1^4 (\tilde{C}_0 + C_2) \right)}{D^3 (\tilde{C}_0 - C_1^2 + C_2^2)^2 (\tilde{C}_0 + C_1 + C_2)}, \quad (52)
$$

284 where $\tilde{C}_0 = 1+C_0$, and the dependence on the angles of the incident and scattered wave is implicit.

285 The solution in 2D is not as straightforward in 3D because of the 3/2 power, which requires a
286 branch cut in the complex plane; however, numerical integration is possible by introducing a small
287 value $\epsilon \to 0$:

$$
\mathrm{Re}(\bar{\sigma}^{M \to Q}(\mathbf{p}_M)) = - \frac{a_1 a_2 \nu^2 \omega^2}{2\pi \rho^2 V_{0M}^2(\theta_p)} \int_0^{2\pi} \frac{I_{M \to Q}^{2D}(\theta_s)}{V_{0Q}^4(\theta_s)} d\theta_s \int_0^{\infty} \frac{\xi^3}{(1 - \xi^2 + \epsilon) (1 + C(\xi))^{3/2}} d\xi. \quad (53)
$$

288 The details of the process developed to obtain the above analytical solution can be found in
289 Appendix A2: Development of analytical formulas for phase velocity.
## 3. Numerical modeling of polycrystalline media with macroscopic anisotropy

As it has already been explained in the theoretical section, for a fiber-textured medium, the homogeneous reference medium is considered to be transversely isotropic. In this section, we discuss some aspects regarding the FE numerical modeling of this type of anisotropic media. As in our previous work, a discontinuous Galerkin method is used, coupled with an explicit Runge-Kutta time discretization scheme. Readers are referred to [29, 25, 26] for more details.

In this work, four 2D models are considered, which are defined in the $e_x - e_y$ plane. Within the framework of the previously defined theoretical setting, the plane spanned by ($e_x, e_y$) is actually the same as the one spanned by ($e_2, e_3$), and the direction of incident wave propagation is $\hat{p} = e_x$. The four models correspond to four different preferred rotations of the fiber axis ($\hat{n}$) with respect to the wave propagation direction $e_x$. The four considered rotations are (a) 90°, (b) 70°, (c) 60° and (d) 45°. Each rotation is defined with respect to the crystallographic axis (001), which is initially aligned with the $e_3$ direction. Around this axis the (100) direction is arbitrarily rotated. To realize the rotation in the plane of interest ($e_2, e_3$) the Euler angles are specified as follows:: the first Euler angle is set to $\phi_1 = 0$, the second Euler angle $\phi$ is set to the required $\theta$ rotation about $e_3$; and finally, a random rotation is applied for the last Euler angle, $\phi_2$.

Two sets of geometrical parameters, denoted as (I) and (II), are considered. Set (I), used for the case (a), consists of a rectangular domain with a length of $l_x = 25.08$ mm and a width of $l_y = 16.34$ mm. Set (II), used for cases (b), (c), and (d), consists of a rectangular domain with a length of $l_x = 33.06$ mm and a width of $l_y = 16.34$ mm. A summary of the two numerical settings is presented in Table 1.

The top and bottom edges of the numerical setting (I) are subjected to symmetric boundary conditions (SBC), as in our previous studies [29], where the untextured case was considered. In contrast, for the setting (II), the periodic boundary conditions (PBC) are essential and must be applied on these edges. Indeed, in the case of $\theta = 90°$, the use of SBC can still yield accurate estimates of both attenuation and phase velocity, due to the alignment between the textured anisotropy axes and the sample axes, which results in a pure plane wave. Otherwise, a free boundary condition is applied to the left and right edges in both cases.

Equiaxed grains are considered with an average grain diameter of approximately $d \approx 447$ µm. A combined Ricker signal with three central frequencies $f_c = 2, 4.5, 10.5$ MHz is employed as the external pressure loading at the left edge. An element size $h = 0.019$ mm is chosen, satisfying at least 20 elements in the average grain's diameter and a CFL = 0.2 condition, additionally satisfying at least 10 elements in the shortest longitudinal wavelength employed. Additionally, recording nodes were placed along the right-most vertical edge, defined by $x = l_x$.

<table><caption>Table 1. Numerical settings (I) and (II) in 2D used in this work, with $l_i$ (mm) the sample length in each dimension, $h$ (mm) the element size, $f_c$ (MHz) the central frequencies of the applied Ricker signal and the covered range of normalized frequency $\bar{x}_0 = \bar{p}_{0qL}d$, and $d$ the average grain diameter.</caption><thead><tr><th>Numerical Setting</th><th>l<sub>x</sub> &times; l<sub>y</sub></th><th>h</th><th>Num. Grains</th><th>f<sub>c</sub></th><th>x&#772;<sub>0</sub></th></tr></thead><tbody><tr><td>(I)</td><td>25.08 &times; 16.34</td><td>0.019</td><td>2608</td><td>(2, 4.5, 10.5)</td><td>0.4 - 7.7</td></tr><tr><td>(II)</td><td>33.06 &times; 16.34</td><td>0.019</td><td>3439</td><td>(2, 4.5, 10.5)</td><td>0.4 - 7.7</td></tr></tbody></table>

For each case, a 2D numerical simulation is also performed on a homogeneous reference medium, defined by a Voigt average stiffness tensor, occupying the same rectangular domain, and subjected to the same boundary conditions. In the present work, Inconel718 polycrystals are considered. The single cubic elastic constants of the Inconel718 are presented in Table 2. The Voigt average tensor exhibits transversely isotropic symmetry and, in its local basis, is determined by five independent constants. For the case (a), where the axis of symmetry – the fiber axis – is aligned with $e_x$, the five constants of the Voigt average for the Inconel718 are: $C_{11}^0 = C_{33}^0 = 268.6$, $C_{22}^0 = 240.9$, $C_{12}^0 = C_{23}^0 = 140.5$, $C_{13}^0 = 112.8$, $C_{44}^0 = C_{66}^0 = 105.7$ and $C_{55}^0 = (C_{11}^0 - C_{13}^0)/2 = 77.9$ GPa. For each tilt of the fiber axis, this Voigt tensor is rotated accordingly.
<table><caption>Table 2. Properties of the considered material (Inconel718): single cubic elastic constants $C_{ij}$ (GPa), density $\rho$ (kg/m$^3$), and Zener's anisotropy factor $A$.</caption><thead><tr><th></th><th>C<sub>11</sub></th><th>C<sub>12</sub></th><th>C<sub>44</sub></th><th>ρ</th><th>A</th></tr></thead><tbody><tr><td>Inconel718</td><td>240.9</td><td>118.6</td><td>83.8</td><td>8193</td><td>2.1</td></tr></tbody></table>

In contrast to elastic wave propagation in an isotropic medium, in an anisotropic one, the polarization direction is generally not parallel to the wave vector, except when the wave vector is aligned with one of the local material axes. More importantly, there is no direct correspondence between the wave front, the wave vector, and the wave phase velocity [7]. In the specific 2D case of the Inconel718 medium defined above, the wavefront is dependent on the fiber axis tilt. Figure 2 shows the theoretical and numerical wavefront calculated for the Voigt homogeneous stiffness tensor for a 70° rotated fiber axis [20]. The particular wavefront prevents the formation of a perfect plane wave when a uniform line excitation is applied in our classical numerical setting, specifically causing wave effects near the top and bottom boundaries. The difference between a propagating wave in a rectangular polycrystalline medium with SBC and PBC can be observed by comparing Figures 3a and 3b.

<figure><img src="image_3.png" /><figcaption>Figure 2. (a) Theoretical and (b) numerical wave fronts for the $qL$ and $qT_2$ modes in a 2D TI medium with its transverse axis rotated 70°.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3. Propagating plane wave in a polycrystal with macroscopic transverse isotropy with the fiber axis oriented $\theta = 70°$ from the propagation direction $e_x$ with (a) applied SBC and (b) PBC.</figcaption></figure>
This is the main reason for choosing PBC that match the top and bottom edges of the sample in the numerical setting (II). The periodicity condition imposes the velocity and stress continuity of both top and bottom degrees of freedom. This matching condition allows for a pure plane wave to propagate even in the anisotropic conditions treated here. This condition does not accurately represent a practical non-destructive test but is a first step that enables a direct comparison with the theoretical estimations and, more importantly, has a clear and sole analysis of the elastic wave scattering. It should be noted that the imposed periodicity means that grains and their crystallographic orientations are duplicated periodically along the top and bottom edges of the model. In the present study the additional error and uncertainty introduced by this assumption were not evaluated. Consequently, the effect of PBCs on the numerical estimation of the TPC spatial function—and, ultimately, on the derived amplitude attenuation and phase velocity dispersion estimates—remains to be investigated. For more realistic simulations, it would be valuable to move beyond the plane-wave excitation considered here and to examine the directivity patterns generated by an ultrasound-like source in a larger computational domain. This extension lies outside the scope of the present work and is left for future research.

For the models (a), (c), and (d), $N = 9$ simulations were performed, and for the model (b), $N = 6$ simulations were performed. The total number of simulations was enough for a coherent error estimation in comparison with the theoretical predictions, for which the results are presented in more detail in Section 5. In this work, the TPC spatial function is obtained from the numerical Neper microstructures with a multi-exponential fitting method previously proposed by Van Pamel et al. [28]. This procedure and the one followed to calculate the amplitude attenuation and phase velocity from the numerical results is described in our previous work [29].

# 4. Influence of the fiber texture on the elastic wave scattering

In this section, the theoretical results of the elastic wave scattering in macroscopically anisotropic media are analyzed. The influence of the fiber type texture, that is, the alignment of the single crystal [001] axis, on the attenuation and phase velocity will be studied. At first, a polycrystalline microstructure with equiaxed grains is considered. The intention is to analyze the anisotropy induced by the texture presence, which is modeled with the introduction of the Green dyadic for a transversely isotropic medium.

## 4.1. Attenuation: propagation direction dependency

The following results are calculated for an equiaxed grain with diameter $d \approx 447$ µm and for the Inconel 718 medium defined above. The variation of the attenuation with the propagation direction is a fascinating phenomenon when considering macroscopic anisotropy. In practical NDT cases, valuable information about the microstructure can be obtained when performing analysis from multiple directions. To better analyze this phenomenon, let us first make two preliminary observations about the reference homogenized medium, which is transversely isotropic.

First, Figure 4a shows the slowness surface for the reference medium. It clearly illustrates the distinct behavior of the three propagating wave modes $qL$, $qT_2$, and $T_1$, with the wave vector and the axis of symmetry lying in the $e_x - e_y$ plane. It can be observed that, within the irreducible periodic range of $0-90^\circ$, there exist two points – $0^\circ$ and $68.71^\circ$ – where the velocities of the two transverse waves coincide. These two points will have a visible influence on the attenuation, as will be explained subsequently.

Second, Figure 4b presents the shift angle $\psi$ between the wave propagation direction and the polarization direction of the $qL$ mode. The plot exhibits a sinusoidal pattern. The shift helps to visualize the directions along which the wave propagation direction aligns with its polarization direction. In the studied case, two points – at $0^\circ$ and $49.3^\circ$ – can be identified where the curve crosses zero. In both directions, one would expect to observe purely propagating $qL$ and $qT_2$ modes.

Figure 5 presents, for six normalized frequencies $\bar{x}_{0qL}$, the attenuation as a function of the propagation direction $\theta$ for each propagating mode: $qL$, $qT_2$, and $T_1$, shown in three separate plots,
<figure><img src="image_5.png" /><figcaption>Figure 4. (a) Slowness surface of the qL, qT2, and T1 velocities for the material Inconel718, and (b) Angle ψ(θ) of the qL mode representing the shift between the wave propagation direction and its polarization direction against the propagation direction θ.</figcaption></figure>

respectively. The attenuation is normalized with the average wavenumber $\bar{p}_{0qL}$. The average wavenumber is calculated with the average of the Voigt effective wave velocity as first proposed by Turner [27], which is calculated as follows:

$$ \bar{V}_{0M} = \frac{1}{2} \int_{0}^{\pi} V_{0M}(\theta) \sin \theta \, d\theta. \qquad (54) $$

It is noteworthy that, since $\alpha_{qL}/\bar{p}_{0qL} = \alpha_{qL}d/\bar{x}_{0qL}$, the normalization used here is actually a normalization of $\alpha_{qL}d$, the per-grain attenuation, with respect to the normalized frequency $\bar{x}_{0qL}$. This normalization is chosen for visualization purposes, which allows the curves for all considered frequencies to be compared on the same figure.

The plots corresponding to the qL (Figure 5a) and qT2 (Figure 5b) modes also includes the 2D attenuation. As the 2D wave propagation is chosen to be restricted in the plane $e_2 - e_3$, the 2D wave polarization is reduced to the qL and qT2 modes.

Figure 5a shows the attenuation of the qL mode increasing from zero at 0° to a maximum at 90°. Zero attenuation is obtained when the wave propagates along the fiber axis, because in each grain, the wave corresponds to the same propagating pure L mode with $\psi = 0°$ (see Figure 4b). Indeed, along the fiber axis, the medium is comparable to a monocrystal: this axis corresponds to a direction with no stiffness variation, and therefore wave propagation in this direction exhibits no dispersion whatsoever. In contrast, at 90°, the attenuation reaches its maximum value, which is proportional to the normalized frequency, even though we have again a pure propagating modes with $\psi = 0°$ (see Figure 4b). This is a common and expected result, known from the macroscopically isotropic case [29, 31, 15]. In this case, the qL wave is in the plane $e_1 - e_2$, which is the plane where the single crystal axes [001] and [010] are randomly oriented. This behavior is comparable, yet different, to a wave propagating in an untextured medium. The 2D attenuation exhibits similar behavior, with the clear difference of having overall lower values than its 3D counterpart.

Figure 5b depicts the attenuation of the qT2 mode varying from zero to a first maximum around 47.27°, then to a second maximum around 68.71°, before decreasing back to zero at 90°. The non-scattering behavior at 0° is explained by the same arguments as for the qL mode. At 90°, on the other hand, the particle polarization again encounters no stiffness variation in its direction, which is parallel to the fiber axis, thus behaving as a non-dispersive mode. The first peak, which is more pronounced at higher frequencies, appears to converge to the direction of 47.27°, a few degrees away from the first angle of 49.3° identified in our preliminary observations. Turner [27] already noticed that this peak varies with frequency within a range of a few degrees. Indeed, he calculated a variation around 47.1°, which he explains by the inclusion of the polarization vectors and the anisotropic Green function. In contrast with the 3D first peak, the 2D peak does not have
<figure><img src="image_6.png" /><figcaption>Figure 5. Variation of the attenuation of the three different propagating wave modes $qL$, $qT_2$, and $T_1$ with the propagation direction $\theta$ for six different average normalized frequencies $\bar{x}_{0qL}$ (in increasing order). The attenuation is normalized with the average wave number $\bar{p}_{0qL}$.</figcaption></figure>

428

429

430

431

432

433

434

435

a visible variation. The peak in the 2D curves seems to converge at around the same direction as in 3D. The second peak is also frequency-dependent and appears to approach the direction of 68.71°, which corresponds to the second angle identified in our preliminary observations as the intersection point of the slowness surface of the $T_1$ and the $qT_2$ modes (see Figure 4a). Therefore, this peak can be traced to the correlation between the $qT_2$ and $T_1$ modes. The 2D attenuation curves help to demonstrate this hypothesis, as these curves do not exhibit the second peak. Indeed, Figure 5b shows that the peak frequency variation evolves as there is a stronger correlation between the transverse modes as the frequency increases.

436

437

438

439

440

441

Figure 5c presents the attenuation of the $T_1$ mode only in 3D. Similar to the $qL$ mode, the $T_1$ mode (Figure 5c) has non-scattering behavior at 0° and reaches maximum attenuation at 90°. At 90°, $T_1$ mode's particle polarization lies within the plane $e_1 - e_2$, thereby encountering stiffness variations caused by the randomly distributed crystallographic axes during propagation. In this plot, the correlation with the $qT_2$ mode can also be observed at around 68.71°, revealed by a subtle slope change on the curve at the high normalized frequencies, more noticeable above $\bar{x}_{0qL} = 20.0$.

442

## 4.2. Attenuation: frequency dependency

443

444

445

446

Here, the attenuation's frequency dependency is analyzed. The attenuation's behavior with varying frequency is already known from previous works concerning the macroscopically isotropic case [29, 31, 14]. In the following, however, we intend to study the variation between the attenuation curve in different propagating directions.
Figure 6 presents three plots, which correspond respectively to the three propagating modes $qL$, $qT_2$, and $T_1$. Each plot displays the attenuation for three propagation directions $\theta = 49^\circ$, $\theta = 68^\circ$ and $\theta = 90^\circ$. The 2D attenuation is also included in the $qL$ and $qT_2$ modes for comparison. The $qT_2$ plot does not include the curve corresponding to $\theta = 90^\circ$ because, in this direction, the propagating wave does not scatter for this mode, as already explained above.

The $qL$ (Figures 6a) and $T_1$ (6c) modes, as observed above, have an attenuation that increases as the propagation direction $\theta$ approaches $90^\circ$. It can also be observed that the gap between curves is higher at higher frequencies, as expected from the previously discussed results. Additionally, the $qL$ mode curves have the characteristic “hump” observed in the master $L$ attenuation curve in the isotropic case. On the other hand, this hump is not very clearly noticeable for the other two $T_1$ and $qT_2$ modes. Looking at the $qT_2$ attenuation, the absence of the second peak (see Figure 5b) in the 2D attenuation curves can be witnessed as both curves for $\theta = 49^\circ$ and $\theta = 68^\circ$ diverge at high frequencies. It is more interesting to note that the 2D curve of $\theta = 49^\circ$ converges with the 3D curves of $\theta = 49^\circ$ and $\theta = 68^\circ$.

<figure><img src="image_7.png" /><figcaption>Figure 6. Variation of the attenuation $\alpha_M/\bar{p}_{0M}$ of the three different propagating wave modes $qL$, $qT_2$, and $T_1$ with normalized frequency $\bar{p}_{0M}d$. For the modes $qL$ and $qT_2$, the 2D attenuation is also plotted in dashed colored lines.</figcaption></figure>

460

### 4.3. Attenuation: scattering mechanisms

We then proceed to study the particular scattering mechanisms of the different wave modes, which may offer a deeper understanding of attenuation behavior.

The scattering mechanisms of the $qL$ mode (see Figure 7) are comparable to those observed in the isotropic case in both 3D and 2D. In the Rayleigh scattering region, energy is predominantly scattered from the longitudinal to the transverse modes. However, this behavior begins to change in
the Rayleigh-to-stochastic (R-to-S) transition region and is fully achieved in the stochastic region, where scattering predominantly occurs within the same $qL$ modes.

<figure><img src="image_8.png" /><figcaption>Figure 7. Contribution to the total $qL$ normalized attenuation of the different scattering modes in (a) 3D and (b) 2D for $\theta = 90^{\circ}$.</figcaption></figure>

The $qT_2$ mode also behaves similarly to the isotropic case. Figure 8 shows that the scattering along the entire frequency range occurs mainly into the transverse modes. Notably, in 3D, the contribution of the scattering into the $T_1$ mode is larger than that into the $qT_2$ mode. However, this contribution becomes negligible in the stochastic scattering region (see Figure 8a). This conversion gives the $qT_2$ attenuation curve a subtle characteristic hump, as already observed in the longitudinal attenuation curve. In contrast, as expected, the 2D curve (Figure 8b) does not display this hump due to the absence of the out-of-plane scattering. It is noteworthy that the 3D scattering mechanisms of the $qT_2$ mode differ at a propagation angle around $\theta = 68^{\circ}$ (see Figure 8c). As previously observed, around this propagation angle, there is a strong correlation between the $qT_2$ and $T_1$ modes. Figure 8c illustrates that, at $\theta = 68^{\circ}$, both $qT_2 \rightarrow T_1$ and $qT_2 \rightarrow qT_2$ mode conversions prevail in both Rayleigh and stochastic scattering regions. This additional contribution from the $qT_2 \rightarrow T_1$ mode conversion in the stochastic region explains the attenuation peak previously observed in Figure 5b at around $\theta = 68^{\circ}$.

As far as the $T_1$ mode is concerned, Figure 9 shows a mechanism similar to the $qT_2$ mode. In the Rayleigh region, the dominant mode conversion is the one into the transverse modes. In the stochastic region, the dominant mode conversion is the one into the same transverse mode, i.e., only the $T_1$ mode.

### 4.4. Phase Velocity: propagation direction dependency

In this subsection, the phase velocity dependency on the wave propagation direction will be studied. Figure 10 presents the phase velocity variation, $V_M/\bar{V}_{0M} - 1$, as a function of the wave propagation direction $\theta$, shown in three separate plots corresponding to the three wave modes: $qL$, $qT_2$, and $T_1$, respectively.

First, for the coupled $qL$ and $qT_2$ modes, Figures 10a and 10b show only small dispersion across frequencies; however, there is a clear directional dependence. For the $qL$ mode, dispersion becomes more noticeable as $\theta$ approaches $90^{\circ}$, whereas for the $qT_2$ mode, it is more pronounced around $\theta = 45^{\circ}$. Notably, the $qT_2$ phase velocity reaches a minimum near this angle, while the $qL$ phase velocity exhibits a peak around $\theta = 49.3^{\circ}$. For these two modes, the Voigt effective phase velocity represents the scattered phase velocity very well, indicating that the dispersion of both modes due to scattering is very low.

For the $T_1$ mode (see Figure 10c), it exhibits the highest dispersion among the three modes, and this dispersion increases as $\theta$ approaches $90^{\circ}$. In this case, it can be seen that the Voigt effective phase velocity represents the dispersive velocity well only at low $\theta$ angles. Interestingly, the $T_1$
<figure><img src="image_9.png" /><figcaption>Figure 8. Contribution to the total $qT_2$ normalized attenuation of the different scattering modes for $\theta = 49^\circ$ in (a) 3D and (b) 2D, and (c) for $\theta = 68^\circ$ in 3D.</figcaption></figure>

501

502

503

504

505

506

mode exhibits a turning point with respect to the propagation direction, occurring at approximately 68.71°. Again, this point highlights the coupling between the $T_1$ and $qT_2$ modes already observed in the attenuation behavior. Indeed, it also corresponds to a turning point in the phase velocity curves of the $qT_2$ mode (see Figure 10b). At this point, the scattered phase velocity transitions from increasing to decreasing with frequency for the $T_1$ mode, and conversely, from decreasing to increasing for the $qT_2$ mode.

507

### 4.5. Phase Velocity: frequency dependency

508

509

510

511

The frequency dependence of the phase velocity for the three wave modes is analyzed here. Figure 11 presents the normalized phase velocity $V_M/\bar{V}_{0M} - 1$ as a function of the normalized frequency $\bar{x}_{0qL} = \bar{p}_{0qL}d$ for each propagating mode $qL$, $qT_2$ and $T_1$. Each plot contains five groups of curves, with each group corresponding to a different propagation direction $\theta$.

512

513

514

515

516

517

518

519

520

521

We begin by considering the coupled $qL$ and $qT_2$ modes in 3D. The $qL$ mode phase velocity (see Figure 11a) exhibits a clear transition zone between two apparently non-dispersive regions, located around $\bar{k}_{qL} \approx 1$. The starting point of this zone shifts to lower frequencies as $\theta$ increases. The non-dispersive behavior in the Rayleigh region is well known from previous studies on the isotropic case [29, 14], whereas in the stochastic region, it may result from the use of the Born approximation. For all the propagation directions considered, the dispersive phase velocity appears to approach a value very close to the Voigt effective phase velocity in the stochastic region. Moreover, as discussed above, the phase velocity becomes more dispersive as $\theta$ approaches 90°, which corresponds to propagation along the cross-fiber direction. The $qT_2$ mode (see Figure 11b) is less dispersive than the $qL$ one. In contrast to the $qL$ mode, for the $qT_2$ mode, the Voigt effective phase velocity
<figure><img src="image_10.png" /><figcaption>Figure 9. Contribution of the different scattering modes to the total normalized attenuation of the $T_1$ mode for $\theta = 90^\circ$.</figcaption></figure>

provides an excellent approximation of the dispersive velocity in the Rayleigh region rather than in the stochastic region. Around $\bar{p}_{0qT_2} \approx 10$, the phase velocity becomes dispersive, in contrast to the $qL$ mode. However, the onset of this dispersion shifts to lower frequencies as $\theta$ increases, similar to the behavior observed for the $qL$ mode. Though this behavior is not clearly visible in Figure 11b, the effect of the turning point at approximately $\theta = 68.71^\circ$, observed for the $T_1$ mode, is still present. At this point, the phase velocity changes from increasing to decreasing with frequency. Finally, in 2D, both the $qL$ and $qT_2$ phase velocities (see Figures 11a and 11b) are almost non-dispersive – at least less dispersive than in 3D – and are very well approximated by the Voigt effective phase velocity across the entire frequency range.

For the $T_1$ mode (see Figure 11c), as noted above, it is the most dispersed of the three modes. Similar to the $qL$ mode, there is a transition around $\bar{p}_{0T_1} \approx 1$, where the phase velocity shifts from non-dispersive to strongly dispersive behavior. However, this dispersive region does not exhibit consistent behavior across different propagation directions except the one related to the turning point at approximately $68.71^\circ$, as reported above. For directions close to $49^\circ$, a similar to that of the $qL$ mode is observed; namely, the establishment of a non-dispersive zone where the dispersive phase velocity is very well approximated by the Voigt effective phase velocity in the stochastic region. This behavior, however, is not observed for directions close to $90^\circ$.

## 5. Numerical results: theoretical validation and analysis

In this section, we present numerically estimated attenuation and phase velocity for the $qL$ mode, using the 2D FE models previously defined in Section 3.

### 5.1. *Theoretical validation*

Figure 12 compares the numerical estimates with theoretical predictions. Regarding the $qL$ attenuation (Figure 12a), overall, excellent agreement is obtained. For $\theta = 90, 70$, and $60^\circ$, the numerical results slightly overestimate the theoretical ones in the Rayleigh region as the frequency decreases, with the highest relative errors $|\alpha_{qL}^{TH} - \alpha_{qL}^{NUM}|/\alpha_{qL}^{NUM}$ being 21.4%, 24.7% and 52.8%, respectively. The case $\theta = 45^\circ$ is notable for its very low attenuation compared with the other three cases. As a result, the numerical curve in the Rayleigh region tends toward zero, diverging from the theoretical prediction. This low attenuation indicates a minimal fluctuation in the amplitude; therefore, accurately capturing such small variations ideally requires averaging over a larger number of simulations or using a larger sample. In our study, averaging nine simulations was sufficient to achieve a good comparison with the theoretical curve. In general, achieving better convergence in the Rayleigh region for low-attenuation cases like $\theta = 45^\circ$ entails a higher computational cost. In contrast, for the four cases, the numerical curves show good agreement with theory in the R-to-S transition and stochastic regions.
<figure><img src="image_11.png" /><figcaption>Figure 10. Phase velocity variation with respect to the average velocity $\bar{V}_{0M}$ of the three different propagating wave modes: $M = (qL, qT_2, T_1)$, against the propagation direction $\theta$, for six different normalized frequencies $\bar{x}_{0qL} = \bar{p}_{0qL}d$. For the modes $qL$ and $qT_2$, the 2D velocity variation is also plotted in dashed colored lines. The phase velocity of the Voigt effective stiffness tensor is also plotted in a black dotted line for reference.</figcaption></figure>

## 5.2. Statistical observations: the directional Young modulus

$$E_{\hat{n}} = \frac{1}{(\hat{n} \otimes \hat{n}) : \mathbf{S} : (\hat{n} \otimes \hat{n})}, \qquad (55)$$

<footnote n="556">556 Regarding the $qL$ wave phase velocity variation (Figure 12b), excellent agreement is obtained between the theoretical and the numerical results, with the highest relative errors found for $\theta = 90, 70, 60$, and $45^{\circ}$, being 1.0%, 0.65%, 0.45%, and 0.22%, respectively.</footnote>

<footnote n="557">557 It is interesting to note that at $\theta = 45^{\circ}$ and $\theta = 60^{\circ}$ the phase velocity is faster than the Voigt averaged phase velocity. There is then a change to a slower phase velocity for $\theta = 70$ and $\theta = 90^{\circ}$.</footnote>

<footnote n="558">558 The numerical results confirm that the phase velocity dispersion becomes stronger as the wave propagation direction approaches the cross-fiber direction $\theta = 90^{\circ}$. This observation is consistent with the already discussed higher attenuation in the same cross-fiber direction. Both results demonstrate that the highest dispersion occurs along the cross-fiber direction, in which the crystallographic orientations are randomly distributed. Such randomness results in greater variability in grain-to-grain stiffness, therefore leading to stronger wave scattering.</footnote>

<footnote n="559">559 In the anisotropic case, the directional Young's modulus – the Young's modulus in a specific direction defined by the unit normal vector $\hat{n}$ – is calculated as follows:</footnote>

<footnote n="560">560 where $\mathbf{S} = \mathbf{C}^{-1}$ is the compliance tensor. The directional Young's modulus provides an insightful way to correlate the wave phase velocity estimated in each model with the specific preferred fiber</footnote>
<figure><img src="image_12.png" /><figcaption>Figure 11. Phase velocity variation with respect to the average Voigt effective velocity $\bar{V}_{0M}$ of the three different propagating wave modes: $M = (qL, qT_2, T_1)$, against the normalized frequency $\bar{p}_{0M}d$ for five different propagation directions $\theta$. For the modes $qL$ and $qT_2$, the 2D velocity variation is also plotted in dashed colored lines. The Voigt effective phase velocity for each propagation direction is also plotted in a dotted colored line for reference.</figcaption></figure>

axis orientation. To study its distribution in the direction of wave propagation, the directional Young's modulus $E_{\hat{n}}$ was calculated at each grain within the 2D rectangular polycrystalline domain of each numerical model, which contains 2608 grains for the numerical setting (I) and 3439 grains for the numerical setting (II).

Figure 13 displays a violin plot of the directional Young's Modulus $E_{\hat{n}}$ in the wave propagation direction $\hat{n} = e_1$ for each model with different fiber axis orientation (a) 90°, (b) 70°, (c) 60°, and (d) 45°. The violin plot is very effective in visualizing several characteristics, such as the probability density and the mean. For instance, the distribution of each case helps to visualize the grain-to-grain stiffness variations encountered by the propagating wave. The density shows two spread-out peaks toward the extremes and a local minimum between these peaks. This density can be represented by a logit-normal distribution. Therefore, the mean for each $\theta$-dependent distribution is calculated assuming that $E_{\hat{n}}$ follows the aforementioned logit-normal distribution. The $\theta = 45^\circ$ plot shows a narrower variation in the modulus. This indicates a less grain-to-grain variation or lower impedance at grain boundaries, which is in agreement with the previously noted low attenuation. Additionally, examining the mean values suggests that, on average, the fastest wave propagation occurs at $\theta = 45^\circ$, while the slowest occurs at $\theta = 90^\circ$. This result is consistent with the obtained wave phase velocity dispersion curves from $\theta = 45^\circ$ to $\theta = 90^\circ$.
<figure><img src="image_13.png" /><figcaption>Figure 12. Comparison between the 2D theoretical and numerical $qL$ (a) attenuation and (b) phase velocity variation for four different fiber axis orientation.</figcaption></figure>

## 6. Conclusions

In this work, the elastodynamic Dyson equation for the average Green's function in a heterogeneous and anisotropic medium was presented. The solution was generalized to fiber-textured polycrystalline materials with equiaxed grains, following the development of Turner [27] for transversely isotropic media. In particular, analytical expressions were derived for the attenuation coefficient and phase velocity in 2D. Numerical simulations were carried out using a space-discontinuous Galerkin FE method, combined with periodic boundary conditions. This BC setup enabled the simulation of pure plane waves in both tilted TI homogeneous media and polycrystalline media exhibiting macroscopic TI anisotropy. The theoretical and numerical approaches were applied to study wave scattering in a fiber-textured Inconel 718 polycrystalline medium with equiaxed grains.

The theoretical results show the directional dependency of wave scattering on the wave propagation direction with respect to the fiber axis direction. The $qL$ wave mode displays a stronger scattering when the wave propagates in the cross-fiber direction, that is, in the direction perpendicular to the fiber axis. This behavior is reflected in higher amplitude attenuation and more phase velocity dispersion. As the propagation direction shifts toward the fiber axis, scattering decreases monotonically, reaching a minimum – non-scattered behavior – A similar pattern is observed for the $T_1$ mode (shear-horizontal mode). However, the $qT_2$ mode (shear-vertical mode) shows a different behavior: it exhibits two distinct scattering peaks in 3D, and a single peak in 2D. These peaks occur between 0 and 90°, indicating that, in contrast with the other modes, the $qT_2$ mode does not scatter when the wave propagates either parallel or perpendicular to the fiber axis.

The numerical results confirmed both the directional and frequency dependency of attenuation and phase velocity in 2D fiber-textured TI media. In terms of phase velocity, the presence of both fast and slow propagation directions was revealed theoretically and numerically. For example, the $qL$ mode was observed to propagate fastest at approximately $\theta \approx 49^\circ$.

Furthermore, through a dimensionality analysis, it was observed that the 2D scattering in macroscopically anisotropic media lacks the contribution of the out-of-plane shear mode, which is essential to accurately predict the 3D behavior. This was demonstrated by analyzing the attenuation mech-
<figure><img src="image_14.png" /><figcaption>Figure 13. Distribution of the directional Young's modulus $E_{\hat{n}}$ in the $\hat{n} = e_1$ direction for a set of 2608 grains for the four different fiber axis rotations.</figcaption></figure>

anisms, showing that in 3D, there is a significant contribution of mode conversion to the $T_1$ mode, which is, in the case of the $qL$ mode, even higher than the mode conversion to the $qT_2$ mode. The present work was a comprehensive study of the amplitude attenuation and phase velocity in fiber-textured media, by comparing the 2D and 3D phenomena, and utilizing numerical approaches to validate the newly derived 2D equations. This work provides insight into the ultrasonic scattering mechanisms, which can be useful for NDT methodologies for WLAM components, by providing knowledge on the directionality dependence of the wave scattering direction.

### CRediT authorship contribution statement

Juan Camilo Victoria Giraldo: Writing – original draft, review & editing, Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Data Curation, Visualization.

Denis Solas: Writing – review & editing, Conceptualization, Validation, Investigation, Supervision, Funding acquisition. Jérôme Laurent: Writing – review & editing, Conceptualization, Validation, Supervision, Funding acquisition. Alain Lhémery: Writing – review & editing, Conceptualization, Validation, Supervision. Bing Tie: Writing – original draft, review & editing, Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Supervision, Funding acquisition.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgments

This work was part of Juan Camilo Victoria Giraldo's doctoral study. This research work was funded/supported by the Agence Nationale de la Recherche, France (No. ANR-21-CE08-0026, Research Project COLUMBO). Computations were performed using computational resources from the "Mésocentre" computing center of Université Paris-Saclay, CentraleSupélec and École Normale Supérieure Paris-Saclay supported by CNRS and Région Île-de-France (https://mesocentre.universite-paris-saclay.fr/).

### References

[1] S. Ahmed and R. B. Thompson. Propagation of elastic waves in equiaxed stainless-steel polycrystals with aligned [001] axes. J. Acoust. Soc. Am., 99(4):2086-2096, April 1996. ISSN 0001-4966, 1520-8524. doi: 10.1121/1.415395. URL https://pubs.aip.org/jasa/article/99/4/2086/761462/Propagation-of-elastic-waves-in-equiaxed-stainless.
[2] Frédéric Allaire, Filippo Gatti, Roberto Miorelli, Vincent Dorval, and Jérôme Laurent. Metamodeling elastic wave propagation using a mixed factorized Fourier encoder-decoder for online laser-ultrasound testing in additive manufacturing. 2025. hal-05011076.

[3] Andrea P. Arguelles and Joseph A. Turner. Ultrasonic attenuation of polycrystalline materials with a distribution of grain sizes. J. Acoust. Soc. Am., 141(6):4347-4353, June 2017. ISSN 0001-4966, 1520-8524. doi: 10.1121/1.4984290. URL https://pubs.aip.org/jasa/article/141/6/4347/915703/Ultrasonic-attenuation-of-polycrystalline.

[4] Andrea P. Arguelles, Christopher M. Kube, Ping Hu, and Joseph A. Turner. Mode-converted ultrasonic scattering in polycrystals with elongated grains. J. Acoust. Soc. Am., 140(3):1570-1580, September 2016. ISSN 0001-4966. doi: 10.1121/1.4962161. URL http://asa.scitation.org/doi/10.1121/1.4962161.

[5] X. Bai, B. Tie, J.-H. Schmitt, and D. Aubry. Comparison of ultrasonic attenuation within two- and three-dimensional polycrystalline media. Ultrasonics, 100:105980, January 2020. ISSN 0041624X. doi: 10.1016/j.ultras.2019.105980. URL https://linkinghub.elsevier.com/retrieve/pii/S0041624X18308096.

[6] Marie Calvet and Ludovic Margerin. Impact of grain shape on seismic attenuation and phase velocity in cubic polycrystalline materials. Wave Motion, 65:29-43, September 2016. ISSN 01652125. doi: 10.1016/j.wavemoti.2016.04.001. URL https://linkinghub.elsevier.com/retrieve/pii/S0165212516300063.

[7] J. M. Carcione, D. Kosloff, and R. Kosloff. Wave-Propagation Simulation in an Elastic Anisotropic (Transversely Isotropic) Solid. Q J Mechanics Appl Math, 41(3):319-346, 1988. ISSN 0033-5614, 1464-3855. doi: 10.1093/qjam/41.3.319. URL https://academic.oup.com/qjam/article-lookup/doi/10.1093/qjam/41.3.319.

[8] D. Dence and J. E. Spence. Wave propagation in random anisotropic media. In Albert T. Bharucha-Reid, editor, Probabilistic methods in applied mathematics, volume 3. Academic Press, New York, 1973. ISBN 978-0-12-095703-3.

[9] Clément Du Burck and Arnaud Derode. Attenuation and dispersion of coherent Rayleigh and head waves on the surface of a polycrystal. J. Acoust. Soc. Am., 159(1):955-973, January 2026. ISSN 1520-8524. doi: 10.1121/10.0042276. URL https://pubs.aip.org/jasa/article/159/1/955/3378305/Attenuation-and-dispersion-of-coherent-Rayleigh.

[10] U. Frisch. Wave propagation in random media. In Albert T. Bharucha-Reid, editor, Probabilistic methods in applied mathematics, volume 1. Academic Press, New York, 1968. ISBN 978-0-12-095703-1.

[11] Cécile Gueudre, Loïc Le Marrec, Joseph Moysan, and Bertrand Chassignole. Direct model optimisation for data inversion. Application to ultrasonic characterisation of heterogeneous welds. NDT & E International, 42(1):47-55, January 2009. ISSN 09638695. doi: 10.1016/j.ndteint.2008.07.003. URL https://linkinghub.elsevier.com/retrieve/pii/S0963869508000753.

[12] Sigrun Hirsekorn. The scattering of ultrasonic waves in polycrystalline materials with texture. J. Acoust. Soc. Am., 77(3):832-843, March 1985. ISSN 0001-4966, 1520-8524. doi: 10.1121/1.392052. URL https://pubs.aip.org/jasa/article/77/3/832/802050/The-scattering-of-ultrasonic-waves-in.

[13] Sigrun Hirsekorn. Directional dependence of ultrasonic propagation in textured polycrystals. J. Acoust. Soc. Am., 79(5):1269-1279, May 1986. ISSN 0001-4966, 1520-8524. doi: 10.1121/1.393706. URL https://pubs.aip.org/jasa/article/79/5/1269/680026/Directional-dependence-of-ultrasonic-propagation.

[14] M. Huang, G. Sha, P. Huthwaite, S. I. Rokhlin, and M. J. S. Lowe. Elastic wave velocity dispersion in polycrystals with elongated grains: Theoretical and numerical analysis. J. Acoust. Soc. Am., 148(6):3645-3662, December 2020. ISSN 0001-4966. doi: 10.1121/10.0002916. URL http://asa.scitation.org/doi/10.1121/10.0002916.

[15] M. Huang, G. Sha, P. Huthwaite, S. I. Rokhlin, and M. J. S. Lowe. Longitudinal wave attenuation in polycrystals with elongated grains: 3D numerical and analytical modeling. J. Acoust. Soc. Am., 149(4):2377-2394, April 2021. ISSN 0001-4966. doi: 10.1121/10.0003955. URL https://asa.scitation.org/doi/10.1121/10.0003955.

[16] Frank C. Karal and Joseph B. Keller. Elastic, Electromagnetic, and Other Waves in a Random Medium. J. Math. Phys., 5(4):537-547, April 1964. ISSN 0022-2488, 1089-7658. doi: 10.1063/1.1704145. URL http://aip.scitation.org/doi/10.1063/1.1704145.

[17] Zichen Kong, Gildas Guillemot, Michel Bellet, and Charles-André Gandin. Multiphysics simulation and microstructure prediction of coaxial wire-laser additive manufacturing process. Materialia, 42:102461, August 2025. ISSN 25891529. doi: 10.1016/j.mtla.2025.102461. URL https://linkinghub.elsevier.com/retrieve/pii/S2589152925001292.

[18] J. Li and S.I. Rokhlin. Propagation and scattering of ultrasonic waves in polycrystals with arbitrary crystallite and macroscopic texture symmetries. Wave Motion, 58:145-164, November 2015. ISSN 01652125. doi: 10.1016/j.wavemoti.2015.05.004. URL https://linkinghub.elsevier.com/retrieve/pii/S0165212515000724.

[19] J Moysan, A Apfel, G Corneloup, and B Chassignole. Modelling the grain orientation of austenitic stainless steel multipass welds to improve ultrasonic assessment of structural integrity. Int. J. Press. Vessel. Pip., 80 (2):77-85, February 2003. ISSN 03080161. doi: 10.1016/S0308-0161(03)00024-3. URL https://linkinghub.elsevier.com/retrieve/pii/S0308016103000243.

[20] Robert G. Payton. Elastic wave propagation in transversely isotropic media, volume 4 of Mechanics of elastic and inelastic solids. Springer Netherlands, Dordrecht, 1983. ISBN 978-94-009-6868-4 978-94-009-6866-0. doi: 10.1007/978-94-009-6866-0. URL http://link.springer.com/10.1007/978-94-009-6866-0.

[21] Anubhav Roy. Elastic Wave Propagation and Scattering in Polycrystals. PhD thesis, The Pennsylvania State University, 2023.

[22] Anubhav Roy and Christopher M. Kube. Multiple scattering of elastic waves in polycrystals. J. Mech. Phys. Solids, 203:106237, October 2025. ISSN 00225096. doi: 10.1016/j.jmps.2025.106237. URL https://linkinghub.elsevier.com/retrieve/pii/S0022509625002133.

[23] Ningyue Sheng, Shahram Khazaie, Mathilde Chevreuil, and Sylvain Fréour. Influence of grain size distribution type on statistical moments of ultrasonic attenuations and phase velocities in random polycrystals. Ultrasonics, 157:107800, September 2025. ISSN 0041624X. doi: 10.1016/j.ultras.2025.107800. URL https://linkinghub.
717 elsevier.com/retrieve/pii/S0041624X25002379.

718 [24] Fred E. Stanke and G. S. Kino. A unified theory for elastic wave propagation in polycrystalline materials. J. Acoust. Soc. Am., 75(3):665-681, March 1984. ISSN 0001-4966. doi: 10.1121/1.390577. URL http://asa.scitation.org/doi/10.1121/1.390577.

719 [25] 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. Comput. Methods Appl. Mech. Eng., 372:113352, December 2020. ISSN 00457825. doi: 10.1016/j.cma.2020.113352. URL https://linkinghub.elsevier.com/retrieve/pii/S0045782520305375.

720 [26] 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:299-332, August 2018. ISSN 00457825. doi: 10.1016/j.cma.2018.04.018. URL https://linkinghub.elsevier.com/retrieve/pii/S004578251631845X.

721 [27] Joseph A. Turner. Elastic wave propagation and scattering in heterogeneous, anisotropic media: Textured polycrystalline materials. J. Acoust. Soc. Am., 106(2):541-552, August 1999. ISSN 0001-4966, 1520-8524. doi: 10.1121/1.427024. URL https://pubs.aip.org/jasa/article/106/2/541/550659/Elastic-wave-propagation-and-scattering-in.

722 [28] A. Van Pamel, G. Sha, M. J. S. Lowe, and S. I. Rokhlin. Numerical and analytic modelling of elastodynamic scattering within polycrystalline materials. J. Acoust. Soc. Am., 143(4):2394-2408, April 2018. ISSN 0001-4966. doi: 10.1121/1.5031008. URL http://asa.scitation.org/doi/10.1121/1.5031008.

723 [29] Juan Camilo Victoria-Giraldo, Bing Tie, Jérôme Laurent, Alain Lhémery, and Denis Solas. Ultrasonic scattering in polycrystalline materials with elongated grains: A comparative 3D and 2D theoretical and numerical analysis. Ultrasonics, 152:107642, August 2025. ISSN 0041624X. doi: 10.1016/j.ultras.2025.107642. URL https://linkinghub.elsevier.com/retrieve/pii/S0041624X25000794.

724 [30] R.L. Weaver. Diffusivity of ultrasound in polycrystals. J. Mech. Phys. Solids, 38(1):55-86, January 1990. ISSN 00225096. doi: 10.1016/0022-5096(90)90021-U. URL https://linkinghub.elsevier.com/retrieve/pii/002250969090021U.

725 [31] L. Yang, O.I. Lobkis, and S.I. Rokhlin. Shape effect of elongated grains on ultrasonic attenuation in polycrystalline materials. Ultrasonics, 51(6):697-708, August 2011. ISSN 0041624X. doi: 10.1016/j.ultras.2011.02.002. URL https://linkinghub.elsevier.com/retrieve/pii/S0041624X11000266.

726 [32] W. Y. Yeoh, B. Lan, and M. J. S. Lowe. Investigation of the influence of macrozones in titanium alloys on the propagation and scattering of ultrasound. Proc. R. Soc. A., 479(2277):20230176, September 2023. ISSN 1364-5021, 1471-2946. doi: 10.1098/rspa.2023.0176. URL https://royalsocietypublishing.org/doi/10.1098/rspa.2023.0176.

727 [33] Wei Yi Yeoh, Bo Lan, and Michael J. S. Lowe. Representative microstructures for two-dimensional computational studies of ultrasonic wave propagation in titanium alloys. Proc. R. Soc. A., 481(2323):20240776, October 2025. ISSN 1364-5021, 1471-2946. doi: 10.1098/rspa.2024.0776. URL https://royalsocietypublishing.org/doi/10.1098/rspa.2024.0776.

## 754 Appendix A1: The Born series

755 Successively iterating the right-hand side of the preceding equation with respect to $G$ gives rise to the following Born series,

$$
\begin{align}
G(x, x') &= G^0(x - x') \nonumber \\
&\quad - \int G^0(x - x_1) \cdot \mathcal{L}^1(x_1) G^0(x_1 - x') dx_1 \nonumber \\
&\quad + \iint G^0(x - x_1) \cdot \mathcal{L}^1(x_1) (G^0(x_1 - x_2) \cdot \mathcal{L}^1(x_2) G^0(x_2 - x')) dx_1 dx_2 \nonumber \\
&\quad - \iiint G^0(x - x_1) \cdot \mathcal{L}^1(x_1) (G^0(x_1 - x_2) \cdot \mathcal{L}^1(x_2) \cdot \mathcal{L}^1(x_3) (G^0(x_2 - x_3) \cdot \mathcal{L}^1(x_3) G^0(x_3 - x')) dx_1 dx_2 dx_3 \tag{56} \\
&\quad + \iiint\int G^0(x - x_1) \cdot \mathcal{L}^1(x_1) (G^0(x_1 - x_2) \cdot \mathcal{L}^1(x_2) (G^0(x_2 - x_3) \cdot \mathcal{L}^1(x_3) (G^0(x_3 - x_4) \cdot \mathcal{L}^1(x_4) G^0(x_4 - x')) dx_1 dx_2 dx_3 dx_4 \nonumber \\
&\quad + \dots \nonumber \\
&\equiv G^0(x - x') + G^{(1)}(x, x') + G^{(2)}(x, x') + G^{(3)}(x, x') + G^{(4)}(x, x') + \dots, \nonumber
\end{align}
$$

757 where $G^{(p)}(x, x')$ is the $p$-th-order perturbation in the Born series.

758 Now, by using the divergence theorem, we perform some development of different terms of the Born series.
• $G^{(2)}(x, x')$:

$$ 
\begin{aligned}
G^{(2)}(x, x') &= \iint G^0(x, x_1) \cdot \mathcal{L}^1(x_1) (G^0(x_1, x_2) \cdot \mathcal{L}^1(x_2) G^0(x_2, x')) \,dx_1dx_2 \\
&= - \iint G^0(x, x_1) \cdot \nabla_{x_1} \cdot ((\delta C(x_1) \otimes \delta C(x_2)) : (\nabla^2_{x_1x_2} G^0(x_1, x_2) \otimes \nabla_{x_2} G^0(x_2, x'))) \,dx_1dx_2
\end{aligned} 
\quad (57) $$

Or with indicial notations:

$$ 
\begin{aligned}
G_{mn}^{(2)}(x, x') &= \iint G_{mi_1}^0(x, x_1) \frac{\partial}{\partial x_{1j_1}} \left( \delta C_{i_1 j_1 k_1 l_1}(x_1) \frac{\partial}{\partial x_{1l_1}} \right. \\
& \qquad \left. \left( G_{k_1 i_2}^0(x_1, x_2) \frac{\partial}{\partial x_{2j_2}} \left( \delta C_{i_2 j_2 k_2 l_2}(x_2) \frac{\partial G_{k_2 n}^0(x_2, x')}{\partial x_{2l_2}} \right) \right) \right) dx_1 dx_2 \\
&= - \iint G_{mi_1}^0(x - x_1) \frac{\partial}{\partial x_{1j_1}} \left( \delta C_{i_1 j_1 k_1 l_1}(x_1) \delta C_{i_2 j_2 k_2 l_2}(x_2) \frac{\partial^2 G_{k_1 i_2}^0(x_1 - x_2)}{\partial x_{1l_1} \partial x_{2j_2}} \frac{\partial G_{k_2 n}^0(x_2, x')}{\partial x_{2l_2}} \right) dx_1 dx_2.
\end{aligned} 
\quad (58) $$

The details of the calculations are:

$$ 
\begin{aligned}
G^{(2)}(x, x') &= \iint G^0(x, x_1) \cdot \mathcal{L}^1(x_1) (G^0(x_1, x_2) \cdot \mathcal{L}^1(x_2) G^0(x_2, x')) \,dx_1dx_2 \\
&= \iint G^0(x, x_1) \cdot \nabla_{x_1} \cdot (\delta C(x_1) : \nabla_{x_1} \\
& \qquad (G^0(x_1, x_2) \cdot (\nabla_{x_2} \cdot (\delta C(x_2) : \nabla_{x_2} G^0(x_2, x'))))) \,dx_1dx_2 \\
&= - \iint G^0(x, x_1) \cdot \nabla_{x_1} \cdot (\delta C(x_1) : \nabla_{x_1} \\
& \qquad (\nabla_{x_2} G^0(x_1, x_2) : \delta C(x_2) : \nabla_{x_2} G^0(x_2, x'))) \,dx_1dx_2 \\
&\quad + \iint G^0(x, x_1) \cdot \nabla_{x_1} \cdot (\delta C(x_1) : \nabla_{x_1} \\
& \qquad (G^0(x_1, x_2) \otimes n_2 : \delta C(x_2) : \nabla_{x_2} G^0(x_2, x'))) \,dx_1dS_2 \\
&= - \iint G^0(x, x_1) \cdot \nabla_{x_1} \cdot ((\delta C(x_1) \otimes \delta C(x_2)) : (\nabla^2_{x_1x_2} G^0(x_1, x_2) \otimes \nabla_{x_2} G^0(x_2, x'))) \,dx_1dx_2 \\
&\quad + \iint G^0(x, x_1) \cdot \nabla_{x_1} \cdot ((\delta C(x_1) \otimes \delta C(x_2)) : ((\nabla_{x_1} G^0(x_1, x_2) \otimes n_2) \otimes \nabla_{x_2} G^0(x_2, x'))) \,dx_1dS_2
\end{aligned} 
\quad (59) $$
Or with indicial notations:

$$
\begin{align*}
G_{mn}^{(2)}(\mathbf{x}, \mathbf{x}') &= \iint G_{mi_1}^0(\mathbf{x}, \mathbf{x}_1) \frac{\partial}{\partial x_{1j_1}} \left( \delta C_{i_1 j_1 k_1 l_1}(\mathbf{x}_1) \frac{\partial}{\partial x_{1l_1}} \right. \\
&\qquad \left. \left( G_{k_1 i_2}^0(\mathbf{x}_1, \mathbf{x}_2) \frac{\partial}{\partial x_{2j_2}} \left( \delta C_{i_2 j_2 k_2 l_2}(\mathbf{x}_2) \frac{\partial G_{k_2 n}^0(\mathbf{x}_2, \mathbf{x}')}{\partial x_{2l_2}} \right) \right) \right) \mathrm{d}\mathbf{x}_1 \mathrm{d}\mathbf{x}_2 \\
&= - \iint G_{mi_1}^0(\mathbf{x}, \mathbf{x}_1) \frac{\partial}{\partial x_{1j_1}} \left( \delta C_{i_1 j_1 k_1 l_1}(\mathbf{x}_1) \frac{\partial}{\partial x_{1l_1}} \right. \\
&\qquad \left. \left( \frac{\partial G_{k_1 i_2}^0(\mathbf{x}_1, \mathbf{x}_2)}{\partial x_{2j_2}} \delta C_{i_2 j_2 k_2 l_2}(\mathbf{x}_2) \frac{\partial G_{k_2 n}^0(\mathbf{x}_2, \mathbf{x}')}{\partial x_{2l_2}} \right) \right) \mathrm{d}\mathbf{x}_1 \mathrm{d}\mathbf{x}_2 \\
&\quad + \iint G_{mi_1}^0(\mathbf{x}, \mathbf{x}_1) \frac{\partial}{\partial x_{1j_1}} \left( \delta C_{i_1 j_1 k_1 l_1}(\mathbf{x}_1) \right. \\
&\qquad \left. \left( n_{2j_2} \frac{\partial G_{k_1 i_2}^0(\mathbf{x}_1, \mathbf{x}_2)}{\partial x_{1l_1}} \delta C_{i_2 j_2 k_2 l_2}(\mathbf{x}_2) \frac{\partial G_{k_2 n}^0(\mathbf{x}_2, \mathbf{x}')}{\partial x_{2l_2}} \right) \right) \mathrm{d}\mathbf{x}_1 \mathrm{d}\mathbf{x}_2 \\
&= - \iint G_{mi_1}^0(\mathbf{x}, \mathbf{x}_1) \frac{\partial}{\partial x_{1j_1}} \left( \delta C_{i_1 j_1 k_1 l_1}(\mathbf{x}_1) \frac{\partial^2 G_{k_1 i_2}^0(\mathbf{x}_1, \mathbf{x}_2)}{\partial x_{1l_1} \partial x_{2j_2}} \frac{\partial G_{k_2 n}^0(\mathbf{x}_2, \mathbf{x}')}{\partial x_{2l_2}} \right) \mathrm{d}\mathbf{x}_1 \mathrm{d}\mathbf{x}_2 \\
&\quad + \iint G_{mi_1}^0(\mathbf{x}, \mathbf{x}_1) \frac{\partial}{\partial x_{1j_1}} \left( \delta C_{i_1 j_1 k_1 l_1}(\mathbf{x}_1) \delta C_{i_2 j_2 k_2 l_2}(\mathbf{x}_2) \frac{\partial G_{k_1 i_2}^0(\mathbf{x}_1, \mathbf{x}_2)}{\partial x_{1l_1}} n_{2j_2} \frac{\partial G_{k_2 n}^0(\mathbf{x}_2, \mathbf{x}')}{\partial x_{2l_2}} \right) \mathrm{d}\mathbf{x}_1 \mathrm{d}\mathbf{s}_2
\end{align*}
$$

where the surface integration resulting from the divergence theorem is zero because the dyadic Green's function can be considered to have compact support.

$$
\bullet G^{(3)}(\mathbf{x}, \mathbf{x}):
$$

$$
\begin{align}
G^{(3)}(\mathbf{x}, \mathbf{x}') &= - \iiint G^0(\mathbf{x}, \mathbf{x}_1) \cdot \mathcal{L}^1(\mathbf{x}_1) (G^0(\mathbf{x}_1 - \mathbf{x}_2) \cdot \nonumber \\
&\quad \mathcal{L}^1(\mathbf{x}_2) (G^0(\mathbf{x}_2 - \mathbf{x}_3) \cdot \mathcal{L}^1(\mathbf{x}_3) G^0(\mathbf{x}_3 - \mathbf{x}')) \mathrm{d}\mathbf{x}_1 \mathrm{d}\mathbf{x}_2 \mathrm{d}\mathbf{x}_3 \nonumber \\
&= - \iiint G^0(\mathbf{x}, \mathbf{x}_1) \cdot \nabla_{\mathbf{x}_1} \cdot (\delta C(\mathbf{x}_1) \otimes \delta C(\mathbf{x}_2) \otimes \delta C(\mathbf{x}_3) : \nonumber \\
&\quad (\nabla^2_{\mathbf{x}_1 \mathbf{x}_2} G^0(\mathbf{x}_1, \mathbf{x}_2) \otimes \nabla^2_{\mathbf{x}_2 \mathbf{x}_3} G^0(\mathbf{x}_2, \mathbf{x}_3) \otimes \nabla_{\mathbf{x}_3} G^0(\mathbf{x}_3, \mathbf{x}')) \mathrm{d}\mathbf{x}_1 \mathrm{d}\mathbf{x}_2 \mathrm{d}\mathbf{x}_3, \tag{59}
\end{align}
$$

Or with indicial notations:

$$
\begin{align}
G_{mn}^{(3)}(\mathbf{x}, \mathbf{x}') &= - \iint G_{mi_1}^0(\mathbf{x}, \mathbf{x}_1) \frac{\partial}{\partial x_{1j_1}} \left( \delta C_{i_1 j_1 k_1 l_1}(\mathbf{x}_1) \frac{\partial}{\partial x_{1l_1}} \left( G_{k_1 i_2}^0(\mathbf{x}_1, \mathbf{x}_2) \frac{\partial}{\partial x_{2j_2}} \left( \delta C_{i_2 j_2 k_2 l_2}(\mathbf{x}_2) \frac{\partial}{\partial x_{2l_2}} \right) \right) \right. \nonumber \\
&\qquad \left. \left( G_{k_2 i_3}^0(\mathbf{x}_2, \mathbf{x}_3) \frac{\partial}{\partial x_{3j_3}} \left( \delta C_{i_3 j_3 k_3 l_3}(\mathbf{x}_3) \frac{\partial G_{k_3 n}^0(\mathbf{x}_3, \mathbf{x}')}{\partial x_{3l_3}} \right) \right) \right) \mathrm{d}\mathbf{x}_1 \mathrm{d}\mathbf{x}_2 \mathrm{d}\mathbf{x}_3 \nonumber \\
&= - \iint G_{mi_1}^0(\mathbf{x}, \mathbf{x}_1) \frac{\partial}{\partial x_{1j_1}} \left( \delta C_{i_1 j_1 k_1 l_1}(\mathbf{x}_1) \delta C_{i_2 j_2 k_2 l_2}(\mathbf{x}_2) \delta C_{i_3 j_3 k_3 l_3}(\mathbf{x}_3) \right. \nonumber \\
&\qquad \left. \frac{\partial^2 G_{k_1 i_2}^0(\mathbf{x}_1, \mathbf{x}_2)}{\partial x_{1l_1} \partial x_{2j_2}} \frac{\partial^2 G_{k_2 i_3}^0(\mathbf{x}_2, \mathbf{x}_3)}{\partial x_{2l_2} \partial x_{3j_3}} \frac{\partial G_{k_3 n}^0(\mathbf{x}_3, \mathbf{x}')}{\partial x_{3l_3}} \right) \mathrm{d}\mathbf{x}_1 \mathrm{d}\mathbf{x}_2 \mathrm{d}\mathbf{x}_3 \tag{60}
\end{align}
$$
The details of the calculations are:

$$
\begin{align*}
G^{(3)}(\mathbf{x}, \mathbf{x}') &= - \iiint G^0(\mathbf{x}, \mathbf{x}_1) \mathcal{L}^1(\mathbf{x}_1) (G^0(\mathbf{x}_1 - \mathbf{x}_2) \\
&\qquad \mathcal{L}^1(\mathbf{x}_2) (G^0(\mathbf{x}_2 - \mathbf{x}_3) \mathcal{L}^1(\mathbf{x}_3) G^0(\mathbf{x}_3 - \mathbf{x}'))) \, d\mathbf{x}_1 d\mathbf{x}_2 d\mathbf{x}_3 \\
&= \iiint G^0(\mathbf{x}, \mathbf{x}_1) \cdot \nabla_{\mathbf{x}_1} \cdot (\delta C(\mathbf{x}_1) : \nabla_{\mathbf{x}_1} (G^0(\mathbf{x}_1, \mathbf{x}_2) \cdot \nabla_{\mathbf{x}_2} \\
&\qquad ((\delta C(\mathbf{x}_2) \otimes \delta C(\mathbf{x}_3)) : (\nabla^2_{\mathbf{x}_2 \mathbf{x}_3} G^0(\mathbf{x}_2, \mathbf{x}_3) \otimes \nabla_{\mathbf{x}_3} G^0(\mathbf{x}_3, \mathbf{x}')))) \, d\mathbf{x}_1 d\mathbf{x}_2 d\mathbf{x}_3 \\
&= - \iiint G^0(\mathbf{x}, \mathbf{x}_1) \cdot \nabla_{\mathbf{x}_1} \cdot (\delta C(\mathbf{x}_1) : \nabla_{\mathbf{x}_1} (\nabla_{\mathbf{x}_2} G^0(\mathbf{x}_1, \mathbf{x}_2) : \\
&\qquad (\delta C(\mathbf{x}_2) \otimes \delta C(\mathbf{x}_3)) : (\nabla^2_{\mathbf{x}_2 \mathbf{x}_3} G^0(\mathbf{x}_2, \mathbf{x}_3) \otimes \nabla_{\mathbf{x}_3} G^0(\mathbf{x}_3, \mathbf{x}')))) \, d\mathbf{x}_1 d\mathbf{x}_2 d\mathbf{x}_3 \\
&= - \iiint G^0(\mathbf{x}, \mathbf{x}_1) \cdot \nabla_{\mathbf{x}_1} \cdot (\delta C(\mathbf{x}_1) \otimes \delta C(\mathbf{x}_2) \otimes \delta C(\mathbf{x}_3) : \\
&\qquad (\nabla^2_{\mathbf{x}_1 \mathbf{x}_2} G^0(\mathbf{x}_1, \mathbf{x}_2) \otimes \nabla^2_{\mathbf{x}_2 \mathbf{x}_3} G^0(\mathbf{x}_2, \mathbf{x}_3) \otimes \nabla_{\mathbf{x}_3} G^0(\mathbf{x}_3, \mathbf{x}'))) \, d\mathbf{x}_1 d\mathbf{x}_2 d\mathbf{x}_3,
\end{align*}
$$

where the surface integration terms are suppressed because they are equal to zero. Or with
indicial notations:

$$
\begin{align*}
G_{mn}^{(3)}(\mathbf{x}, \mathbf{x}') &= - \iint G_{mi_1}^0(\mathbf{x}, \mathbf{x}_1) \frac{\partial}{\partial x_{1j_1}} \left( \delta C_{i_1 j_1 k_1 l_1}(\mathbf{x}_1) \frac{\partial}{\partial x_{1l_1}} \left( G_{k_1 i_2}^0(\mathbf{x}_1, \mathbf{x}_2) \frac{\partial}{\partial x_{2j_2}} \left( \delta C_{i_2 j_2 k_2 l_2}(\mathbf{x}_2) \frac{\partial}{\partial x_{2l_2}} \right. \right. \\
&\qquad \left. \left. \left( G_{k_2 i_3}^0(\mathbf{x}_2, \mathbf{x}_3) \frac{\partial}{\partial x_{3j_3}} \left( \delta C_{i_3 j_3 k_3 l_3}(\mathbf{x}_3) \frac{\partial G_{k_3 n}^0(\mathbf{x}_3, \mathbf{x}')}{\partial x_{3l_3}} \right) \right) \right) \right) \right) d\mathbf{x}_1 d\mathbf{x}_2 d\mathbf{x}_3 \\
&= \iint G_{mi_1}^0(\mathbf{x}, \mathbf{x}_1) \frac{\partial}{\partial x_{1j_1}} \left( \delta C_{i_1 j_1 k_1 l_1}(\mathbf{x}_1) \frac{\partial}{\partial x_{1l_1}} \left( G_{k_1 i_2}^0(\mathbf{x}_1, \mathbf{x}_2) \frac{\partial}{\partial x_{2j_2}} \left( \delta C_{i_2 j_2 k_2 l_2}(\mathbf{x}_2) \frac{\partial}{\partial x_{2l_2}} \right. \right. \\
&\qquad \left. \left. \left( \frac{\partial G_{k_2 i_3}^0(\mathbf{x}_2, \mathbf{x}_3)}{\partial x_{3j_3}} \delta C_{i_3 j_3 k_3 l_3}(\mathbf{x}_3) \frac{\partial G_{k_3 n}^0(\mathbf{x}_3, \mathbf{x}')}{\partial x_{3l_3}} \right) \right) \right) \right) d\mathbf{x}_1 d\mathbf{x}_2 d\mathbf{x}_3 \\
&= \iint G_{mi_1}^0(\mathbf{x}, \mathbf{x}_1) \frac{\partial}{\partial x_{1j_1}} \left( \delta C_{i_1 j_1 k_1 l_1}(\mathbf{x}_1) \frac{\partial}{\partial x_{1l_1}} \left( G_{k_1 i_2}^0(\mathbf{x}_1, \mathbf{x}_2) \frac{\partial}{\partial x_{2j_2}} \left( \delta C_{i_2 j_2 k_2 l_2}(\mathbf{x}_2) \delta C_{i_3 j_3 k_3 l_3}(\mathbf{x}_3) \frac{\partial^2 G_{k_2 i_3}^0(\mathbf{x}_2, \mathbf{x}_3)}{\partial x_{2l_2} \partial x_{3j_3}} \frac{\partial G_{k_3 n}^0(\mathbf{x}_3, \mathbf{x}')}{\partial x_{3l_3}} \right) \right) \right) d\mathbf{x}_1 d\mathbf{x}_2 d\mathbf{x}_3 \\
&= - \iint G_{mi_1}^0(\mathbf{x}, \mathbf{x}_1) \frac{\partial}{\partial x_{1j_1}} \left( \delta C_{i_1 j_1 k_1 l_1}(\mathbf{x}_1) \frac{\partial}{\partial x_{1l_1}} \left( \frac{\partial G_{k_1 i_2}^0(\mathbf{x}_1, \mathbf{x}_2)}{\partial x_{2j_2}} \delta C_{i_2 j_2 k_2 l_2}(\mathbf{x}_2) \delta C_{i_3 j_3 k_3 l_3}(\mathbf{x}_3) \frac{\partial^2 G_{k_2 i_3}^0(\mathbf{x}_2, \mathbf{x}_3)}{\partial x_{2l_2} \partial x_{3j_3}} \frac{\partial G_{k_3 n}^0(\mathbf{x}_3, \mathbf{x}')}{\partial x_{3l_3}} \right) \right) d\mathbf{x}_1 d\mathbf{x}_2 d\mathbf{x}_3 \\
&= - \iint G_{mi_1}^0(\mathbf{x}, \mathbf{x}_1) \frac{\partial}{\partial x_{1j_1}} \left( \delta C_{i_1 j_1 k_1 l_1}(\mathbf{x}_1) \delta C_{i_2 j_2 k_2 l_2}(\mathbf{x}_2) \delta C_{i_3 j_3 k_3 l_3}(\mathbf{x}_3) \frac{\partial^2 G_{k_1 i_2}^0(\mathbf{x}_1, \mathbf{x}_2)}{\partial x_{1l_1} \partial x_{2j_2}} \frac{\partial^2 G_{k_2 i_3}^0(\mathbf{x}_2, \mathbf{x}_3)}{\partial x_{2l_2} \partial x_{3j_3}} \frac{\partial G_{k_3 n}^0(\mathbf{x}_3, \mathbf{x}')}{\partial x_{3l_3}} \right) d\mathbf{x}_1 d\mathbf{x}_2 d\mathbf{x}_3
\end{align*}
$$

• $G^{(4)}(\mathbf{x}, \mathbf{x}')$:

$$
\begin{align}
G^{(4)}(\mathbf{x}, \mathbf{x}') &= \iiint\int G^0(\mathbf{x} - \mathbf{x}_1) \cdot \mathcal{L}^1(\mathbf{x}_1) (G^0(\mathbf{x}_1 - \mathbf{x}_2) \cdot \mathcal{L}^1(\mathbf{x}_2) (G^0(\mathbf{x}_2 - \mathbf{x}_3) \cdot \nonumber \\
&\qquad \mathcal{L}^1(\mathbf{x}_3) (G^0(\mathbf{x}_3 - \mathbf{x}_4) \cdot \mathcal{L}^1(\mathbf{x}_4) G^0(\mathbf{x}_4 - \mathbf{x}')))) \, d\mathbf{x}_1 d\mathbf{x}_2 d\mathbf{x}_3 d\mathbf{x}_4 \nonumber \\
&= - \iiint\int G^0(\mathbf{x}, \mathbf{x}_1) \cdot \nabla_{\mathbf{x}_1} \cdot (\delta C(\mathbf{x}_1) \otimes \delta C(\mathbf{x}_2) \otimes \delta C(\mathbf{x}_3) \otimes \delta C(\mathbf{x}_4) : \nonumber \\
&\qquad (\nabla^2_{\mathbf{x}_1 \mathbf{x}_2} G^0(\mathbf{x}_1, \mathbf{x}_2) \otimes \nabla^2_{\mathbf{x}_2 \mathbf{x}_3} G^0(\mathbf{x}_2, \mathbf{x}_3) \otimes \nabla^2_{\mathbf{x}_3 \mathbf{x}_4} G^0(\mathbf{x}_3, \mathbf{x}_4) \otimes \nabla_{\mathbf{x}_4} G^0(\mathbf{x}_4, \mathbf{x}'))) \nonumber \\
&\qquad d\mathbf{x}_1 d\mathbf{x}_2 d\mathbf{x}_3 d\mathbf{x}_4, \tag{61}
\end{align}
$$
Or with indicial notations:

$$
\begin{align}
G_{mn}^{(4)}(\mathbf{x}, \mathbf{x}') &= -\iiint\int G_{mi_1}^0(\mathbf{x}, \mathbf{x}_1) \frac{\partial}{\partial x_{1j_1}} \left( \delta C_{i_1 j_1 k_1 l_1}(\mathbf{x}_1) \frac{\partial}{\partial x_{1l_1}} \left( G_{k_1 i_2}^0(\mathbf{x}_1, \mathbf{x}_2) \frac{\partial}{\partial x_{2j_2}} \left( \delta C_{i_2 j_2 k_2 l_2}(\mathbf{x}_2) \frac{\partial}{\partial x_{2l_2}} \right. \right. \right. \right. \nonumber \\
&\qquad \left. \left. \left( G_{k_2 i_3}^0(\mathbf{x}_2, \mathbf{x}_3) \frac{\partial}{\partial x_{3j_3}} \left( \delta C_{i_3 j_3 k_3 l_3}(\mathbf{x}_3) \frac{\partial}{\partial x_{3l_3}} \right. \right. \right. \right. \right) \left. \left. \left( G_{k_3 i_4}^0(\mathbf{x}_3, \mathbf{x}_4) \frac{\partial}{\partial x_{4j_4}} \left( \delta C_{i_4 j_4 k_4 l_4}(\mathbf{x}_4) \frac{\partial G_{k_4 n}^0(\mathbf{x}_4, \mathbf{x}')}{\partial x_{4l_4}} \right) \right) \right) \right) \right) \mathrm{d}\mathbf{x}_1 \mathrm{d}\mathbf{x}_2 \mathrm{d}\mathbf{x}_3 \mathrm{d}\mathbf{x}_4 \nonumber \\
&= -\iint G_{mi_1}^0(\mathbf{x}, \mathbf{x}_1) \frac{\partial}{\partial x_{1j_1}} \left( \delta C_{i_1 j_1 k_1 l_1}(\mathbf{x}_1) \delta C_{i_2 j_2 k_2 l_2}(\mathbf{x}_2) \delta C_{i_3 j_3 k_3 l_3}(\mathbf{x}_3) \delta C_{i_4 j_4 k_4 l_4}(\mathbf{x}_4) \right. \nonumber \\
&\qquad \left. \frac{\partial^2 G_{k_1 i_2}^0(\mathbf{x}_1, \mathbf{x}_2)}{\partial x_{1l_1} \partial x_{2j_2}} \frac{\partial^2 G_{k_2 i_3}^0(\mathbf{x}_2, \mathbf{x}_3)}{\partial x_{2l_2} \partial x_{3j_3}} \frac{\partial^2 G_{k_3 i_4}^0(\mathbf{x}_3, \mathbf{x}_4)}{\partial x_{3l_3} \partial x_{4j_4}} \frac{\partial G_{k_4 n}^0(\mathbf{x}_4, \mathbf{x}')}{\partial x_{4l_4}} \right) \mathrm{d}\mathbf{x}_1 \mathrm{d}\mathbf{x}_2 \mathrm{d}\mathbf{x}_3 \mathrm{d}\mathbf{x}_4 \tag{62}
\end{align}
$$

• $G^{(p)}(\mathbf{x}, \mathbf{x}')$:

$$
\begin{align}
& G^{(p)}(\mathbf{x}, \mathbf{x}') \nonumber \\
&= (-1)^p \underbrace{\int \cdots \int}_{p} G^0(\mathbf{x}, \mathbf{x}_1) \mathcal{L}^1(\mathbf{x}_1) (G^0(\mathbf{x}_1, \mathbf{x}_2) \mathcal{L}^1(\mathbf{x}_2) \cdots (G^0(\mathbf{x}_{p-1} - \mathbf{x}_p) \mathcal{L}^1(\mathbf{x}_p) G^0(\mathbf{x}_p - \mathbf{x}'))) \mathrm{d}\mathbf{x}_1 \cdots \mathrm{d}\mathbf{x}_p \nonumber \\
&= - \underbrace{\int \cdots \int}_{p} G^0(\mathbf{x}, \mathbf{x}_1) \nabla_{\mathbf{x}_1} \cdot (\delta\mathcal{C}(\mathbf{x}_1) \otimes \cdots \otimes \delta\mathcal{C}(\mathbf{x}_p)) : \nonumber \\
& \qquad (\nabla^2_{\mathbf{x}_1 \mathbf{x}_2} G^0(\mathbf{x}_1, \mathbf{x}_2) \otimes \cdots \otimes \nabla^2_{\mathbf{x}_{p-1} \mathbf{x}_p} G^0(\mathbf{x}_{p-1}, \mathbf{x}_p) \otimes \nabla_{\mathbf{x}_p} G^0(\mathbf{x}_p, \mathbf{x}')) \mathrm{d}\mathbf{x}_1 \cdots \mathrm{d}\mathbf{x}_p, \tag{63}
\end{align}
$$

Or with indicial notations:

$$
\begin{align}
G_{mn}^{(p)}(\mathbf{x}, \mathbf{x}')
&= (-1)^p \iiint\int\int G_{mi_1}^0(\mathbf{x}, \mathbf{x}_1) \frac{\partial}{\partial x_{1j_1}} \left( \delta C_{i_1 j_1 k_1 l_1}(\mathbf{x}_1) \frac{\partial}{\partial x_{1l_1}} \left( G_{k_1 i_2}^0(\mathbf{x}_1, \mathbf{x}_2) \frac{\partial}{\partial x_{2j_2}} \left( \delta C_{i_2 j_2 k_2 l_2}(\mathbf{x}_2) \frac{\partial}{\partial x_{2l_2}} \right. \right. \right. \nonumber \\
& \qquad \left. \left. \left( G_{k_{p-1} i_p}^0(\mathbf{x}_{p-1}, \mathbf{x}_p) \frac{\partial}{\partial x_{pj_p}} \left( \delta C_{i_p j_p k_p l_p}(\mathbf{x}_p) \frac{\partial G_{k_p n}^0(\mathbf{x}_p, \mathbf{x}')}{\partial x_{pl_p}} \right) \right) \right) \right) \mathrm{d}\mathbf{x}_1 \cdots \mathrm{d}\mathbf{x}_p \nonumber \\
&= - \iint G_{mi_1}^0(\mathbf{x}, \mathbf{x}_1) \frac{\partial}{\partial x_{1j_1}} \left( \delta C_{i_1 j_1 k_1 l_1}(\mathbf{x}_1) \cdots \delta C_{i_p j_p k_p l_p}(\mathbf{x}_p) \right. \nonumber \\
& \qquad \left. \frac{\partial^2 G_{k_1 i_2}^0(\mathbf{x}_1, \mathbf{x}_2)}{\partial x_{1l_1} \partial x_{2j_2}} \cdots \frac{\partial^2 G_{k_{p-1} i_p}^0(\mathbf{x}_{p-1}, \mathbf{x}_p)}{\partial x_{p-1 l_{p-1}} \partial x_{pj_p}} \frac{\partial G_{k_p n}^0(\mathbf{x}_p, \mathbf{x}')}{\partial x_{pl_p}} \right) \mathrm{d}\mathbf{x}_1 \cdots \mathrm{d}\mathbf{x}_p \tag{64}
\end{align}
$$

## Appendix A2: Development of analytical formulas for phase velocity

One can define the inner integral in the complex plane as:

$$ I_C = \frac{1}{2} \int_{-\infty}^{\infty} \frac{z^4}{(1 - z^2)(1 + C_0 - zC_1 + z^2 C_2)^2} dz. \qquad (65) $$

The first thing that can be noticed is that the function have poles at $z = 1$ and at the roots of the

polynomial $1 + C_0 - zC_1 + z^2 C_2$, which are second order poles. Also, the integrand is convergent

at $z \to \infty$. For this specific type of integrals one can use an upper half plane contour (UHP) in

the real axis. After several calculations, it was observed that the correct solution was obtained for

the following form of the integral:

$$ I_C = \frac{1}{2C_2^2} \int_{-\infty}^{\infty} \frac{z^4}{(1 - z^2) \left( \tilde{C}_0 / C_2 - zC_1 / C_2 + z^2 \right)^2} dz, \qquad (66) $$
782 where $\tilde{C}_0 = 1 + C_0$. For this contour one can use the residue theorem:

$$I_C = \oint_C f(z) dz = 2\pi i \sum_k^n \operatorname{Res}_{z=z_k} f(z=z_k), \qquad (67)$$

783 where $z_k$ are the poles inside the contour. The first pole is a simple pole $z_1 = 1$. The residue for a simple pole in the real axis can be shown to be $\operatorname{Res}_{f(z=1)} = \pi/2$. The residue for the second order pole is calculated as:

$$\operatorname{Res}_{z=z_k} f(z_k) = \frac{1}{(m-1)!} \frac{\mathrm{d}^{m-1}}{\mathrm{d}z^{m-1}} [(z-z_k)^m f(z)] \Big|_{z=z_k} \qquad (68)$$

786 Before proceeding with the derivation of the function one can simplify it by writing the polynomial in terms of the complex poles as:

$$\frac{z^4}{(1-z^2)(1+C_0-zC_1+z^2C_2)^2} = \frac{z^4}{(1-z^2)(z-z_2)^2(z-z_2)^2}, \qquad (69)$$

788 where $\bar{z}_2$ corresponds to the complex conjugate of the second order pole. The derivative then can be found to be:

$$\begin{aligned} \frac{\mathrm{d}}{\mathrm{d}z} \frac{(z-z_2)^2 z^4}{(1-z^2)(z-z_2)^2(z-\bar{z}_2)^2} &= \frac{\mathrm{d}}{\mathrm{d}z} \frac{z^4}{(1-z^2)(z-\bar{z}_2)^2} \\ &= \frac{2z^5}{(1-z^2)^2(z-\bar{z}_2)^2} - \frac{2z^4}{(1-z^2)(z-\bar{z}_2)^4} + \frac{4z^3}{(1-z^2)(z-\bar{z}_2)^2} \end{aligned} \qquad (70)$$

790 The complex pole is:

$$\begin{aligned} z_2 &= \frac{C_1 \pm iD}{2C_2}, \\ D &= \sqrt{4\tilde{C}_0 C_2 - C_1^2}. \end{aligned} \qquad (71)$$

791 After manipulation, the solution can be finally written for $\chi$ (see Equation (52)), and with this 792 solution, the real part of the mass operator can be obtained (see Equation (51)).