Finite-element and semi-analytical study of elastic wave propagation in strongly scattering polycrystals

This work studies scattering-induced elastic wave attenuation and phase velocity variation in three-dimensional untextured cubic polycrystals with statistically equiaxed grains using the theoretical second-order approximation (SOA) and Born approximation models and the grain-scale finite-element (FE) model, pushing the boundary towards strongly scattering materials. The results for materials with Zener anisotropy indices A > 1 show a good agreement between the theoretical and FE models in the transition and stochastic regions. In the Rayleigh regime, the agreement is reasonable for common structural materials with 1 < A <  3.2 but it deteriorates as A increases. The wavefields and signals from FE modelling show the emergence of very strong scattering at low frequencies for strongly scattering materials that cannot be fully accounted for by the theoretical models. To account for such strong scattering at A > 1, a semi-analytical model is proposed by iterating the far-field Born approximation and optimizing the iterative coefficient. The proposed model agrees remarkably well with the FE model across all studied materials with greatly differing microstructures; the model validity also extends to the quasi-static velocity limit. For polycrystals with A < 1, it is found that the agreement between the SOA and FE results is excellent for all studied materials and the correction of the model is not needed.


Introduction
Elastic waves scatter as they travel through inhomogeneous media and thus exhibit scatteringinduced attenuation and phase velocity dispersion. The problem of wave propagation and scattering in inhomogeneous media has received extensive study in the fields of e.g. seismology [1,2] and non-destructive evaluation [3,4]. A subject of particular interest in both fields is to study wave propagation in polycrystalline media to facilitate the detection and characterization of inhomogeneities like faults, defects and grains within the media. The first complete theoretical treatment of the scattering-induced attenuation and velocity dispersion in polycrystals was conducted by Stanke & Kino [5] based on the multiple scattering theory developed by Karal and Keller [6,7]. This model was formulated for untextured polycrystals with statistically equiaxed grains of cubic symmetry. An equivalent model was later offered by Weaver [8] using the Dyson equation and the first-order smoothing approximation [9][10][11]. A variety of extensions of this later model have since been performed for various grain structures and crystal symmetries, see e.g. [12][13][14][15] or our earlier work [16][17][18][19] for an overview.
The theoretical models of the Stanke & Kino [5] and Weaver [8] type take the statistical information of the polycrystals as input and for a given wave modality they produce the scattering-induced attenuation coefficient and phase velocity as outputs. The models (before the Weaver model [8] invokes the Born approximation) are accurate when the second-order degree of inhomogeneity is small and hence we collectively call both models the secondorder approximation (SOA) following Stanke & Kino [5]. The SOA results exhibit three specific frequency regions that are known as the Rayleigh, stochastic and geometric regimes with increased scattering intensities. All three regimes can also be predicted if the far-field approximation (FFA) [20] is invoked in the SOA model but the strongly scattering geometric regime vanishes if the single-scattering Born approximation is employed. The validity of these model approximations has recently been evaluated by three-dimensional grain-scale finiteelement (FE) simulations which are capable of accurately describing the interaction of waves with grains [16][17][18][19][21][22][23]. These comparative studies demonstrated that the SOA, FFA and Born models agree very well with the FE results in the simulated Rayleigh and stochastic regimes. The SOA and FFA models are mostly indistinguishable from each other and have a better agreement with the FE results than the Born model. These studies showed that the theoretical models are valid for polycrystals with the spatial two-point correlation (TPC) of either scalar type for statistically equiaxed grains [16,17,[21][22][23] or direction-dependent form for statistically elongated grains [18,19]. It was also shown that the elastic scattering factors [20], as combinations of elastic constants, are representative of the degree of inhomogeneity in the theoretical models for polycrystals of the highest cubic [16,18,19] and lowest triclinic [17][18][19] crystal symmetries; for cubic symmetry, the elastic scattering factors are related to the Zener anisotropy index A.
Despite the excellent performance of the SOA model, the comparative studies of the theoretical and FE models also revealed that, for cubic polycrystals with A > 1, the theoretical models start to deviate from the FE results at low frequencies as A increases [12][13][14][15]18], while the agreement remains good in the transition and stochastic regimes. This is somewhat unexpected because it is reasonable to assume that the theoretical models perform less satisfactorily at high frequencies rather than at low frequencies since the degree of scattering increases with frequency. This finding has led us to hypothesize that strong scattering arises at low frequencies in strongly scattering polycrystals that is not fully considered by the theoretical models since they only account for a subset of scattering events. This work aims to further investigate this finding by studying a variety of equiaxed materials of cubic symmetry with greatly differing A; some studied materials have significantly larger A than previously considered [5,[16][17][18][19][21][22][23]. We shall see that the validity

Finite-element model
We use the three-dimensional FE method to simulate the propagation of plane longitudinal waves in polycrystals with statistically equiaxed grains. We have reported the details of this method in our prior work [16,17,21,22,28] (see [18,19] for polycrystals with elongated grains) so only the essential steps are summarized below.
In the three-dimensional FE method, we use a cuboid, composed of densely packed and fully bonded convex grains, to represent a polycrystal, as can be seen from figure 1a. The geometric models created for this work are detailed in table 1. The selection of model dimensions ensures that there are more than 10 grains and 10 wavelengths in the z-direction of wave propagation. Each model is deployed in three separate forms, involving three microstructures with the same number of grains per unit volume but different grain uniformities, as illustrated in figure 1b. The three microstructures are generated using the Neper program [29] by the Laguerre, Poisson Voronoi and centroidal Voronoi tessellations [29][30][31], as abbreviately called Laguerre, PVT and CVT hereafter. The PVT creates uniformly random seeds in the model space of a polycrystal, with each seed being enclosed by a convex grain within which all points are closer to the enclosed seed than to any other [16][17][18][19]21,22,28]. The equivalent spherical radii of the PVT grains are normally distributed as demonstrated in [21][22][23] and shown in figure 1c. In comparison to the PVT, the seeds of the Laguerre tessellation are weighted [30,31] to create the type of microstructure as commonly found in applications [32,33], in which case the equivalent grain radii follow the lognormal distribution as illustrated in figure 1c. The CVT iteratively changes the locations of the seeds to achieve a uniform distribution of grains and the resulting equivalent grain radii follow a much narrower normal distribution than that of the PVT [23], as can be seen from figure 1c. As discussed below and shown in figure 1d, the three microstructures have greatly differing TPC functions and thus induce distinctive scattering behaviours. Each generated polycrystal model is then discretized in space and time. The spatial discretization uses a structured mesh to divide the model space into identical eight-node linear 'brick' elements as illustrated in figure   and the size is chosen to be smaller than λ/10 and D/10 to satisfactorily suppress the numerical error of the simulation and to well represent the microstructure of the model [21,22,28], where λ is the wavelength and D is the cubic root of the average grain volume representing the average grain size. The time-stepping solution is based on the central difference scheme [28,34] using a time step of t = 0.8h/V 0L that satisfies the Courant-Friedrichs-Lewy condition, where V 0L is the longitudinal Voigt velocity. Each polycrystal model is subsequently assigned material properties. We consider singlephase polycrystals in this work, so the individual grains of each model have the same mass density and elastic constants. The grains within each model are defined with uniformly randomly oriented crystallographic axes, making the model macroscopically homogeneous and isotropic (untextured) [21]. The polycrystalline materials used in this work are recorded in table 2. The materials have the same cubic crystal symmetry but different anisotropy indices A = 2c 44 /(c 11 − c 12 ), with eight materials having A > 1 (four naturally occurring and four fictitious [16]) and four having A < 1 [35].
We simulate the propagation of plane longitudinal waves in the z-direction.  Table 2. Polycrystalline materials with cubic crystal symmetry. Elastic constants c ij (GPa), density ρ (kg m −3 ), Voigt velocities V 0L/T (m/s), Zener anisotropy index A and elastic scattering factors Q L→T/L .   Together with the SBCs on the y = 0 and y = d y surfaces, an infinitely wide model would be formed across the x-y plane to sustain plane longitudinal waves. We note that the SBC is not an exact representation for the behaviour at an arbitrarily located slice through a polycrystal, but it is the best option we have available. The alternative periodic boundary condition could be used in this case, but it substantially increases the complexity of the FE solution, and in any case it too is approximate. Further discussion is given in [22,28]. To excite the desired z-direction plane longitudinal wave, a z-direction force in the form of a three-cycle Hann-windowed toneburst is uniformly applied to every node on the z = 0 surface. Each polycrystal model is solved in the time domain using the GPU-accelerated Pogo program [34]. z-direction displacements are monitored during the time-stepping solution, and example results are provided in figure 2 for a single realization of the model N115200 with the PVT microstructure, simulated at a centre frequency of 1 MHz (2k 0L a ≈ 1). Figure 2a shows the displacement wavefields on an arbitrary cross-section at an arbitrary time. Figure 2b displays the signals at individual nodes as thin grey lines and shows the respective coherent signals as thick lines that are averaged over all nodes on the monitoring boundaries. Figure 2c presents the normalized wavefront fluctuations on the receiving boundary at a normalized frequency of 2k 0L a = 1. To find the value of such a fluctuation, the time domain displacement field at z = d z , namely u z (x, y; t), is Fourier transformed to the frequency domain, u z (x, y; f ). The amplitude of the resulting field is normalized by the coherent amplitude to get the fluctuation as u f (x, y; f ) = |u z (x, y; f )|/ |u z (x, y;f )| x,y − 1, where · x,y represents the coherent average over all x and y nodes. The fluctuations can be alternatively given in the time domain but they are not provided here because they essentially offer the same information as the frequency domain ones, see [28] for details. Since the fluctuation is normalized by the coherent amplitude, its root-mean-square royalsocietypublishing.org/journal/rspa Proc. R. Soc. A 478: All figure panels show stronger scattering as A increases from aluminium through Inconel to lithium, leading to weaker coherent signals on the receiving boundary. However, the recorded coherent signals are strong and clear even for lithium. The coherent waves are only affected by the multiple scattering events arriving simultaneously with the main coherent beam. In addition, the incoherent waves are diminished after averaging over a significant number (approx. 60 000 in the least case) of monitoring nodes in each case; this was discussed in-depth in Sec. IV.3 of [19]. Therefore, the coherent signals have a high degree of signal-to-noise ratio (SNR) and the SNR is further improved in this work by averaging over as many as 15 repeated realizations of the model domain as will be discussed below. Such coherent signals are well suited to calculating scattering-induced attenuation and velocity change. This is achieved by applying appropriate time windowing to the emitted and received coherent signals and then comparing the resulting main wave pulses in the frequency domain; see [21,22,28] for details. The SNR would decline as frequency increases towards the transition to the geometric regime, but all FE results reported in this work have achieved a good SNR by limiting them to a reasonable frequency range [19]. We note that we use all three models in table 1 for the simulation of each material in table 2. The individual models are simulated with different centre frequencies for the applied toneburst force, allowing each material to be simulated in a wide frequency range as reported later in this work; example frequency ranges covered by individual models are provided in table 1 for aluminium. The combination of different polycrystal models with various centre frequencies delivers a very high degree of numerical accuracy across the whole-frequency range [28]. Besides, the simulation results presented in this work have achieved a very good statistical convergence by taking the average of 15 realizations for each case (each model-frequency combination). The different realizations of each case use the same polycrystal model but the crystallographic orientations of the grains are randomly reshuffled for each realization [22]. Alternatively, these multiple realizations can also be based on polycrystal models of different grain arrangements [22]. It was shown [22] that this approach is equivalent to, but more computationally intensive than, that used here. The applicability of our approach is also supported by the FE results presented in §4 and §5 that show consistency across different models at their overlapping frequencies for all studied materials of different microstructures.
The TPC statistics w(r), describing the probability of two points separated by a distance of r falling into the same grain, are numerically measured from the generated polycrystal models and the resulting data points are shown in figure 1d for the three polycrystal microstructures. Here we treat the TPC as a scalar function because it is direction independent. Taking the model N11520 with the PVT microstructure as an example, the TPC curves measured in 30 randomly chosen directions have a mean correlation length (i.e. integral of TPC curve) of 0.23 mm and a very small standard deviation of 7.19 × 10 −4 mm. This supports the direction independence of the TPC and also substantiates the statistically equiaxed nature of the grains. Since all three polycrystal models in table 1 have a large number of grains, they possess nearly the same statistical characteristics for a given microstructure. Taking the PVT microstructure as an example, the TPCs of the three models are indistinguishable [28] and their average correlation length is 0.23 mm, with a standard deviation of 5.84 × 10 −5 mm. For this reason, we have selected the model N11520 to determine the TPC for each microstructure. We note that the SBCs effect, which doubles the sizes of the symmetry boundary grains and makes these grains larger on average than the grains located within the body of the models [28], is not considered when measuring the TPC; this will be further discussed in §4a. To incorporate the measured statistics into the theoretical models, they are fitted into generalized TPC functions, w(r) = i A i e −r/a i , which are displayed in figure 1d as solid lines; the A i and a i coefficients are provided in the electronic supplementary material, table S1 for the three microstructures. Detailed TPC measurement and fitting procedures are reported in [16,17,28]. As can be found in figure 1d, the agreement of the generalized TPC curves with the points is less satisfactory at the tail than at the origin. Improving the agreement at the tail would increase the accuracy of the volumetric characteristic of the polycrystal, as described by the effective grain volume V [5,8,17] in the electronic supplementary material, table S2. However, the accuracy improvement is very limited and is thus not pursued here because the TPC at the tail is very small in probability. In the electronic supplementary material, table S2, we also provide the other two characteristic parameters of the TPC as will be used below, including the mean line intercept a = −1/w (r = 0) = 1/ i (A i /a i ) [5,17,36] and correlation length a CL = ∞ 0 w(r) dr = i A i a i [5,37,38]. Note that the three microstructures have approximately the same mean line intercept (i.e. the same slope at the origin).

