Role of a fractal shape of the inclusions on acoustic attenuation in a nanocomposite

Heterogeneous architectured materials are of particular interest in engineering applications.¹ These are man-made structural materials that have been developed for obtaining ad hoc properties, which cannot be generally found in nature. They are obtained by engineering at different scales the mixing of different materials, either in a random spatial distribution (composite materials) or with the artificial repetition of regular patterns (metamaterials). Depending on their application, the length scale of such patterns can span from the nanometer to the macroscopic range, being smaller than, or comparable to, the wavelength of the phenomena that the material is meant to affect.

This concept has, indeed, been largely exploited for efficiently manipulating long-wavelength acoustic phonons, which assure sound propagation at low frequency, through the introduction in materials of periodic interfaces on a macroscopic scale [phononic crystal (PC)].^2–5 In recent years, they have attracted increasing attention among the scientific community due to their extraordinary acoustic and elastic wave propagation performances obtained by designing the phase gradient at the sub-wavelength scale, such as negative refraction, waveguiding, cloaking, and band gaps.^6–9 

Thanks to both engineering and theoretical progress, they have also been introduced in thermal science, with a microstructure in the nanometer scale. Thermal transport is intimately related to the sound propagation (acoustic transfer) in materials because in insulators and semiconductors, the main heat carriers are acoustic phonons.¹⁰ Specifically, at room temperature, heat is mainly carried by phonons with THz frequencies and nanometric wavelengths. As such, a periodic nanostructuration has proved to be promising for affecting phonon dispersions and their ability in transporting heat.¹¹ It is important to remind here that phonons participate to thermal transport through two different contributions: a propagative one, which depends on the phonon mean-free path, heat capacity, velocity, and vibrational density of states (DOS),¹² and a diffusive one, involving the phonon diffusivity rather than the mean-free path and velocity,^13–15 the latter dominating at higher energies and smaller wavelengths. The propagative contribution can be reduced through the presence of interfaces, which scatter or eventually trap phonons. This contribution is directly related to the coherent phonon transport, which can be understood on the basis of the wave nature of phonons. Much work has been recently done for investigating these coherent phenomena.^16–18 Depending on the relative importance of the diffusive contribution, in fact, the reduction but also the improvement of thermal conductivity have been reported in different nanostructured systems.^19–21 Numerous atomistic simulations of out-of-equilibrium phonon transport have also been reported in different systems.^22–30 Recent works have looked at the phonon dynamics to get a better insight into transport properties and found exotic behaviors such as an energy localization between pores,³¹ asymmetric transport (rectification),³² or the filtering of high-frequency phonons.³³ 

It is, thus, clear that, depending on the length scale at play, the design of the phononic crystal structure contributes directly to the performance of filtering, hindering, and guiding the propagation of acoustic waves (phonons), responsible for the sound propagation when their wavelength is macroscopic and for the thermal transport at room temperature when their wavelength is nanometric.^26,33–36

Many theoretical studies have tried to shed light on the key parameters of the phononic crystals that determine the effect of the interfaces on acoustic and thermal transport. It is worth mentioning that research studies on mechanical and geometrical properties of the nano-interfaces are at the core of the design of the structured materials, involving, for example, the shape and dimensions of the inclusion and the contrast of properties between the materials on the two sides of the interfaces. Recently, we have shown that circular interfaces in a 2D nanophononic material differently affect phonons of different wavelengths, and it is necessary to look at all phonons relevant for heat transport at a certain temperature and their perturbed dynamics, for being able to understand thermal conductivity in such nanocomposites.^29,30 In that work, we could highlight that the rigidity contrast between the two phases is a determinant parameter controlling the strength of the scattering, which affects the phonons with a wavelength comparable with the nanostructuration length scale, significantly anticipating the propagative-to-diffusive crossover in an amorphous matrix. A better understanding of the mechanisms at play and the role of the different parameters could be acquired with a more complete and systematic parametric study on the combined effect of rigidity contrast, interface density, nanostructure length scale, and phonon wavelength. Our results allowed us to identify different transfer regimes (propagative, diffusive, localized, or mixed regimes) in an elastic nanophononic 2D crystal and to show that softer inclusions are more efficient for energy attenuation, but rigid inclusions are also able to pin the vibrational energy at specific frequencies.³⁷ In addition, the existence of an optimal radius of circular inclusion clearly shows that, instead of monotonously increasing the volume fraction until the percolation effect is dominant, other geometrical parameters of the interface are relevant, such as shape, asymmetry, inclusion inter-distance, and interface-to-volume ratio. Finally, when the size of the nanostructure becomes comparable to the wavelength of the excitation, various complex behaviors may occur that affect the acoustic attenuation with non-monotonous frequency dependence. Among such phenomena are the acoustic resonances of the inclusions, multiple reflections between the inclusions, or the interfacial modes.