Theoretical models
As with the above FE method, the same wave propagation problem of elastic wave propagation in polycrystals with statistically equiaxed grains is addressed here from a theoretical perspective. The theoretical models considered here use the statistical information of the FE models as input, enabling a direct comparison of both methods. These models are briefly introduced below; readers are referred to [17] for details.

(a) Second-order approximation
The equiaxed polycrystals considered in this work are macroscopically homogeneous and isotropic; the spatially varying elastic tensor of a polycrystal can thus be expressed as c ijkl (x) = c 0 ijkl + δc ijkl (x), with the Voigt average c 0 ijkl = c ijkl (x) representing the homogeneous reference medium and δc ijkl (x) denoting the elastic fluctuation. An incident wave scatters on the elastic fluctuation and the wavenumber is thus perturbed as the wave propagates. The perturbed wavenumber k satisfies the dispersion equation [8,16,17,20] where k = kp is the wavevector; the unit vector p represents the wave propagation direction, which is arbitrary due to the macroscopic isotropy of the polycrystal. ω = 2π f is the angular frequency and f is the frequency. V 0M represents the Voigt phase velocity of the wave M in the homogeneous reference medium. m M = N=L,T m M→N is the spatial Fourier transform of the mass operator describing the random scattering events occurring in the polycrystal. The component m M→N , denoting the scattering of the wave M into N, is given below [16,17,20] by using the first-order smoothing approximation [8,9] (equivalent to the Bourret approximation [9][10][11]) where the mass density ρ is constant for a single-phase polycrystal considered in this work. k 0N denotes the wavenumber of the wave N in the reference medium. P.V. represents the Cauchy principal value and ξ is a dimensionless variable. The coefficient η is 1 and 2 for longitudinal (M = L) and transverse (M = T) propagating waves, respectively. The factor f M→N in equation (3.2) describes the TPC of the elastic fluctuation and is given by [16,17,20] where the terms in the summation symbol correspond to the spectral representation of the spatial TPC function, w(r) = i A i e −r/a i , and the rest of the terms in the parentheses represent the elastic part of the TPC. For longitudinal waves in cubic polycrystals, the A MN , . . . . coefficients in equation 44 is the invariant anisotropy coefficient [16] (c = 0; A = 1 for isotropy); those coefficients for arbitrary crystal symmetries can be found in our prior work [17]. One obtains the perturbed wavenumber k for the propagating wave M by numerically solving the dispersion equation, equation (3.1). Consequently, the attenuation coefficient and phase velocity are calculated from the imaginary and real parts of the perturbed wavenumber by α M = Imk and V M = ω/Rek.