In contrast to the most work focusing on the periodic circular in 2D (spherical in 3D) interfaces, there are only few investigations conducted with geometrical variations. First, the specific non-aligned arrangement of circular holes (where holes can be considered as the limit of soft inclusion) is reported to have possibly stronger phonon attenuation than the perfectly periodic arrangement. Similarly, the emerging of the gradient-index PCs is meant to independently control the size of each unit hole.³⁸ The impact of disorder on the coherent phonon transport has also been investigated, showing the strong phonon confinement and localization induced by randomness.^16,39 Then, recent experimental studies suggested that the non-circular holes might be more efficient in the thermal conductivity reduction and acoustic attenuation compared to conventional circular ones.^40,41 For example, pacman-shaped holes showed a 40% reduction in thermal conductivity compared to circular ones at room temperature, due to the high surface-to-volume ratio and possible additional resonances.⁴¹ It has also been proposed to optimize the inclusion shape in phononic structures using the homogenized model of strongly heterogeneous elastic composites;⁴² however, this model is limited to simple shapes and topological isomorphism. Topology optimization of metamaterials seems to be a promising alternative^43,44 but is still limited to the current level of nanotechnology manufacturing.

In this work, we propose to investigate the impact of a complex inclusion shape on acoustic propagation and attenuation at nanometric wavelengths. To this purpose, we have chosen to work on a realistic nanocomposite, as can be naturally obtained using vitrification of some metallic glasses. Specifically, we will work on a Ti₄₅Zr₂₅Nb₆Cu₅Be₁₇Sn₂ bulk metallic glass (BMG), which, once produced by casting techniques, exhibits dendrite-phase titanium precipitates at the micrometric scale, as shown in Fig. 1. Bulk metallic glasses are promising structural materials because of their excellent properties such as high yield strength, excellent corrosion resistance, and low stiffness.^45–48 However, at the same time, BMGs lack ductility and always fail in an apparently brittle manner, which seriously limits their applications.^49–52 Impressively, a larger dendrite-phase dimension offers higher ductility but decreases the yield strength of the composites. It has been proposed that the dendrite inclusions actually suppress the catastrophic failure due to the propagation of shear bands and, thus, enhance the global plasticity.^53,54 In addition, the alloy compositional design could be employed to modulate the mechanical properties of the BMG, indicating a potential of tuning, for example, the stiffness ratio between the matrix and the dendrite.⁵⁵ The multi-branch, tree-like geometry of these inclusions is of particular interest for understanding the role of a complex geometry in acoustic attenuation. Compared to the simple shapes, its fractal property determines a remarkably high interface-to-volume ratio, leading to potentially good sound attenuation performance. It is, thus, an ideal sample for our study.

To have a deep understanding of the role of the dendritic shape inclusion in acoustic attenuation and dynamic properties of the phononic crystal, we have performed finite element simulations of out-of-equilibrium acoustic wave-packet propagation in 2D nanocomposite periodically distributed circular and dendritic inclusions. This has allowed us to establish a direct comparison between simple and complex geometries.

This work is organized as follows: in Sec. II, we give a brief introduction on the in situ formed dendrite phase, and we describe how we extract a cluster of dendritic structures from a SEM image; in Sec. III, we describe how the finite element calculations are performed; in Sec. IV, we compare the acoustic attenuation properties in two media containing softer or stiffer inclusions with dendritic and circular shapes; finally, we discuss the results, and we conclude in the last part.

In this section, we will present the preparation of the Ti-based metallic glass containing a dense dendrite phase.^45,55 Based on the SEM images of the samples, instead of looking at the whole dendrite phase, we use a representative cluster of the dendritic structure as the elementary pattern periodically distributed in our model phononic crystal. As for the mechanical properties of the model, we have chosen to use the same properties as in Ref. 37, referring thus to a well-known system. This choice is motivated by the possibility to compare wave-packet propagation results previously obtained on circular inclusions with the ones on dendritic inclusions and, thus, get a direct understanding of the impact of a fractal-like interface shape.

Ingots with the nominal composition Ti₄₅Zr₂₅Nb₆Cu₅Be₁₇Sn₂ are prepared by arc-melting the mixture of high purity elements (>99.9 wt. %) under a Ti-gettered argon atmosphere. The ingots were re-melted at least five times to ensure the homogeneity. Plate samples (5 × 20 × 60 mm³) were prepared by casting into a water-cooled copper mold. The dendrite phase (light gray regions) was found to distribute uniformly within the featureless glass matrix (dark gray regions). Volume fractions of the dendrite phase were analyzed by the Image-Pro Plus software, resulting in a volume fraction of 66% ± 2%. The dendrites have an average size (measured as the diameter of the circular approximation to its shape) of about 30 μm, with a typical single-branch diameter of 3 μm. The size and volume fractions are sensitive to the cooling rate and the alloy composition.⁵⁶ Some mechanical properties are given in  Appendix A (Tables IV and V), and more details on sample preparation and materials properties can be found in Refs. 45 and 55.

In the following, the exact size of the inclusions will not matter much. We will keep only their shape and volume fractions. Thanks to the scalability of our numerical calculations, lengths will be expressed in units of L (the average distance between the inclusions) and frequency in units of ω₀ = 2πc_L/L, with c_L being the longitudinal wave velocity. The inclusions with a reported diameter of 30 µm are not supposed to affect the thermal transport at room temperature, but only acoustic wave propagation at GHz frequencies. However, playing with the exact value of L, it will be possible to also consider the effect of the inclusion shape on thermal transport: for example, by choosing L ≈ nm, the effect on thermal transport at room temperature can be investigated.

In this work, we focus on a cluster of dendritic structures as shown in Fig. 2(1), which is extracted from Fig. 1. It is interesting to focus on this representative zone: its global shape seems to be comparable to a circular inclusion, while its internal tree-like structure may induce different acoustic features. For this aim, a non-dimensional analysis is expected to be carried out by scaling the above cluster of dendrites and the circular inclusion to a quasi-equivalent dimension. As shown in Fig. 2, this cluster of dendrites is encapsulated inside a square block. The square block containing a dendritic inclusion is then used as an elementary brick in the finite element simulation.

Referring to the labels in Fig. 2, the imaging procedure from the SEM image to a finite element mesh is: (1) The region of interest (ROI) is selected and extracted from the original SEM image. (2) Using Matlab toolbox Image Labeler, pixel-level labeling is manually performed in which pixels belonging to either the dendritic inclusion or to the matrix are labeled accordingly. Of course, there are human factors in the steps of determining the two areas, but we have found that the manual method provides a better result than the tested automatic imaging methods such as Canny edge detector,⁵⁷ level sets,⁵⁸ region growing,⁵⁹ and watershed.⁶⁰ (3) The binary labeling information is then transformed into a binary image in which redundant parts are removed. The ROI is the white zone, and the matrix is the black zone. (4) Contour detection gives accurate interfaces between the dendritic cluster and the matrix, and the width of the interfaces is one pixel. (5) The interfaces are segmented into separate zones, and each zone consists of a closed curve. Here, we have 17 independent zones. (6) To form the final geometry, we sequentially import pixel coordinates of the 17 zones into COMSOL Multiphysics to create interpolation curves. Each closed curve creates a part of the dendritic inclusion inside which mechanical properties are homogeneous. The surrounding zone forms the matrix whose mechanical properties are different from the inclusion. Finally, P-1 triangle elements are employed to generate the displayed elementary mesh including both the inclusion and the matrix. It is essential that the number of nodes on the four boundaries of the square are defined a priori as the same and the nodes are equally spaced, for the reason of compatibility, given that we will copy and arrange this elementary mesh in the horizontal direction and implement the periodic boundary condition on the upper and lower boundaries, as shown in Fig. 3.

In this work, we prepared two models: one with dendritic inclusions and another with circular inclusions used for comparison. Since we are interested in the role of inclusion shape, we must ensure the best possible equivalence of the circular and dendritic inclusions, apart from the interface shape. As discussed previously, all the parameters involved in the material constitutive laws (the elastic moduli here) are scale-invariant, resulting in the same invariance for the whole model. We, thus, give all lengths in terms of trivial length L, the distance between the inclusions. Equivalently, the frequency unit is chosen as ω₀ = 2πc_L/L, which is used to define the unit time t₀ = L/c_L in the equation of motion with c_L being the longitudinal wave velocity. The relative size of the inclusions is then chosen in order to reproduce the data already obtained in Ref. 37 for circular inclusions. As in Ref. 37, the apparent diameter of the inclusion is, thus, chosen here as 5/6 L. It corresponds to the exact diameter of the circular inclusions and to the longest axis of the dendritic inclusion. Therefore, the outer contour length is comparable between the two types of inclusions, which can be considered as the primary interface. However, the volume fraction for the dendritic inclusion (shown in Fig. 2) and in 2D corresponding to a surface fraction is then measured as 28% and 35%, while for the circular inclusion, it is 54% and 54%, respectively. The latter, which is close to twice the inclusion area of the former, intuitively allows for a more efficient scattering according to the results in Ref. 37. The geometrical and materials parameters used are summarized in Table I, and materials parameters are shown in Table II.

In the following, we will analyze the role of inclusion shape in acoustic attenuation in periodically arranged nanocomposites. As we are primarily interested in thermal transport, it is possible to scale the dendritic diameter to the nanometric length scale, pertinent for thermal transport at room temperature. For L = 6 nm, for example, as in Ref. 37, the reference frequency will be ω₀ = 8.34 THz, while for L = 30 µm, ω₀ = 1.67 GHz.

We used finite element numerical calculations to study the vibrational properties of a 2D semi-infinite elastic system with dendritic and circular inclusions positioned along a cubic lattice. The computational model consists of nine squares, aligned in the horizontal direction. There is no initial inclusion-free block in this model. The size of each square is defined as L, thus determining the distance between inclusions as L/6. The wave packet is generated imposing a displacement on the left side of the first square around t = 0,

where U₀ is a constant value, ω is the frequency of this quasi-monochromatic excitation, and $t0=3πω$ is the half period of the excitation. A displacement parallel to the boundary corresponds to a transverse excitation, while the one perpendicular to the boundary corresponds to a longitudinal excitation. For the sake of simplicity, we will consider here only longitudinal excitations.

As shown in Fig. 3, periodic boundary conditions (PBCs) are implemented along the vertical direction at the top and bottom of the sample. Perfect Matched Layers (PMLs) are applied on the right side to limit the wave reflection as much as possible. The technical details about the boundary conditions and the time integration scheme can be found in the Appendix in Ref. 37.

Like in our previous works performed on a medium with periodic circular inclusions,⁶¹ the matrix material is linearly elastic with isotropic homogeneous elastic behavior characterized by a typical Young’s modulus of E_m = 92.25 GPa, a mass density of ρ = 2303 kg/m³, and a Poisson ratio of ν = 0.347.⁶² For the inclusions, Poisson’s ratio is supposed to be the same, while Young’s modulus E_i is taken as another control variable and defined as $Ei=Em×EiEm$⁠, the latter being the stiffness ratio of 0.2 or 10. Table II summarizes the values of the parameters used in this work:

A set of transient simulations of longitudinal wave-packet propagation are done using FEM for both media with dendritic and circular inclusions. From the results, we analyze the envelope of the kinetic energy in order to identify the attenuation regime. In addition, the penetration length and diffusivity are calculated to compare the attenuation ability for the two types of inclusions. Finally, long-wavelength and instantaneous sound speeds are estimated for some representative cases. In the following, the kinetic energy E_k is normalized by the maximum value at x = 0.

As said before, the wave packet is created by imposing a displacement on the left side of the sample. Its propagation is then followed along the sample, in the x direction. Due to the presence of interfaces and related spatial inhomogeneities, the wave-packet wave-vector k does not remain constant, the wave packet being scattered by the inclusions. To understand how such scattering affects the energy transfer, we measure the envelope of the kinetic energy induced in the system by the propagation of the wave packet, summed over the y direction. The energy envelope is defined for each excitation frequency ω as

where E_k(x, t) is the instantaneous kinetic energy supported by the frame located in x with a width of Δx = L/60.

In Fig. 4, we present the envelopes of the kinetic energy as a function of position on a semi-log scale and for different normalized frequencies ω/ω₀ with $EiEm=0.2$⁠. Surprisingly, for low frequencies, up to $ωω0=1.2$ included, the circular inclusions strongly attenuate the wave packet, which is almost unaffected by the presence of the dendritic inclusions. The strong effect of circular inclusions has been previously understood as a consequence of the localization of the energy due to a resonance of the inclusions that keeps the acoustic energy.³⁷ However, the irregularly shaped dendritic inclusion does not lead to such resonance phenomena. As such, we can conclude that for low frequencies and long wavelengths, the high interface density of the dendritic inclusion is not dominant in determining the wave-packet attenuation. The situation is reversed above $ωω0$ = 1.56. From this frequency, the dendritic inclusions become more efficient in attenuating the wave packet, and such attenuation strongly increases with increasing frequency, i.e., decreasing wavelength. In the circular case, the attenuation is almost constant for frequencies above $ωω0$ = 2.16, suggesting a saturation of the attenuation effect of circular inclusions in the matrix. We can understand the change in the attenuation efficiency regime as due to the major importance of the tree-like interface in dendritic inclusions at wavelengths comparable with the dendritic structure length scale, i.e., of order of ≈L/10.

We can also identify the different regimes directly from the envelopes of the kinetic energy. Indeed, in the case of weak scattering (propagative regime), a global exponential attenuation similar to the Beer–Lambert law can be observed,^37,63,64

while for strong scattering, for example, due to the larger rigidity contrasts as in this work (E_i/E_m = 0.2 or 10), the algebraic attenuation of the envelope,

is the signature of a diffusive process (diffusive regime).^37,64 Since Fig. 4 is given on a semi-log scale, the propagative regime that follows an exponential decay gives a straight line with a negative slope. In Fig. 4, as exemplified by the green dotted line for ω/ω₀ = 0.60 and 0.96 in the case of dendritic inclusions, an evident exponential decay can be recognized at low frequencies. At higher frequencies, we find 1/x behavior for both the cases of dendritic (ω/ω₀ = 1.20) and circular inclusions (ω/ω₀ = 2.16 and 4.80), which can be then verified on the log–log plot, where the 1/x behavior corresponds to a straight line with a fixed slope, as shown in Fig. 5. Interestingly, for the dendritic case with ω/ω₀ = 1.56, 2.16, and 4.80, the envelope follows neither the B–L fit nor the diffusive fit. If we use a compressed exponential function, such as P ∝  exp(−(x/Λ)^β) with β > 1, the best fit value of β is 1.32, 1.34, and 1.44, respectively. Such a compressed exponential behavior marks a strong difference from the case of circular inclusions, where the diffusive regime is well followed, with a reduction in the oscillation amplitude, as better visible in Fig. 5. In the case of dendritic inclusions, instead, the novel compressed exponential regime exemplifies, in fact, the succession of two regimes: if at the beginning the diffusive law seems to be followed, this leaves space very rapidly to a much stronger attenuation. In both regimes, oscillations are present, similar to the circular inclusion case, with a possibly larger amplitude at high frequency, indicating that the energy is pinned within the inclusions.

In order to get insight into the appearance of this new attenuation regime in the dendritic sample, we report in Fig. 6 the snapshots of the displacement field for low frequency (ω/ω₀ = 0.6), medium frequency (ω/ω₀ = 1.2), and high frequency (ω/ω₀ = 4.8) at the half time (4200 dt) and at the final time (8400 dt) for both samples. For ω/ω₀ = 0.6, the energy is localized and pinned to the first inclusions in the case of circular inclusion but spreads rapidly in the dendritic ones without huge dispersion, meaning that the wavefront propagates ahead followed by an energy tail. At the medium frequency ω/ω₀ = 1.2, the energy is dispersed in space, and the attenuation length in the two samples is quite comparable. At the high frequency ω/ω₀ = 4.8, the energy is scattered violently at the inclusion–matrix interfaces. In the case of circular inclusions, energy is only scattered few times when crossing the circular interface. In addition, due to the large curvature of the circle compared with the short wavelength, there is little but existing energy that keeps spreading ahead along the x direction. On the contrary, for dendritic inclusions, energy is totally scattered due to the random orientations of the normals at the interface and high interface density.