(b) Far-field approximation
The FFA model does not involve the complex calculation of the Cauchy integral as in the SOA model, and it has an explicit expression for the mass operator as given below for polycrystals with statistically equiaxed grains [17,20] .  [17,20], with V R M being the Rayleigh velocity limit given below. The solution of the FFA model is mostly indistinguishable from that of the SOA model across the wholefrequency range but the former is more computationally efficient [17]. An important advantage of the FFA model for our purpose is that it explicitly relates the attenuation and velocity dispersion to the elastic scattering factors Q M→N [20]. Two factors, [20], exist for macroscopically isotropic polycrystals with grains of arbitrary symmetry, and they have simple expressions of Q L→L = 4c 2 /(525 c 11 2 ) and Q L→T = c 2 /(105 c 11 c 44 ) in a cubic polycrystal [20].

(c) Born approximation
Closed-form solutions can be found for the SOA and FFA models by invoking the Born approximation. Since the resulting numerical solutions are nearly the same [17], here we only provide the Born approximation of the FFA model. To obtain the solution, we substitute (ω/V 0M ) 2 − k 2 by 2k 0M (k 0M − k) in equation (3.1) and replace k with pk 0M in equation (3.4). This leads to the following solution for a longitudinal propagating wave [17,20] where 2(Q * LL + Q L→T )k 0L is added to consider the above-mentioned velocity correction. The imaginary and real parts of the solution k L are further obtained as and Rek where

(d) Rayleigh asymptotes
At the low-frequency Rayleigh limit, the attenuation and phase velocity asymptotes are given by [17] where N = M. V g eff is the effective grain volume defined by the volumetric integral of the TPC function [5,8,17], electronic supplementary material, table S2.
is an elastic factor for simplifying the equation. For longitudinal waves, since Q * LL is generally negligible and V 3 0L /V 3 0T just differs slightly among most structural materials, we can obtain from equation (3.8) that α R L ∝ Q L→T . Appending the fact that Q * LL and Q L→T are far smaller than unity, it follows from equation (3.8) that the phase velocity variation is

(e) Stochastic asymptotes
At the high-frequency stochastic limit, the attenuation and phase velocity asymptotes are given by [17] , (3.9) in which N = M. a CL is the correlation length defined by the integral of the TPC function [37], is an elastic factor introduced for simplifying the equation. For longitudinal waves, it follows from equation (3.9) Thus, Q L→L is the major elastic factor determining the level of scattering at the stochastic limit.

Comparison of finite-element and theoretical models
Now we present and discuss the numerical FE model (FEM) results and theoretical predictions. The first three subsections focus on the eight cubic materials with anisotropy indices A > 1, while the last subsection on the four materials with A < 1. Both cases consider the PVT microstructure only.
(a) Dependence of attenuation and velocity on frequency The FEM points of each material are obtained using all three polycrystal models in table 1, modelled with multiple centre frequencies, to cover the shown frequency range. A combination of 15 realizations is used for each modelling case to achieve statistically meaningful results. The average of the multiple realizations is shown as the points in the figure while the corresponding standard deviation (i.e. error bar) is not provided because it is smaller than the size of the FEM point markers. The relative differences between the theoretical curves and the FEM points are provided in figure 4. A small discontinuity (sudden jump) can be observed as we look along sequential FEM points for each material, which is more evident from the phase velocity points for A = 5.0 and lithium at around 2k 0L a = 1 (those jumps are especially visible for the relative differences in figure 4). Such a discontinuity occurs because two different FE material models are used for calculating its left-and right-side FEM points; e.g. N115200 and N11520 for the left and right sides of 2k 0L a ≈ 1. The models on the two sides have different mesh sizes and thus different numbers of elements per wavelength at their overlapping frequency of calculation. This causes different numerical errors [28] that exhibit a discontinuity of the FEM results. The discontinuity is observable in the phase velocity results because the FE scheme used in this work is more prone to numerical phase errors, namely numerical dispersion [28]. Also, this discontinuity is more evident for highly scattering materials because numerical errors depend on material anisotropy [28]. The discontinuities essentially define the bound of the numerical error, which is mostly one order of magnitude smaller than the difference between the theoretical and numerical FE results, as can be more evidently seen from figure 4. For the extreme case, i.e. the phase velocity of lithium, the numerical error is about 0.4%, whereas the SOA-FEM difference is 2.4%. In essence, the FE results have achieved a very high degree of numerical accuracy and statistical convergence (also see our prior studies [21,22,33] and especially [28] for a more detailed assessment of these two aspects), and therefore they are used as reliable references below to evaluate the approximations of the theoretical models.
The theoretical SOA and Born curves are produced by incorporating the generalized TPC function of the FE material models, electronic supplementary material, table S1. The FFA model results are indistinguishable from the SOA curves [17] and are thus not provided. We   note that the SBCs effect is not considered in the measured and generalized TPC functions and is accordingly not accounted for in the theoretical SOA and Born curves. Since the SBCs enlarge the symmetry boundary grains and slightly increase the scattering intensity in the FE simulations (such an increase occurs to both the L-L and L-T scattering components but the L-T component may increase more at low frequencies where the L-T component is dominant) [28], the theoretical predictions would underestimate the level of attenuation and phase velocity variation in comparison to the FE results. Nonetheless, our prior work [28] has shown that this underestimation is very small and the same estimation for all materials addressed in this work reveals that this underestimation does not depend on material anisotropy. Also, it can be seen from figure 2a that   For attenuation, the theoretical results approach the Rayleigh and stochastic asymptotes in the low-and high-frequency regions for all materials, and therefore they have fourth-and secondpower dependences on frequency in the two regions, as follows from equations (3.8) and (3.9); the FEM results display the same frequency dependencies. The 'hump' in the transition frequency regime is a unique feature for longitudinal waves, which is attributed to the transition of the dominant L → T scattering in the Rayleigh regime to the dominant L → L in the stochastic regime [5,17,23]. Also, for a more strongly scattering material [5,17,19], the attenuation results show a less pronounced stochastic region and an earlier transition to the geometric regime (seemingly flat part appearing at the far right of the SOA curve for lithium). For phase velocity, the FEM and theoretical results tend to be non-dispersive in the Rayleigh regime while having different static limits for high anisotropies. Between the Rayleigh and stochastic asymptotes the phase velocity is dispersive due to transition of the dominant scattering mechanisms. The velocity results also show a shorter stochastic regime for a more anisotropic material, and this regime can hardly be observed from the FEM points because (i) for weakly scattering materials (aluminium, A = 1.5, A = 1.8 and A = 2.4), the velocity decreases slowly with frequency due to the numerical dispersion in the FE simulation [28], while (ii) for strongly scattering materials, the velocity increases with frequency potentially as a result of the numerical dispersion and the complete disappearance of the stochastic regime due to the strongly scattering nature of the materials.
Considering both the low-and high-frequency regimes, the attenuation and velocity variation ranges in figure 3 increase with the material anisotropy. This is due to the increase of scattering with material anisotropy, as can be observed from the wavefields, signals and wavefront fluctuations of different materials in Following the increase of scattering with material anisotropy, the theoretical SOA and Born curves in figure 3 show a deteriorated agreement with the FEM points as material anisotropy increases. As aforementioned, this is because the theoretical models only account for a subset of scattering events whereas the FEM points accurately incorporate all possible scattering. A further quantitative evaluation is provided in figure 4, showing the relative differences between A distinctive observation from figure 4 is that the relative differences in attenuation and phase velocity between the theoretical and FE models are relatively larger at low frequencies. Taking the frequency at 2k 0L a = 1 as an example, the relative difference in attenuation between the SOA and FEM results increases from −10% for aluminium to −66% for lithium, while the relative difference in phase velocity increases from 4 × 10 −4 % for aluminium to 2.4% for lithium. Such a large SOA-FEM difference is somewhat unexpected since it was previously believed that scattering is small in the low-frequency range and can thus be appropriately accounted for by the SOA model. However, it can be observed from figure 2 that very strong scattering arises at low frequencies, especially for strongly scattering lithium. Also, a strong scattering signal is visible in figure 2b long after the main pulse. Our interpretation here is that we are seeing multiple scattering after the arrival of the main coherent pulse; however, its effect on the main signal's amplitude and phase is unproven; it is possible that at low frequencies both could be solely affected by strong single scattering in the strongly anisotropic materials. One of the limitations of the SOA model is that it considers multiple scattering only partially, and the Born approximation accounts for only single-scattering events and thus deviates even more greatly from the FEM. However, we note that the Born-SOA difference is much less than the SOA-FEM difference; for attenuation at 2k 0L a = 1, the former is −0.3% for aluminium and −14.9% for lithium.
In contrast with the low-frequency range, the middle-frequency range exhibits mostly smaller differences, and interestingly, the attenuation differences nearly overlap across the materials in the range of 2k 0L a = 5 − 12 in figure 4a. At very high frequencies, the differences in both attenuation and phase velocity tend to grow with frequency and material anisotropy. It might be valuable to see how the differences progress with frequency in the future when computation power allows. Overall, the SOA model has a reasonable agreement with FEM in the middle-and high-frequency regions, and the single-scattering Born model agrees just as well with FEM even for strongly scattering polycrystals. This may suggest that multiple scattering is not strong even in these regions.