To conclude, different regimes can be identified: (1) exponential attenuation is observed mainly in the low frequency range, which relates directly to the propagative contribution; (2) a diffusive regime is observed at higher frequencies, which can be identified through the P ∝ 1/x behavior and is related directly to the diffusive contribution; (3) a localized (mixed diffusive-localized) regime can be found at some specific frequencies, for example, ω/ω₀ = 0.6 for circular inclusions, which most effectively prevents the energy transport; (4) a compressed exponential attenuation regime is found in the high-frequency range for dendritic inclusions, where the initially diffusive regime seems to be replaced by a stronger attenuation, leading to a global compressed exponential behavior. This results from the combined effect of periodicity and complex interface shape of the dendritic structure.

In Sec. IV A, by comparing the kinetic energy envelopes, we have given the first picture of the attenuation ability in the solids with two types of soft inclusions. Now, we will quantify this attenuation by calculating the penetration length and the diffusivity as a function of frequency in two extreme cases of rigidity contrast, being $EiEm$ = 0.2 and 10.0.

First, we look at the long-time penetration length, defined as the traveled length above which the energy per unit length remains always smaller than the maximum excitation energy per unit length divided by e. An example is reported in Fig. 4 for ω/ω₀ = 0.96, where the abscissa of the intersection between the envelope and the A/e, with A the maximum excitation energy (black solid line), determines the attenuation length.

The results of the normalized penetration length l_p/L are shown in Fig. 7 for both dendritic (red triangles) and circular inclusions (blue circles). In both cases of $EiEm$⁠, at low frequencies, the dendritic inclusions do not exhibit their high-interface-density advantage. Still, as wavelength decreases, the dendrite shape begins showing a better performance to attenuate energy transfer with the reduced crossover frequency $ωcω0$ ranging between 1 and 1.5, depending on the stiffness ratio. Generally, for ω > ω_c, values of the penetration length are systematically shorter in the medium with dendritic inclusions, which indicate a high potential for such samples for reducing the propagative contribution to thermal transport. Interestingly, except for the soft inclusion case at the highest frequency (ω/ω₀ = 4.8), almost all the penetration lengths are greater than the characteristic length L. Finally, it is worth noting that for E_i/E_m = 10, a local minimum appears in the penetration length for both types of inclusions, being around ω/ω₀ = 1 in the circular case and ω/ω₀ = 0.8 in the dendritic case. This effect has already been reported and ascribed to the collective resonance of the nanoparticles in the low frequency range and is related to the effective acoustic impedance of the composite, its frequency position being proportional to $Ei/Em$⁠.³⁷ The slight dependence of this position on the inclusion shape is related to the effective stiffness of the medium. Indeed, the ratio of the two positions is 1.25, which is very close to the ratio between the square roots of the effective stiffnesses (E_eff) of the two media that we calculate in Sec. IV C and can be found in Table III. This suggests the existence of a mixing law of the effective mechanical properties for the bi-phase composites.

In order to get better insight into the strong differences between mean-free paths in the case of circular vs dendritic inclusions, localized at low frequency for soft inclusions and high frequency for hard inclusions, we have calculated using Comsol Multiphysics FEM the phonon band structure, which is reported for the two composites in Figs. 8 and 9 in our frequency window (for ω/ω₀ ≤ 2.5). The periodicity of the nanostructure induces a new Brillouin Zone, with the folding of the intrinsic phonon branches of the matrix material, leading to a raise in a large number of optic modes and the possible appearance of forbidden gaps. Indeed, we can identify some gaps in our investigated frequency window, both in the direction of propagation, corresponding to the Γ − X direction, than in the other high-symmetry directions, X − M and M −Γ. In Figs. 8 and 9, we have reported the band gaps as red stripes and the first five investigated frequencies in green, where each frequency becomes a band due to the spectral broadening intrinsic to our simulation. Indeed, since the phonon wave packet has a finite coherence time (t₀ = 3π/ω), the frequency spectrum is centered at ω with a full width at half maximum of 0.24ω. In the circular case, soft inclusions, gaps are at normalized frequencies smaller than 1, and specifically, our investigated frequency ω/ω₀ = 0.6 falls within a bandgap in the Γ − X direction. Still, this particular mode can propagate in the other directions, indicating that its propagation in the Γ − X direction is, in fact, the projection in this direction of a movement in other directions due to multiple scattering. For E_i/E_m = 10, several gaps appear, in correspondence to our investigated frequencies below 1.2. Interestingly, the gap at ω/ω₀ = 0.6 is now extended to all high-symmetry directions with only slight frequency modulations, indicating that the propagation of this wave packet will be effectively strongly hindered, and indeed, we observe a strong minimum in the mean-free path. Concerning the dendritic case, gaps are generally narrower than those in the circular case and with little coincidence with our investigated frequencies. The lack of the large bandgap at 0.6 for E_i/E_m = 10 in all directions explains the absence of the strong minimum in the mean-free path. At high frequency, ω/ω₀ > 1.2, and no gaps are present in the propagation direction, indicating that if propagation is mainly affected by the modification of the band structure at low frequency, in this range, the dominant effect is the multiple scattering from the fine interface structure of the dendrite. This is, indeed, confirmed by the absence of a wave front in Figs. 6(b) and 6(c).

We have discussed here how the nanostructure affects the wave-packet propagation. As previously mentioned, when it comes to energy transport, there exists a diffusive contribution as well. In order to investigate it in the presence of different inclusions and, thus, understand how thermal transport could be affected, we first identified the situation where the diffusive regime can be recognized from the time evolution of the wave-packet average position ⟨x⟩(t). For this aim, we have calculated this position as

where x_i is the position of the ith frame with width Δx in the x direction and E_k(i, t) is the instantaneous total kinetic energy supported by that frame. In the case of a diffusive process, the squared deviation σ²(t) is proportional to the time t, with a slope related to the one-dimensional diffusivity,

where D is the diffusivity and the spreading σ reads

As an example, σ² of random longitudinal wave packets normalized by L² as a function of time step is shown in Fig. 10 for four frequencies. In the high-frequency range, in both samples, we clearly find a σ² ∝ t relation, indicating the diffusive transport of energy, from which the diffusivity D can be fitted. It is worth underlining that in the sample with dendritic inclusions, we find a diffusive time evolution of the wave packet even at frequencies for which the anomalous stronger-than-diffusive attenuation appears. Concerning the sample with circular inclusions, interestingly, we observe the presence of a short plateau for ω/ω₀ = 0.96, which can be understood as the signature of localization. Indeed, in the localized regime, the wave packet is pinned so that ⟨x⟩ = constant as well as σ², thus giving a plateau. This localized regime is similar to the localization of phonon wave packets observed in previous MD simulations of glasses and disordered phononic crystals (PnCs), for example, in Refs. 39 and 62 

The results of diffusivity are shown in Fig. 11. We can draw similar conclusions as for the penetration length, that is, at high frequencies, the dendritic shape of the inclusions leads to a lower diffusivity than that in the circular case. If this conclusion seems to be universal, independent on the stiffness ratio, it is interesting to note that for $EiEm=10$⁠, the diffusive contribution is reduced by the dendritic shape of the inclusions at all wavelengths. If now we compare the diffusivity between $EiEm=0.2$ and $EiEm=10$⁠, we remark that the overall diffusivity in the soft inclusion case is lower than that of the hard inclusion case, which is coherent with the conclusions in Ref. 37.

In composite nanophononic materials, the wave propagation velocity will be mainly influenced by two factors: the effective rigidity of the medium and the properties of the interfaces. In the following, we will investigate two types of wave speed: (1) effective wave speed at low frequencies and (2) the instantaneous wave speed. The latter is especially useful when it is difficult to get a stationary speed (independent of time), due to the strong elastic heterogeneity of the medium.

Long-wavelength speed for longitudinal waves can be calculated as

where E_eff is effective Young’s modulus of the composite estimated by the Reuss model, which states that the elastic modulus of a composite can be expressed as

where E_i (E_m) is Young’s modulus of the inclusion (matrix) and Φ_i (Φ_m) is the volume fraction of the inclusion (matrix).

We summarize the material properties, effective Young’s modulus, and effective longitudinal wave speeds in Table III. From the estimated wave speeds, it is clear that the long-wavelength wave speed increases compared to the homogeneous solid with E_i/E_m = 1 in the case E_i/E_m > 1 and decreases in the case E_i/E_m < 1 whatever the shape of inclusion is. We find that, both for a more rigid or a softer inclusion, the long-wavelength speed is more strongly affected for a larger volume fraction, i.e., in the case of the circular inclusions. The same conclusion can be drawn from the phonon band diagram, from which we have estimated the phonon density of states (DOS) by counting the number of modes at each frequency interval (Δω/ω₀ = 0.1) over all k-values. It is worth reminding that, following the Debye prediction, the DOS is proportional to the frequency (2D) with a slope inversely proportional to the square of the velocity. As shown in Fig. 12, the slope for E_i/E_m = 0.2 is larger than that of E_i/E_m = 10, indicating a smaller effective sound velocity. It is noted that the difference in slope between two cases for E_i/E_m = 0.2 is also larger than that of E_i/E_m = 10, meaning that the slope of DOS and, thus, the sound velocity for the composites with soft inclusions are more sensitive to the volume fraction. Using the Debye approximation, we can estimate the effective sound velocity for all the investigated elastic contrasts starting from the softest inclusion case and using the scaling law,

where A₀ and c₀ = 4466.1 are the slope and effective sound velocity for the circular inclusions and E_i/E_m = 0.2 and A and c are the slope and effective sound velocity for any other case. For the circular inclusions with E_i/E_m = 10, we, thus, obtain c = 10 842. For the dendritic inclusions, c = 5346 and 8862 for E_i/E_m = 0.2 and 10, respectively, in good agreement with the values in Table III, close to the estimation given by the effective medium approximation. The speed at such wavelengths is, thus, essentially determined by the elastic moduli of the phases and the volume fraction of the secondary phase, independent of the inclusion shape.