(b) Dependence of attenuation and velocity on elastic scattering factors
It has been demonstrated by figure 3 that attenuation and velocity variation increase with material anisotropy. A similar increase has also been found in figure 4 for the relative difference between the theoretical and FE models, revealing a more prominent anisotropy dependence at relatively low and high frequencies. For this reason, two normalized frequencies, 2k 0L a = 1 and 2k 0L a = 12, are chosen in these two ranges to quantitatively evaluate this dependence. At these two frequencies, the normalized attenuation and phase velocity points are extracted from figure 3 and plotted versus the elastic scattering factors in figure 5. Since the two frequencies roughly fall into the Rayleigh and stochastic regimes, the elastic scattering factors Q L→T and Q L→L are employed to characterize the degrees of scattering respectively, see §3 for the selection of these factors. We note that the factor ε 2 = (k L − k 0L ) 2 /k 2 0L defined in [5,39] is identical to the factor Q L→L (thus refer to the Q L→L column of table 2 for the values of 2 ) that describes the degree of inhomogeneity for the longitudinal-to-longitudinal scattering in the stochastic regime. We use the factor Q L→T to describe the degree of inhomogeneity for the longitudinal-to-transverse scattering dominating in the Rayleigh regime [17,19].
At the low frequency 2k 0L a = 1, the numerical FE points show a distinct quadratic relationship with the elastic scattering factor Q L→T for both attenuation and phase velocity. Quadratic fits are generated for the data points and are plotted in the figure, and   with the goodness-of-fit of R 2 = 0.997 and R 2 = 0.999, respectively. By comparison, the theoretical SOA and Born predictions have a rather different dependence of linear order on the scattering factor Q L→T , which is represented by the linear fits in the figure. Also, these predictions are smaller in magnitude than the FEM points due to the aforementioned reason that the theoretical models only consider a subset of scattering events whereas the FEM considers all. At the high frequency 2k 0L a = 12, all results suggest a linear dependence on the scattering factor Q L→L , as is more evident from the attenuation results shown in figure 5b. Nonetheless, we emphasize that this assertion of linear dependence may not be appropriate. This is due to the possibility that not all evaluated materials are consistently in the stochastic regime at 2k 0L a = 12 because weakly scattering materials are yet to enter the stochastic regime while the strongly scattering ones are already transiting into the geometric regime.
(c) Quasi-static velocity limit The quasi-static limit of longitudinal phase velocity is related to the effective elastic constant C 11 of the polycrystal medium by V L = √ C 11 /ρ. Since this limit can be determined to a high degree of accuracy using three-dimensional FEM [17], its FEM results can be used to evaluate the suitability of effective medium theories. For this reason, we calculate the quasi-static FEM velocities for the eight cubic materials with A > 1 and provide the results in figure 6a as a function of the scattering factor Q L→T . Each point is statistically converged by taking the average of 30 realizations using the model N115200 (table 1), and the respective error bar is not shown because it is smaller than the size of the markers. The effective medium theories considered here include the Hashin-Shtrikman (HS) bounds [40][41][42] and the self-consistent (SC) theory [43,44]. The lower and upper HS bounds, denoted, respectively, as LHS and UHS, prescribe the limiting range for the quasi-static velocity, while the SC theory provides a unique estimation of quasi-static velocity satisfying the continuity of stress and strain throughout the polycrystal. The bounds in figure 6a are calculated with the scripts by Brown [42] and the SC results with the scripts by Kube and De Jong [43]. Additionally, the Rayleigh asymptote of the SOA model, equation Consistently for all materials, the FEM, SC and SOA points lie well between the LHS and UHS bounds, whose range becomes wider as the scattering factor increases. The SOA estimates perform well for weakly scattering materials, but the agreement deteriorates with the increase of grain anisotropy. The SC theory shows an excellent agreement with the FEM results even for lithium, while the SOA model is less satisfactory in this regard. All results show a quadratic relationship between the normalized quasi-static velocity V L /V 0L − 1 and the elastic scattering factor Q L→T . This sole dependence on the elastic scattering factor is a new finding and is simpler than the previously observed dependence on both the Poisson's ratio and anisotropy index [16]. The differences between the theoretical and FEM results also exhibit a quadratic dependence on Q L→T .