This global measurement does not reflect, however, the diversity of behaviors that we have seen in Secs. IV A and IV B and, specifically, the facts that the dendritic shape can strongly affect both propagation and energy diffusion at some wavelengths and the wave packet can be trapped between inclusions. Such differences will affect the wave-packet speed at wavelengths comparable with the nanostructure, i.e., when the wave packet is actually perturbed by the presence of the inclusions in its propagation behavior. For this reason, we address in the following the instantaneous speed.

When the wavelength approaches the nanostructure length scale (size and inter-distance between scatterers), the phonon-interface scattering becomes more important and wave speed begins to deviate from the long-wavelength value. As the wave packet moves in a highly heterogeneous medium, its velocity is not homogeneous in space or in time, the wave packet being scattered in different directions, backward included. Therefore, we need to calculate the instantaneous speed defined as

where 〈x〉 has been defined in Eq. (5) and needs to be smoothed, because the energy oscillates back and forth in the medium, due to the multiple reflections from the interfaces. This can be clearly seen in  Appendix B, Fig. 14. The real, non-smoothed, instantaneous speed will necessarily sharply fluctuate and even assume negative values for back reflections. By smoothing 〈x〉, Eq. (11) will, instead, give us the average speed of energy transport in a relatively short time. In this work, we use the Bezier interpolation as detailed in  Appendix B.

We select four frequencies from low to high: ω/ω₀ = 0.6, 0.96, 1.56, and 4.8, with E_i/E_m = 0.2. The corresponding instantaneous wave speeds are calculated and shown in Fig. 13. Note the initial increase due to the establishment of the wave packet inside the sample, at earlier times for higher frequencies. Only the part after this initial increase must be considered in the following. First, for ω/ω₀ = 0.6 and dendritic inclusions, c_ins exhibits a plateau showing a quasi-constant instantaneous speed of ≈0.45 c_L = 3580 m/s. This speed is always lower than the estimated long-wavelength value reported in Table III, which is 5453 m/s for E_i/E_m = 0.2, indicating a slight attenuation of energy transport. Conversely, this plateau does not exist for the circular inclusion case because it is already in the diffusive-localized regime where energy is pinned in the first inclusions. The instantaneous velocity in the circular inclusion case is always lower than the one in the dendritic inclusion case, meaning that energy moves more slowly in the medium with circular inclusions at all time. This is in agreement with the smaller penetration length found at this frequency in the circular inclusion sample than that in the dendritic one, as shown in Fig. 7(a). For ω/ω₀ = 0.96, the instantaneous speeds in the two samples are quite similar to each other, which, once again, is in good agreement with our findings on the penetration length. ω/ω₀ = 0.96 can be considered as a separation point, because at the following two frequencies ω/ω₀ = 1.56 and 4.8, the instantaneous speed in the circular inclusion samples is systematically larger than that in the dendritic case, meaning that at high frequency, the dendritic sample is more efficient in slowing down and reducing energy transport.

To conclude, the analysis of the instantaneous speed gives a more obvious picture of the crossover of the attenuation performance from the circular inclusion case at low frequencies to the dendritic inclusion case at high frequencies for E_i/E_m = 0.2, especially when the wave-vector no longer exists. Interestingly, this velocity is not stationary in the frequency range studied here.

To conclude, we have shown the effect of a complex shape in the periodic pattern of a 2D nanophononic crystal. We have compared the sound attenuation performance between two shapes of inclusions: circular and dendritic. Our results show that the multi-branching tree-like form of dendrites enhances phonon–interface scattering and phonon attenuation specifically for wavelengths comparable with the dendritic structure length scale regardless of the rigidity ratio. Unlike the circular inclusion, which has only one characteristic length, the sub-interfaces inside the dendritic inclusion provide a continuous source of scattering leading to an increase in sound attenuation. This leads to a stronger reduction in both the penetration length and the apparent diffusivity in samples with dendritic inclusions when the wavelength becomes smaller than the first characteristic length even for far smaller volume fractions of inclusions. Moreover, the instantaneous wave speed is also globally affected, being much reduced at high frequencies by the dendritic fine structure.

It is important to note that the better performance of materials with dendritic vs circular inclusions does not hold at low frequencies, where the larger volume fraction that characterizes the equivalent circular inclusions gives the major contribution to attenuation, besides resonance effects at specific frequencies. If, on the one side, this could suggest to use dendritic inclusions with a larger volume fraction for an optimized nanocomposite, it is worth reminding that the increase in the average inclusion size will translate into an increase in the affected wavelengths and, thus, a decrease in the corresponding frequencies. A compromise between the frequency range that one aims at affecting and the extent of the attenuation, thus, needs to be found.

Concerning the effect of the stiffness ratio between inclusions and the matrix, for soft inclusions, the propagative contribution is reduced also at low frequencies, while for hard inclusions, it is the diffusive that is reduced at all frequencies. Understanding how these modifications can affect thermal transport is not trivial. Indeed, as said, thermal transport has both contributions, but their relative weight depends on the intrinsic properties of the materials and on the temperature at which the study is done.^30,32 Generally, for crystalline materials with little intrinsic phonon attenuation, the propagative contribution dominates at all temperatures, but this is not seen in amorphous or complex crystalline materials. Still, at low temperature, only long-wavelength phonons contribute to thermal transport so that the propagative contribution is dominant. Furthermore, thermal conductivity should be calculated in both contributions (including diffusivity) as a function of temperature, but this is out of the scope of this work.