(d) Cubic materials with anisotropy indices A < 1
Here, we present the results for the four cubic materials with anisotropy indices smaller than 1 (A < 1, table 2  In contrast with the A > 1 case, a prominent finding is that the theoretical SOA predictions agree excellently with the FEM points. The normalized root-mean-square deviations (RMSD) of the SOA results from the FEM results are 0.94%, 4.57%, 5.48% and 6.88% for RbF, RbCl, RbBr and RbI for attenuation, while the respective numbers for phase velocity are 0.15%, 0.36%, 0.45% and 0.51%. We note that RbI has approximately the same Q L→T as lithium and its universal anisotropy index [45] is of the same order of magnitude as that for lithium (2.64 for RbI and 8.70 for lithium). The two materials also show very similar levels of scattering, as can be observed from their FE wavefields, signals and wavefront fluctuations in the electronic supplementary material, figure S1. Despite these similarities, the SOA-FEM difference for RbI is an order of magnitude smaller than that for lithium. An in-depth study shows that the FEM results can be fitted to second-order polynomials of Q L→T for both the A < 1 and A > 1 cases, but the second-order term for the A < 1 case is negligible while that for the A > 1 case is comparable to its linear-order term. By contrast, the SOA model predicts a linear dependence on Q L→T for both cases. This contrasting result seems to be the reason why the SOA model performs well only for the A < 1 case while not for the A > 1 case.
Both A > 1 and A < 1 cases in figures 6 and 7e exhibit an exceptional agreement between the SC and the quasi-static FEM results. This fact further supports the generality of the above conclusions on the dependences of Q L→T in the low-frequency range. We further confirmed this finding by analysing an extra set of 441 cubic materials (153 with A < 1 and 288 with A > 1, using elastic constants from [46]) by comparing our quasi-static results with those obtained from the SC estimate [43] and the Reuss bound; the details of this study will be reported elsewhere. In particular, we found from the simple analytical expression of the Reuss bound [47] that the quasistatic velocity bound has both linear and quadratic Q L→T terms for either case; however, the quadratic term is very much smaller than the linear term for the A < 1 case, whereas the two terms are comparable for the A > 1 case. Therefore, there seem to be some general grounds for the difference between the A < 1 and A > 1 cases; however, it is not clear to us why, physically, the coefficient for the second-order Q L→T term is so small for A < 1.

Semi-analytical model for strongly scattering polycrystals with A > 1
The above attenuation and velocity results show a reasonably good agreement between the theoretical SOA and numerical FE models at high frequencies beyond the Rayleigh regime. In the low-frequency Rayleigh regime, however, the results exhibit a contrasting feature for the A < 1 and A > 1 cases. This is because the FEM results have different dependences on Q L→T , being nearly linear, at A < 1, and quadratic, at A > 1, whereas the SOA model exhibits a linear order for both cases. This leads to an excellent agreement between the SOA and FE models for the A < 1 case but a less satisfactory agreement for the A > 1 case that deteriorates rapidly with the increase of A. For this reason and the fact that most structural materials have anisotropy indices greater than unity, our focus of this section will be on the A > 1 case only.
The SOA model is inherently approximate; among those approximations the most important are: (1) The SOA model involves a major approximation by replacing a discrete polycrystal with a continuous random medium with fluctuating elastic tensor and statistical representation of the polycrystal by the TPC function [5,[8][9][10]. This replacement is intuitively applicable to materials of weak anisotropy but may introduce non-negligible errors for strongly anisotropic materials. (2) The SOA model uses the first-order smoothing approximation, as described by Weaver for polycrystals [8]. In the equivalent diagram perturbation series method, this relates to accounting for a subset of the scattering diagrams in the solution of the exact Dyson equation [9,10]. This approximation is also equivalent to the Keller approximation [7] as applied by Stanke & Kino [5] and in [6] to solids. The neglected scattering events may be negligible for weakly scattering materials but become increasingly important as material anisotropy gets stronger. (3) The model assumes the validity of factorizing the TPC function into the elastic and geometric terms [5,8]. There is some numerical support for the validity of this factorization [48]. (4) The high orders of the scattering diagrams depend on the multi-point correlation functions [9,10]. The effect of the additional statistics on scattering was not addressed in the literature.
The effect of those approximations on the obtained solution is not yet known even for the scalar case due to the lack of exact solutions. Therefore, numerical methods are the only alternative at this time to evaluate the quality of obtained solutions for polycrystals. The fact that the FEM results depend quadratically on the elastic scattering factor Q L→T indicates that an iterative approach may be applied to the theoretical SOA model (with a linear dependence on Q L→T ) to add a higher order term on the scattering factor. Following Rytov et al. [10] (pages 139-141), we may produce the iteration series for the SOA model by obtaining an initial effective wavenumber from the dispersion equation, then using it as the wavenumber for the reference medium to get from the dispersion equation the next iteration for the effective wavenumber. Further repeating this one obtains higher iteration solutions. Evidently, one iteration is sufficient to introduce a quadratic term of the elastic scattering factors into the solution. However, we note that even an infinite iteration will only consider a summation of a subset of the infinite types of scattering diagrams [10], and the contribution of the unaccounted diagrams seems significant as revealed from comparison with the FEM results. Thus, it is not feasible for this approach to fully take into account all scattering events. Also, note that the continuous random medium approximation (i) is done even before the Dyson equation can be derived by this approach. Alternatively, we propose a semi-analytical model by using the first iteration of the far-field Born approximation that results in a corrective second-order term on the scattering factor and then significantly increasing the coefficient of the corrective term for the model prediction to match the FEM results. Surprisingly, this empirical coefficient is nearly π 3 , and as will be seen below this semi-analytical model works very well for various cubic materials with different microstructures and also for cubic materials with elongated grains [49]. We start the iteration by assuming that the effective wavenumber k L of the Born approximation, equation (3.5), becomes the wavenumber of the reference medium [10] and thus substituting this wavenumber into equation (3.5) to form a new effective wavenumber. Since we are dealing with the low-frequency range where the L → T scattering is dominant, we use only the L → T term in equation (3.5) for the iteration. We do the iteration separately for the attenuation and phase velocity based on equations (3.6) and (3.7); this results in the expressions: and The terms 4p Im i Q L→T and 2p Re i Q L→T directly come from the iteration. The coefficients π 3 and π 3 /2 are obtained by best matching equations (5.1) and (5.2) with the FEM results at 2k 0L a = 1, equations (4.1) and (4.2). The iterative factors are given originally by p Im i = p Re i = 1/([1 + (k 0T a i ) 2 (η 2 LT − 1)] 2 + 4(k 0T a i ) 2 ), but our parametric study indicates that they need to be modified as follows to improve the transition of the semi-analytical model into the stochastic regime We note that the empirical coefficient π 3 /2 for Rek L (phase velocity) is half of that for attenuation, π 3 , and there is also an extra 1/2 constant in p Re i as compared to p Im i . This systematic difference between the correction coefficients for attenuation and phase velocity is not yet understood. It is worth pointing out that removing the summations over i in equations (5.1) and (5.2) would deliver a semi-analytical model for polycrystals with the TPC given by a single exponential, w(r) = e −r/a . Also, we emphasize that the empirical coefficients are nearly unchanged when the SBCs are accounted for because they only affect the TPC function and the effect is small, see §4a and [28].
At the Rayleigh limit, the attenuation and phase velocity asymptotes obtained from equations (5.1) and (5.2) are given in the especially simple form The stochastic attenuation and phase velocity asymptotes for the semi-analytical model are the same as those given by equation (3.9) because the model has retained the same stochastic behaviours as the Born approximation. Below, by comparing with the FEM results, we evaluate the applicability of the semi-analytical model, equations (5.1) and (5.2), to different materials and microstructures and also assess its accuracy for the quasi-static velocity limit, equation (5.4).

(a) Applicability of the semi-analytical model to cubic polycrystals with various anisotropy indices
First, we evaluate the applicability of the semi-analytical model to the eight cubic materials with the same PVT microstructure but greatly differing anisotropy indices. The semi-analytical model predictions are compared with the FEM and SOA results in figure 3. The figure shows that the semi-analytical model mostly overlaps with the FEM at low frequencies for both attenuation and phase velocity across all materials, demonstrating its greatly improved accuracy in comparison to the SOA model for high-grain anisotropies. Table 3 summarizes the normalized RMSD of the models from the FEM (the FEM values as the reference). It reveals that the semianalytical model mostly performs an order of magnitude better than the SOA model in the low-frequency range; the difference between the semi-analytical and FE models barely shows dependence on material anisotropy, especially for attenuation. Although not shown, we note that the Rayleigh asymptotes of the semi-analytical model also excellently represent the attenuation and velocity behaviours at the low-frequency limit. In the transition region, the semi-analytical model exhibits a slightly better agreement with the FEM results than the SOA model, and table 3 suggests that the bettering of the agreement is more evident for materials of a stronger scattering. The semi-analytical model approaches the SOA model in the stochastic range for all shown cases.

(b) Applicability of the semi-analytical model to different polycrystal microstructures
The excellent agreement between the semi-analytical model and the FEM discussed above are for the datasets for which the matching coefficients π 3 and π 3 /2 were determined in the model, equations (5.1) and (5.2). Nonetheless, it is also important to compare the model for unrelated microstructures. Here we evaluate the applicability of the semi-analytical model to different polycrystal microstructures with greatly contrasting TPC. Among the eight materials with A > 1, the above analysis illustrates that lithium most critically challenges the existing SOA and Born models. For this reason, the lithium polycrystal is used here for the evaluation, and it is additionally simulated with the Laguerre and CVT microstructures, §2. The resulting FEM points and those for the PVT microstructure (already shown in figure 3) are plotted in figure 8, compared with the predictions of the SOA and semi-analytical models. In contrast with the SOA model, the semi-analytical model has a remarkably better agreement with the FEM in both the low-frequency and transition regions for all three microstructures. In the lowfrequency region, in particular, the RMSD for the attenuation, given in table 3, decreases from 60% for the SOA model to 5% for the semi-analytical model, while that of phase velocity reduces from 2% to 0.3%. The figure also shows that the semi-analytical model practically overlaps with the SOA model in the stochastic range but starts to deviate from the latter because the model is developed based on the Born approximation. Essentially there is no difference in the model performance for those three microstructures, indicating the independence of the semianalytical model on the TPC of polycrystals; however, obviously, the TPC should be accurately measured.
(c) Applicability of the semi-analytical model to quasi-static velocity limit Finally, we appraise the applicability of the semi-analytical model to the quasi-static velocity limit, which may be of particular interest in developing effective medium theories. As shown in figure 9a,   achieving an order of magnitude improvement in accuracy in comparison to the SOA model. The normalized RMSD between the semi-analytical and FEM results over the eight materials is 0.04% whereas that between the SOA and FEM points is 0.81%. As shown in figure 9a, the SC estimate has the same excellent agreement with the FEM. Considering the additional results shown in figures 6 and 7e, it is reasonable to assume that the SC estimate would perform well for any cubic polycrystal. There are other reasons to believe that the SC method provides accurate estimates of homogenized elastic moduli of polycrystals: (i) iterative convergence of high-order bounds to the SC values [43] and (ii) FEM confirmation of the SC results [50]. For this reason, we use the SC estimate to calculate the quasi-static velocities for the same 288 materials with A > 1 (A = 9.81 is the largest) [46] as in §4d, and then we use the SC results as the reference to further evaluate the applicability of the semi-analytical model. It is found that the semi-analytical results nearly overlap with the SC estimations, figure 9b, with the normalized RMSD over the 288 materials being 0.04% (this is the same as that between the semi-analytical and FEM results in figure 9a). This further substantiates the validity of the semi-analytical model and also indicates the generality of the model constants, π 3 and π 3 /2, for cubic materials with A > 1. We note that the SC theory in general needs to be iteratively solved [43], whereas the semi-analytical model has a simple explicit expression given by equation (5.4).  [46]. The inset in (b) displays the respective relative difference, which is within ±0.2%; the normalized RMSD over all materials is 0.04%. Lines are linear/quadratic fits. (Online version in colour.)

Summary and conclusion
This work uses three-dimensional FE and theoretical models to study the scattering-induced attenuation and phase velocity variation of plane longitudinal waves in untextured cubic polycrystals with statistically equiaxed grains. The study is predominantly performed for materials with anisotropy indices greater than unity (A > 1). The results of such materials exhibit a good agreement between the SOA and FE models in the transition and stochastic regimes, even for very strongly scattering lithium (A = 9.14). This agreement also holds for the single-scattering Born approximation, thus indicating the possibility that the effect of multiple scattering on the coherent wave is weak in these regions. In the low-frequency Rayleigh regime, the theoretical models agree reasonably well with the FEM for common structural materials with A < 3.2: the largest difference in attenuation between the SOA and FEM is −10%, −35% and −37% for aluminium, Inconel and copper (the figures are slightly larger for the Born-FEM difference). However, the relative difference can reach the level of −70% for more strongly scattering materials like lithium. The emergence of such unsatisfactory agreement in the Rayleigh regime for cubic materials with A > 1 is somewhat unexpected. A study on materials with A < 1 is also conducted in the Rayleigh regime. It shows an excellent agreement between the SOA and FE models, with an attenuation difference of smaller than 7% for RbI that has nearly the same elastic scattering factor Q L→T as lithium. Further analysis reveals that this excellent agreement is due to the nearly linear dependence of the FEM results on Q L→T (proportional to the elastic covariance), which is the same as for the SOA results. By contrast, the FEM results for the A > 1 case are described by a quadratic polynomial on Q L→T with a significant dependence on the Q 2 L→T term, while the SOA model results are still linear to Q L→T (irrespective of A being larger or smaller than one). This contrasting dependence leads to the unsatisfactory SOA-FEM agreement for the A > 1 case. In addition to the FEM evidence, the linear and quadratic dependences of the A < 1 and A > 1 cases on Q L→T are supported by the quasi-static velocity limit, particularly by the SC estimate and Reuss bound. The SOA model is inherently approximate, and its disagreement with the FEM may be attributed to the replacement of the polycrystal by a continuous random medium and by approximations in the solution of the Dyson equation, mainly by limiting the order of perturbations, thus not accounting for all scattering events.
To consider strongly scattering materials with A > 1, we have proposed a semi-analytical model by iterating the far-field Born approximation and optimizing the coefficient of the secondorder term on the scattering factor Q L→T to achieve the best fit of the model to the FE results. We have demonstrated that the semi-analytical model works remarkably well for all materials considered in this work and for different polycrystal microstructures with largely differing TPCs. The largest difference in attenuation between the semi-analytical and FE models is within a reasonable ±15% range for all evaluated materials and microstructures. The semi-analytical model also delivers a very accurate prediction for the quasi-static velocity limit obtained by the FEM. In addition to the FE evidence supporting the semi-analytical model to predict the quasi-static velocity limit, an excellent agreement is observed between the semi-analytical model and the SC estimate for 288 materials with A > 1 (with a normalized RMSD of 0.04%, which is within the accuracy of the FE method). This finding substantiates the generality of the empirical semi-analytical model coefficients π 3 and π 3 /2 for cubic materials with A > 1.
The applicability of the proposed model is demonstrated for materials of cubic crystal symmetry, but we expect that the model may be applicable to polycrystals of lower symmetries and general inhomogeneous materials after the adjustment of the iterative coefficients. We hope that the simple form of the proposed semi-analytical model and its exceptional performance against the FE simulations will stimulate more rigorous theoretical developments.
Data accessibility. The FE simulations in this work were performed using the open-source program Neper (available at https://github.com/rquey/neper) for polycrystal model generation and meshing, and the proprietary software Pogo (demo version available at http://www.pogo.software/) for high-speed wave propagation solution. The simulation results presented in this work are provided in the electronic supplementary material [51].