In addition, as presented in the Introduction, the dendrite phase is formed inside the glassy matrix. In addition, amorphous materials have been used as building blocks for metamaterials for thermal engineering. It was shown recently that an effective elasto–viscous law must be considered in this case, in place of the elastic material behavior used in this work to take into account the intrinsic acoustic attenuation in glasses.^65,66

Finally, it is important to note that the wave dynamics in nanocomposites with a complex geometrical shape such as the dendritic one is much more complex than a simple transition from ballistic to diffusive energy transfer. Indeed, we have shown that above a critical frequency ω_c, the nanocomposite with circular inclusions gives rise to a clear diffusive attenuation combined with a reduction in oscillations, while in the same frequency range, the attenuation appears to be only initially diffusive and then much stronger for the dendritic inclusions while presenting similar oscillations. This is different from the anomalous diffusion of acoustic waves reported in 2D periodic media, characterized by the occurrence of heavy-tailed distribution (as opposed to 1/x decay), interpreted as a consequence of the hybridization of the ballistic and diffusive transport.^32,67 In fact, it can be fitted by a compressed exponential function with β > 1. Therefore, this is truly another type of attenuation process due to the combined effect of periodicity and complex interface shape, giving a stronger attenuation than normal diffusion. The stretched (compressed) exponential character of this attenuation could be related, for example, to the superposition of Beer–Lambert laws with competing mean-free paths due to different attenuation scales.⁶⁸ 

All these results indicate that the use of a complex sub-structure of the interfaces in a phononic material allow us to realize novel optimized materials for acoustic attenuation, leading to applications as high-frequency acoustic filters or thermal insulation, depending on the length scale of the micro(nano)-structure. Indeed, all our findings, reported in normalized units, can be easily scaled at larger or smaller frequencies depending on a smaller or larger nanostructure. As such, our work is more general and gives insights into the universal effect of a complex shape onto acoustic attenuation at all length scales.

See the supplementary material for the videos of wave-packet propagation in Fig. 6.

H.L. was financed by the French Ministry of Research. V.M.G. acknowledges funding from Lyon Idex Breakthrough (IPPON) and the ANR (Grant No. ANR-20-CE05-0046) for research on phonon propagation in nanocomposites. Y.R. acknowledges the NPU Short-term Overseas Visit and Study and the Innovative Experimental Project for postgraduates. A.G. , A.T., and H.L. acknowledge ECOS Chile, Award/Contract Number C17E02.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Mechanical and intrinsic properties of the Ti₄₅Zr₂₅Nb₆Cu₅Be₁₇Sn₂ BMG composites. Young’s modulus (E) and Hardness (H) of the dendrite phase and glass matrix in the Ti₄₅Zr₂₅Nb₆Cu₅Be₁₇Sn₂ BMG composites measured by the nanoindentation (Tables IV and V).

Before calculating the instantaneous speed, we need to pre-treat the data of 〈x〉. Since when waves pass through a deeply heterogeneous medium, 〈x〉 (t) oscillates sharply, and the calculated wave speed can be unreal and meaningless, as illustrated by the yellow dashed line shown in the right panel of Fig. 14. Two smoothing methods are considered here: the first one is nearest neighbor smooth (kernel smoother) defined as

where P_j is the energy envelope at j and n is the number of the nearest neighbors. The second method is Bezier interpolation smoothing. Bezier curves can be defined for any degree n,

where $Cmi$ is equal to the binomial coefficient and m + 1 is equal to the length of the array P. It is reported that the Bezier based smooth curve gives smaller fluctuations and curvatures than other regular smoothing methods in Ref. 69, meaning that it can effectively reduce the oscillations of the first derivative of the 〈x〉.

We used both methods for smoothing 〈x〉 (t) for dendritic inclusions with $EiEm$ = 0.2 and $ωω0$ = 0.6, in which case a wave front can still be identified, as shown in Fig. 6; thus, a well-defined sound speed should be given by a quasi-constant instantaneous speed. The smoothed data of 〈x〉 (t) are shown in the left panel of Fig. 14 by using the two smoothing methods. The Bezier curve is clearly much closer to the real data. In the right panel, the derivative of every 〈x〉 (t) shows that the Bezier interpolation gives the most stable result of wave speed, while the result from the unsmoothed data is useless with such a huge oscillation and the nearest neighbors smoother result is still quite noisy. In addition, a plateau is observed for t ∈ [8, 13] in the case of Bezier curve, which gives a quasi-constant value of the wave speed (around 3500 m/s), confirming the prediction of the existence of a well-defined wave-vector. However, the very beginning and end of the Bezier curve should be ignored, because the Bezier curve must begin and end at given points, i.e., endpoint interpolation property, causing a much sharper oscillation than that with the nearest neighbor smoother at the two ends. Except for those extreme points, the initial stage of acceleration before the plateau corresponds to the establishment step of the wave packet whose duration depends inversely on the wave-packet frequency. Compared to the nearest neighbor smoother, Bezier interpolation gives a clearer presentation and interpretation of the instantaneous speed. Therefore, in the following work, we have chosen the Bezier curve to smooth 〈x〉 to get the instantaneous wave speed.