Occurrence of anomalous diffusion and nonlocal response in highlyscattering acoustic periodic media
Abstract
We investigate the occurrence of anomalous diffusive transport associated with acoustic wave fields propagating through highlyscattering periodic media. Previous studies had correlated the occurrence of anomalous diffusion to either the random properties of the scattering medium or to the presence of localized disorder. In this study, we show that anomalous diffusive transport can occur also in perfectly periodic media and in the absence of disorder. The analysis of the fundamental physical mechanism leading to this unexpected behavior is performed via a combination of deterministic, stochastic, and fractionalorder models in order to capture the different elements contributing to this phenomenon. Results indicate that this anomalous transport can indeed occur in perfectly periodic media when the dispersion behavior is characterized by anisotropic (partial) bandgaps. In selected frequency ranges, the propagation of acoustic waves not only becomes diffusive but its intensity distribution acquires a distinctive Lévy stable profile having pronounced heavytails. In these ranges, the acoustic transport in the medium occurs according to a hybrid transport mechanism which is simultaneously propagating and anomalously diffusive. We show that such behavior is well captured by a fractional diffusive transport model whose order can be obtained by the analysis of the heavy tails.
I Introduction
In recent years, several theoretical and experimental studies have shown that field transport processes in nonhomogeneous and complex media can occur according to either hybrid or anomalous mechanisms. Some examples of these physical mechanisms include anomalous diffusive transport (such as nonFourier povstenko2013fractional (); borino2011non (); ezzat2010thermoelectric (), or nonFickian diffusion benson2000application (); benson2001fractional (); cushman2000fractional (); fomin2005effect () with heavytailed distribution) or hybrid wave transport (characterized by simultaneous propagation and diffusion mainardi1996fractional (); mainardi1996fundamental (); mainardi1994special (); mainardi2010fractional (); chen2003modified (); chen2004fractional ()). Simultaneous hybrid and anomalous transport has also been observed, particularly in wave propagation problems involving random scattering media. Electromagnetic waves traveling through a scattering materialyamilov2014position () such as fog belin2008display () or murky water zevallos2005time () are relevant examples of practical problems where such transport process can arise.
A distinctive feature of anomalous transport is the occurrence of heavytailed distributions of the representative field quantities benson2001fractional (). In this case, the diffusion process does not follow a classical Gaussian distribution but instead is characterized by a highprobability of occurrence of the events associated with large variance (i.e. those described by the ”heavy” tails).
This behavior is typically not accounted for in traditional field transport theories based on integer order differential or integral models. Purely numerical methods, such as Monte Carlo or finite element simulationshuang1991optical (); ishimaru2012imaging (); mosk2012controlling (); sebbah2012waves (); gibson2005recent (), can capture this response but are very computationally intensive and do not provide any additional insight in the physical mechanisms generating the macroscopic dynamic behavior. The ability to accurately predict the anomalous response and to retrieve information hidden in diffused fields remains a challenging and extremely important topic in many applications. Acoustical and optical imaging, nonintrusive monitoring of engineering and biomedical materials are just a few examples of practical problems in which the ability to carefully predict the field distribution is of paramount importance to achieve accurate and physically meaningful solutions. Nevertheless, in most classical approaches, information contained in the heavy tails is typically discarded because it cannot be properly captured and interpreted by integerorder transport models.
Hybrid and anomalous diffusive transport mechanisms are pervasive also in acoustics. This type of transport can arise when acoustic fields propagate in a highly scattering medium such as a urban environment albert2010effect (); remillieux2012experimental (), a forest aylor1972noise (); tarrero2008sound (), a stratified fluid (e.g. the ocean) baggeroer1993overview (); dowling2015acoustic (); casasanta2012fractional (), or a porous medium benson2001fractional (); schumer2001eulerian (); fellah2003measuring (); fellah2000transient ().
From a general perspective, classical diffusion of wave fields occurs within a range where the wavelength is comparable to the size of the scatterers, the socalled Mie scattering regime. Any deviation from classical diffusion, being either subdiffusion metzler2000random (); goychuk2012fractional () (typically linked to Anderson localization) or superdiffusion (typically linked to Lévyflights) barthelemy2008levy (); bertolotti2010multiple (), still arises within the same regime. The two dominant factors are either the relation between the transport mean free path and the wavelength, or the statistical distributions of the scattering paths in presence of disorder. When a wave field interacts with scattering elements, it undergoes a variety of physical phenomena including reflection, refraction, diffraction, and absorption that significantly alter its initial characteristics. Depending on the quantity, distribution, and properties of the scatterers the momentum vector of an initially coherent wave can become quickly randomized. For most processes, the Central Limit Theorem (CLT) guarantees that the distribution of macroscopic observable quantities (e.g. the field intensity) converges to a Gaussian profile in full agreement with the predictions from classical Fourier diffusion. At the same time, the transition to a macroscopic diffusion behavior leads to an inevitable coexistence of diffusive and wavelike processes at the meso and macroscales.
There are numerous physical processes in nature whose basin of attraction is given by the normal (Gaussian) distribution. On the other hand, when the distribution of characteristic steplength has infinite variance, the diffusion process no longer follows the standard diffusion theory, but rather acquires an anomalous behavior with a basin of attraction given by the socalled stable Lévy distribution. In the latter case, the unbounded value of the variance of the steplength distribution is due to the nonnegligible probability of existence of steps whose lengths greatly differ from the mean value; these are usually denoted as Lévy flights. The distinctive feature of the stable Lévy distributions is the occurrence of heavy tails having a powerlaw decay of the form . This characteristic suggests that transport phenomena evolving according to Lévy statistics are dominated by infrequent but very long steps, and therefore their dynamics are profoundly different from those predicted by the random (Brownian) motion. Many of the complex hybrid transport mechanisms mentioned above fall in this category, and therefore cannot be successfully described in the framework of classical diffusion theory.
In addition, these complex transport mechanisms are typically not amenable to closedform analytical solutions therefore requiring either fully numerical or statistical approaches to predict the field quantities under various input conditions. Typical modeling approaches rely on random walk statistical models metzler2000random (); bouchaud1990anomalous () or on semiempirical corrections to the fundamental diffusive transport equation via renormalization theory asatryan2003diffusion (); cobus2016anderson (). These approaches imply a considerable computational cost and do not provide physical insight in the operating mechanisms of the anomalous response. A few studies have also indicated that, for this type of processes, the macroscopic governing equation describing the evolution of the wave field intensity could be described by a generalization of the classical diffusion equation using fractional derivatives bertolotti2010multiple (); metzler2000random (); bertolotti2007light ().
Todate, the occurrence of anomalous diffusion of wave fields has been connected and observed only in random and disordered media barthelemy2008levy (); burresi2012weak (); bouchaud1990anomalous (); asatryan2003diffusion (); cobus2016anderson (). In this study, we show theoretical and numerical evidence that anomalous behavior can occur even in presence of perfectly periodic media and in absence of disorder. We present this analysis in the context of diffusive transport of acoustic waves although the results could be generalized to other wave fields. In particular, we investigate the specific case of propagation of acoustic waves in a medium with identical and periodically distributed hard scatterers. We develop a theoretical framework for multiple scattering in superdiffusive periodic media. We first show, by full field numerical simulations, that under certain conditions, acoustic waves propagating through a periodic medium are subject to anomalous diffusion. Then, we propose an approach based on a combination of deterministic and stochastic methodologies to explore the physical origin of this unexpected behavior. Ultimately, we show that fractional order models can predict, more accurately and effectively, the resulting anomalous field quantities. More important, we will show that the analysis of the heavy tails provide a reliable means to extract the equivalent fractional order of the medium.
Ii Anomalous diffusion in acoustic periodic media: overview of the method
We consider the generic problem of an acoustic bulk medium made of periodicallydistributed cylindrical hard scatterers in air (Fig. 1). We assume a monopolelike acoustic source, located in the center of the lattice, which emits at a selected harmonic frequency chosen within the scattering regime.
The main objective is to characterize the propagation of acoustic waves in such medium based on different regimes of dispersion. As previously anticipated, in selected regimes the propagation of acoustic waves will exhibit anomalous diffusive transport properties. The reminder of this study will be dedicated to investigating the causes leading to the occurrence of such phenomenon. In order to identify the fundamental mechanisms at the origin of this behavior, we have designed a multifolded approach capable of characterizing the different processes contributing to the macrosopic anomalous response.
The approach consists of the following components. First, we investigate via fullfield numerical simulations the propagation of acoustic waves in either a 1D or a 2D periodic scattering medium. The numerical results will allow making important observations on the different propagation mechanisms occurring in the two systems and on the corresponding diffusive processes. Then, the radiative transfer theory will be applied to interpret the evolution of the wave intensity distribution and analyze the nature of the diffusive phenomena in the context of a renormalization approach.
In order to identify the physical mechanism at the origin of the anomalous diffusion, a multiple scattering analysis based on the multipole expansion method will be applied in order to characterize the interaction between different scatterers. In particular, this approach was intended to identify and quantify possible longrange interactions between pairs of scatterers. Based on the results of the multiple scattering analysis, a Monte Carlo model is used to confirm that the anomalous transport is in fact originated by the longrange interactions between different directions of propagation in the lattice.
Finally, we show that the behavior of the lattice can be effectively described in a homogenized sense, by a fractional continuum diffusion model whose fractional order can be identified by fitting an stable distribution to the heavy tails of the wave intensity. This approach can be seen as an equivalent fractional homogenization of the medium. Of particular interest is the fact that the fractional (homogenized) model allows a closedform analytical solution the agrees very well with the numerical predictions.
Iii Scattering and diffusive transport
From a general perspective, it is possible to identify four different wave propagation regimes in scattering media which are classified based on the relative ratio of quantities such as the transport mean free path , the wavelength of the propagating field , and the characteristic size of the scattering domain . The four regimes are:

The homogenized regime: it occurs when the wavelength of the incident wave field is much larger than the typical characteristic size of the scatterer, that is .

The diffusive regime: it occurs when the wavelength satisfies the relation . In this regime the wave intensity can be approximated by the diffusion equation.

The anomalously diffusive regime: it occurs when the interference of waves causes the reduction of the transport mean free path and consequently a renormalization of the macroscopic diffusion constant . In this regime, the transport mean free path varies according to the size of the cluster and to the degree of disorder.

The localization regime: it occurs in the range and corresponds to a diffusion constant tending to 0.
As mentioned in the classification above, there are regimes in which the intensity of the wave field can be properly described by the diffusion approximation, that is it varies in space as prescribed by the field evolution in a diffusion equation. In particular, when the incident wave has a wavelength smaller than the lengthscale characterizing the material and/or of the geometric variations of the physical medium, the wave field undergoes multiple scattering with a consequent randomization of its phase and direction of propagation. In order to characterize this phenomenon a statistical description based on random walk models is typically employed. These models rely on phenomenological quantities such as the scattering and the transport mean free paths. From a physical perspective, represents the average distance between two successive scattering events, while is the mean distance after which the wave field loses memory of its initial direction and becomes randomized van1999multiple (). When the filling factor (which describes the density of scatterers) is low, and are defined as ishimaru1978wave ():
(1)  
where is the total scattering cross section, is the anisotropy factor and is the scatterers concentration. Note that the relations in Eq. (1) are valid only for low filling factors, approximately in the range . For increased values of the filling factor, the scattering cross section needs to be rescaled. The rescaling factor in the range is given by , while higher filling factors require a more elaborated rescaling procedure ishimaru1978wave ().
The scattering cross section plays a crucial role in the characterization of multiple scattering phenomena and in two dimensions takes the form:
(2) 
The integrand is the differential scattering cross section defined as:
(3) 
.
In Eq. (3), the term is the scattered power flux density at a distance from the scatterer in the direction caused by an incident power flux density . The azimuthal angle is the angle between the incident () and the scattered wave fields ().
The scattering phase function is obtained by normalizing the differential scattering cross section with respect to :
(4) 
and represents the probability that a wave field impinging on the scatterer from a given direction will be scattered by an angle . The mean value of the previous probability distribution defines the anisotropy factor:
(5) 
This factor varies between 0 and 1, and it accounts for the existence of preferential scattering directions. For all the scattering directions have the same probability and the scattering is isotropic. As approaches 1, the forward scattering becomes the most probable event. These quantities will be used in the following analyses in order to identify the different scattering regimes.
Iv Onedimensional medium
Consider a onedimensional bulk scattering medium composed of hard cylindrical scatterers equally distributed in an air background (Fig. 3).
This system can be interpreted by all means as a classical 1D acoustic metamaterial. The radius of the individual scatterer is , where indicates the distance between two neighboring cylinders. The filling fraction for this particular cluster is .
The waveguide is excited by a monochromatic acoustic monopole that replaces the center cylinder. The response of the system is obtained numerically by means of a commercial finite element software (Comsol Multiphysics) and using symmetric boundary conditions on the top and bottom edges and perfectly matched layers (PML) on the left and right edges. The frequency of excitation is selected in the first bandgap (see Fig. 8 for the general dispersion properties of this waveguide) and has a nondimensional value . In this excitation regime the diffusive behavior is expected. Remember that, in the absence of disorder or trapping mechanisms and in the range of excitation frequencies where the diffusion approximation holds, the variance of the steplength distribution characterizing the multiple scattering process of the acoustic field is expected to be finite and, if the steps are independent (by virtue of the Central Limit Theorem) the limit distribution should be the Normal distribution as predicted by the standard diffusion model.
The resulting normalized magnitude of the acoustic pressure field generated in the waveguide is shown in Fig. 4(a) in terms of a contour plot and in Fig. 4(b),(c) in terms of the intensity profile along the midline of the waveguide as defined later in §VI. From Fig. 4(b),(c) a characteristic exponential decay of the type (consistent with the BeerLambert law) is very well identifiable. This trend represents the decay of the coherent part of the intensity and corresponds to the squared absolute value of the Green’s function solution. In more general terms, this result shows a solution which is perfectly consistent with the classical diffusion behavior. This is a well expected result and it is reported here only for comparison with the results that will be presented below.
V Twodimensional medium
v.1 Radial lattice
The immediate extension of the previous scenario to a twodimensional system corresponds to a radial distribution of equally distributed hard scatterers. As in the 1D case, the system is excited by a monochromatic monopole source located in the center of the 2D lattice at point . The source is monochromatic and it is actuated at the nondimensional frequency , that belongs to the first bandgap. The normalized magnitude of the acoustic pressure field for this system is numerically calculated and shown in Fig. 6(a). Fig. 6(b) provides a closeup view of the field around the source (the area within the black dashed line).
The acoustic intensity profile along the axis direction shows an exponential decay as illustrated by Fig. 7. All radial directions (not shown) exhibit an identical response as expected due the azimuthal symmetry of the system.
As in the 1D case, this linear decay of the intensity distribution was expected and confirms that, in systems with a high degree of symmetry, a classical diffusion behavior should be recovered. From a practical perspective, this radial lattice could be seen as a radial arrangement of 1D waveguides.
v.2 Rectangular lattice
The dynamic behavior of the lattice changes quite drastically when the axialsymmetry is removed. Consider the square lattice of scatterers schematically illustrated in Fig. 1. Assume each scatterer having an individual radius of , where is the distance between two neighboring cylinders. The filling fraction for this periodic cluster is .
v.2.1 Dispersion analysis
In order to understand the dynamic behavior of this lattice and interpret the results that will follow, we start analyzing the fundamental dispersion structure of the square lattice. The dispersion was calculated using finite element analysis and the band structure is plotted along the irreducible part of the first Brillouin zone, as shown in Fig. 8.
The results highlight the existence of anisotropy in terms of directions of propagation. These directions are connected to the existence of a partial bandgap in the direction between the nondimensional frequencies and . When the system is excited at a frequency within the bandgap, the propagation acquires an anisotropic distribution (see § V.2.2 and Fig. 9), because propagation can only occur in the direction. This is not an unexpected result and, in fact, it is fully consistent with the propagation behavior expected in square periodic lattice. However, we will show that these dispersion characteristics play a key role in the occurrence of anomalous behavior.
v.2.2 Forced response
The forced response of the lattice was also numerically evaluated. In this case, the lattice is excited by a monochromatic acoustic monopole that replaces the center cylinder. As for the previous two lattices, the total acoustic pressure field is calculated numerically using the finite element method and reported in Fig. 9. More specifically, Fig. 9(a) presents the response to an excitation outside the first bandgap, while Fig. 9(b) reports the case just inside the first bandgap. Note that due to symmetry considerations, only a quarter of the domain was solved.
As the acoustic wave fronts propagate through the medium in the radial directions and interact with the scattering particles, the rays are scattered in multiple directions. In both cases it is evident that the propagation is strongly anisotropic and occurs mostly along the diagonal directions of the lattice.
The response of the medium is shown in Fig. 9 in terms of the normalized magnitude of the acoustic pressure distribution. Contrarily to what observed for the radial lattice, in this case the intensity distribution does not decay linearly. This behavior is very evident by performing a numerical fit of the simulation data, as shown in Fig. 10. These results suggest the occurrence of an unexpected mechanism of diffusion despite the lattice periodicity.
This is a remarkable departure compared with available results in the literature that, todate, have highlighted the occurrence of anomalous diffusion only in connection with random distributions of geometric or material properties.
Vi Radiative transport approach
The results presented above illustrated that in case of anisotropic propagation a departure from the classical diffusive behavior is observed. In this section, we use a traditional radiative transport approach with renormalization to show that this observed behavior can be mapped to anomalous diffusion.
We investigate the presence of anomalous diffusion for wavelength ranges in the passband and in the bandgap. As already pointed out, within the regime the diffusion approximation applies and the spatial evolution of the wave amplitude can be predicted by a diffusion equation for the wave intensity.
Starting from a cluster of particles, as schematically illustrated in Fig. 11, and applying the diffusion approximation the 2D diffusion equation for harmonic excitation and lossless scatterers is given by:
(6) 
where is the intensity of the acoustic wavefield, is the total emitted acoustic power, and are the position vectors of the source and of a generic point , respectively. The average acoustic intensity of a monochromatic monopole source can be obtained as , where is the pressure field, and is the complex conjugate of the velocity field. The diffusion equation Eq. (6) requires the following boundary conditions at the edge of the domain to be solved:
(7) 
where is the unit inward normal. These boundary conditions are obtained by the requirement of zero inward flux at the boundariesishimaru1978wave (). The numerical value of this boundary condition on the intensity was obtained by the previous finite element model.
By enforcing this boundary condition, Eq. (6) can be solved analytically:
(8) 
where is the value of the intensity at the boundary of the cluster of scatterers and is the size of the computational domain.
In order to be able to solve Eq. (8), we need to estimate the parameters and and characterize the specific regime of propagation. To achieve this result, we first plot and versus the wavelength as shown in Fig. 12. These curves were numerically determined using the model presented in §V.2 and the Eqs. (1).
The transport mean free path is always expected to be greater than and to converge asymptotically to for large wavelengths. In fact, for long wavelengths the wavefield is marginally affected by the presence of the scatterers. In the short wavelength limit, tends to 1 because the wave fronts are highly directional (this is the range of validity of ray acoustics approximation) and is approximately given by the average distance between two neighboring scatterers. Fig. 13 shows a detailed view of the previous curves in the frequency range corresponding to the first bandgap and within the diffusive regime. The labels and indicate the nondimensional wavelengths corresponding to the excitation conditions analyzed in the following sections. Note that these curves provide the foundation to investigate the different regimes of propagation and to implement the renormalization approach.
vi.0.1 Renormalization and anomalous diffusion
Fig. 14 shows the acoustic intensity distribution along the axis for the two excitation conditions identified by the labels A and B.
The red circles show the numerical solution obtained by the FE model and provide a onedimensional section of the data in Fig. 9 along the axis. The continuous blue line is the analytical solution of the diffusion equation Eq. (8) after having rescaled the transport mean free path. In particular, for excitation wavelengths in the first passband the value was rescaled to (label in Fig. 13), while for the first bandgap the value was rescaled to (label in Fig. 13).
These results show that, in order to be able to predict the numerical data by using the diffusion approximation, a renormalization of the transport mean free path (and consequently of the diffusion coefficient) must take place. The renormalization requires smaller values of the transport parameters which is a clear indication of superdiffusive anomalous behavior.
Vii Causes of anomalous diffusion
In the previous sections we showed the occurrence of anomalous diffusion of acoustic waves in perfectly periodic square lattices and suggested that the possible origin of this mechanism is linked to the anisotropy of the dispersion properties (i.e. to the anisotropy of the bandgaps).
In this section we will present theoretical and numerical models with the intent of uncovering the physical mechanism leading to this unexpected propagation modality. It is anticipated that the occurrence of anomalous diffusion will be connected to the existence of long range interactions between different directions of propagation governed by either bandpass or stopband behavior. We will use a combination of both deterministic and stochastic methods in order to quantify the longrange interactions and to demostrate that they are at the origin of the macroscopic anomalous diffusion mechanism.
More specifically, we will use a scattering matrix approach to quantify the interaction between different scatterers in different regimes. Then, we will use a discrete random walk diffusion model (which uses probability density functions obtained from the scattering model) to show that, under these assumptions, the anomalous diffusion process matches well with the numerically predicted behavior.
vii.1 The scattering matrix
In order to evaluate and quantify the strength of the interaction between different scatterers in the lattice, we use a multiple scattering approach based on the multipole expansion method. According to this method, after applying the Jacobi’s expansion and the Graf’s addition theorem, the general solution of the wave field can be expressed as:
(9) 
where represents the wave vector, is the angle of the impinging wave field with respect to the axis, is the position vector of the scatterer’s center with respect to the origin of the system of reference, is the position vector of a generic point with respect to the scatterer’s center , , are the Bessel functions of the first kind, are the Hankel functions of the first kind.
To determine the unknown amplitude coefficients the boundary conditions at the surface of the th cylinder must be enforced. The result is a linear set of equations as followslinton2005multiple (); martin2006multiple (); kafesaki1999multiple ():
(10) 
where specifies the Neumann boundary conditions on the surface of the cylinders. The unknown amplitude coefficients can be determined by solving the infinite system of algebraic equations with inner sum truncated at some positive integer . The information about the relative energy exchange between the scatterers can be obtained by rearranging Eq. (10) in matrix form as:
(11) 
where is the unit matrix, is the block diagonal impedance matrix, the vector represents the unknown expansions of scattered waves, and the vector stands for the expansion vector of incident waves on all the scattering cylinders. Finally the matrix is the so called combined translation matrix that can be expressed as follows:
(12) 
where the matrix is defined as follows:
(13) 
The matrix represents the translation matrix between the th and the th cylinder, representing therefore the incident wave on the th cylinder caused by the scattered wave off the th cylinder. The elements of the translation matrix can be obtained from the addition theorem of cylindrical harmonics also known as Graf’s theorem.
The generic term quantifies the portion of the acoustic intensity scattered by the cylinder capable of reaching the cylinder . Equivalently, it represents the fraction of the acoustic intensity reaching the cylinder due to the wave scattered by the cylinder .
This approach was applied to model both the 1D and the 2D waveguides. The normalized scattering coefficients for the 1D waveguide are shown in matrix form in Fig. 16. Each block of this matrix has a size , where is the total number of spherical harmonics used in the multipole series expansion. The total size of the matrix is , where is the total number of cylinders in the cluster. The main diagonal represents the coefficient , that is the scattering of a given cylinder towards itself, and therefore these terms are all zero.
The analysis of Fig. 16 shows, as expected, that in a 1D waveguide the scatterers only interact with their closest neighbors. In other terms, there is no evidence of longrange interaction in 1D periodic waveguides. This is not surprising because we had already found from the fullfield simulations in Fig. 4(a) that the diffusive transport was following a purely Gaussian distribution (hence dominated by nearest neighbor interactions).
In a similar way, the analysis can be repeated for the 2D waveguide with square lattice structure. The resulting scattering coefficients are shown in Fig. 17. Contrarily to the 1D example, this scattering matrix has the appearance of a tridiagonal matrix that highlights the substantial interactions between distant neighbors. In other terms, the rectangular lattice show strong evidence of longrange interactions. These results provide a first important observation concerning the cause of anomalous diffusion in periodic rectangular lattice, that is the anisotropy of the dispersion bands gives rise to longrange interactions that ultimately alter the diffusion process.
vii.2 Discrete random walk models: approximate acoustic intensity
The previous analysis is not yet sufficient to provide conclusive evidence that the longrange interactions due to the bandgap anisotropy are the main cause of the anomalous wave diffusion. In order to identify this further logical link, we developed a discrete random walk (DRW) model capable of simulating the diffusion process resulting from the multiple scattering of the acoustic waves.
The interaction between the different elements of a DRW model is typically represented by probability density functions (pdf). In the following, the pdfs are synthesized based on the coefficients of the scattering matrix. The model can then be numerically solved in order to predict the approximate acoustic intensity resulting from the scattered field.
vii.2.1 1D discrete random walk model
The DRW model for a 1D waveguide is composed of a series of boxes (see Fig. 18), each one representing a scatterer. This model can be seen as the direct discrete equivalent of the 1D waveguide in Fig. 3. The dots in each box represent the different acoustic rays impinging on a given scatterer and being refracted towards different (scattering) elements. This model follows a ray acoustic approximation which is a reasonable assumption in the range of wavelength we have been considering. To simulate the monopole acoustic source located at the center of the waveguide, the center box (labeled ) contains a source term that serves as an omnidirectional source of rays. In the 1D model, the rays emitted from the center box can be scattered both to the left and to the right according to the associated pdf synthesized based on the elements of the scattering matrix.
At every time increment, the rays ”jump” into another box following a Markovian process and a pdf proportional to the coefficients extracted from the scattering matrix. The equilibrium condition needed to solve the DRW model and simulate the evolution of the acoustic intensity upon scattering is given by imposing the conservation of rays:
(14) 
where is the box index, is the time index, is the number of rays at time entering the box (i.e. impinging on the scatterer ), is the source term, and and represent the number of boxes on the left and right side, respectively.
The previous equation can be rearranged as follows:
(15) 
The comparison between the intensity distributions obtained with the FE model and by the equivalent 1D DRW model is given in Fig. 19.
The direct comparison of the results shows a very good agreement between the two models. Note that the DRW is a diffusive model therefore the comparison between the intensity distributions is meaningful only in the tail region. As expected the tails evolve according to a Gaussian distribution. The comparison with the 1D waveguide was provided to illustrate the validity of the proposed approach and to confirm that, under the given assumptions, the results from the DRW converge to the fullfield simulations.
vii.2.2 2D discrete random walk model
The same approach illustrated above for the 1D waveguide can be applied to the analysis of the 2D square lattice. In this case, the DRW model is composed of a 2D distribution of boxes simulating the scatterers. The interactions between different boxes are again expressed in terms of pdfs that are synthesized based on the scattering coefficients obtained from the 2D multipole expansion model (Fig. 17). The equilibrium condition for the 2D DRW model is given by:
(16) 
where and are the box indices, is the time index, is the number of particles at time in the box , is the source term and , , and represent the number of boxes on the left, right, up and down sides, respectively.
The comparison between the intensity distributions obtained by numerical FE simulations and by equivalent 2D DRW model is reported in Fig. 20.
Also in this case, the DRW model is in very good agreement with the FE simulations and, most important, is perfectly capable of capturing the anomalous (powerlaw) decay of the tails of the distribution. This result provides the conclusive proof that the anomalous behavior observed in the square lattice is in fact the result of longrange (Lévy flights) interactions due to scattering events occurring along different directions of propagation that are characterized by anisotropic dispersion.
Viii stable distributions and fractional diffusion equation
The renormalization criterion used in section §VI.0.1 to determine the existence of the anomalous diffusion regime is theoretically wellgrounded but it does not allow a convenient approach to classify the anomalous regime. This classification typically requires the analysis of the time scales involved in the evolution of the moments of the distribution metzler2000random (). Here we suggest a different approach that, not only provides a more direct classification based on the available data, but opens new routes for an analytical treatment of the resulting diffusion problem.
The intensity distributions reported in Fig. 14 suggest a powerlaw behavior of the tails. Recent studies mainardi1996fractional (); mainardi1995fractional (); benson2000application (); benson2001fractional () have shown that, for physical phenomena exhibiting this characteristic distribution of the field variables, the governing equations are generalizations to the fractional order of the classical diffusion equation. Powerlaw distributions, associated with infinite variance random variables (the so called Lévy flights), are in the domain of attraction of stable random variables also called Lévy stable densities (their properties are summarized in Appendix A). On the other hand finitevariance random variables are in the Normal domain of attraction that is a subset of Lévy stable densities. This suggests that the trend of the tails carries information about the stable order of the underlying distribution.
In order to show that this situation occurs also in the present case, we performed numerical fits of the acoustic intensity profiles (Fig. 14) using stable distributions.
The four parameters defining the stable distributions were obtained by numerically solving a nonlinear optimization problem. The most important parameter is the characteristic exponent (also called the index of stability) that is also connected to the slope of the tails. For the square lattice distribution, the values of determined with the optimization procedure are and for the passband and bandgap excitation wavelengths, respectively.
In order to show that the order of the stable distribution effectively describes the anomalous diffusive dynamics of the system, we use a generalized fractional diffusion equation mainardi2007fundamental ():
(17) 
where , , are real parameters always restricted as follows:
(18) 
In Eq. (17), is the field variable, is the RieszFeller space fractional derivative of order , and is the Caputo timefractional derivative of order . The fractional operator in this equation exhibits a nonlocal behavior which makes it ideally suited to model dynamical systems dominated by longrange interactions.
Mainardi mainardi2007fundamental () reported the Green’s function for a Cauchy problem based on the spacetime fractional diffusion equation. The selfsimilar nature of the solution allows the application of a similarity method that separates the solution into a space dependent (the reduced Green’s function ) and a time dependent term. In our system, we use a harmonic (constant amplitude) source and we analyze the steady state response, that is we consider a selfsimilar problem. In other terms, the reduced Green’s function proposed by Mainardi coincides with the normalized solution of the forced fractional diffusion equation governing our problem:
(19) 
Note that this solution is valid in the case . In our case, to model a space fractional diffusion equation.
Fig. 22 shows the comparison between the normalized acoustic intensity from the FE numerical data and the result from the reduced Green’s function (Eq. (19)) calculated for the order obtained by the previous stable fits.
The above results clearly show that the fractional diffusion equation is able to capture the heavytailed behavior of the intensity distribution with good accuracy. They also confirm that the use of stable fits provides a reliable approach to classify the anomalous behavior and to extract the corresponding fractional order of the operator. The above numerical results provide also further confirmation that the observed dynamic behavior from the full field numerical simulations is in fact dominated by anomalous diffusion. These results are particularly relevant if seen in a perspective of developing predictive capabilities for transport processes in highly inhomogeneous systems. As an example, fractional models would provide an excellent framework for the solution of inverse problems in imaging and remote sensing through highly scattering media. The ability to properly capture a mixed transport behavior, such as partially propagating and diffusive, would allow extracting more information from the measured response therefore improving the sensitivity and resolution of these approaches. From a broader perspective, this methodology has general applicability and could be extended to a variety of applications involving wavelike field transport such as those mentioned in the introduction.
Ix Conclusions
In this paper, we investigated the scattering behavior of sound waves in a perfectly periodic acoustic medium composed of a square lattice of hard cylinders in air. From a general perspective, the most remarkable result consists in the observation of the occurrence of anomalous hybrid transport in perfectly periodic lattice structures without disorder or random properties. This result is particularly relevant because the anomalous response of a scattering system was previously observed only in systems with either stochastic material or geometric properties. By using a combination of theoretical and numerical models, both deterministic and stochastic, it was determined that the existence of longrange interactions associated with the anisotropy of the dispersion bands was the driving factor leading to the occurrence of the anomalous transport behavior. The resulting diffused intensity fields were characterized by heavytails with marked asymptotic powerlaw decay, that were well described by stable distributions. It was also shown that the stable nature of the dynamic response provided a reliable approach for the classification and characterization of the nonlocal effect via the intrinsic parameters of stable distributions.
Observing that stable distributions represent the fundamental kernel for the solutions of fractional continuum models, we showed that a space fractional diffusion equation having the order predicted by the stable fit of the acoustic intensity was capable of capturing very accurately the characteristic features of the anomalous transport process. From a general perspective, this approach can be interpreted as a fractional order homogenization of the periodic medium which is capable of mapping the complex inhomogeneous system to a (fractional) governing equation that still accepts an analytical solution.
This latter characteristic is particularly remarkable if seen from a practical application perspective because it could open the way to accurate and noniterative inverse problems that play a critical role in remote sensing, imaging, and material design, just to name a few. Another key observation concerns the strong deviation of the tails of the acoustic intensity from the Gaussian distribution which highlights that much information is still contained in the tails. This aspect is particularly relevant for imaging and sensing in scattering media because traditional analytical methodologies typically assume a Gaussian distribution of the measured intensity field hence leading to two main drawbacks: 1) the loss of important information about the internal structure of the medium which is contained in the tails, and 2) the lack of a proper model capable of extracting and interpreting this information from measured data.
Acknowledgments
The authors gratefully acknowledge the financial support of the National Science Foundation under the grants DCSD CAREER and of the Air Force Office of Scientific Research under Grant No. YIP FA95501510133.
*
Appendix A Appendix A
This appendix summarizes some basic properties of the stable distributions that have been used to analyze and interpret the simulation data in the paper. The family of stable distributions are defined by the Fourier transform of their characteristics functions that can be written in the explicit form asherranz2004alpha (); benson2002fractional ():
(20)  
where , , , and . The parameters , , and uniquely and completely identify the stable distribution.

The parameter is the characteristic exponent, or the stability parameter, and it defines the degree of impulsiveness of the distribution. As decreases the level of impulsiveness of the distribution increases. For we recover the Gaussian distribution. A particular case is obtained for and that corresponds to the Cauchy distribution. For the inverse Fourier transform is not positivedefinite and hence is not a proper probability density function.

The parameter is the symmetry, or skewness parameter, and determines the skewness of the distribution. Symmetric distributions have , whereas and correspond to completely skewed distributions.

The parameter is the scale parameter. It is a measure of the spread of the samples from a distribution around the mean.

The parameter is the location parameter and corresponds to a shift in the axis of the pdf. For a symmetric distribution, is the mean when and the median when .
The characteristic functions described in Eq. (20) are equivalent to a probability density function and do not have analytical solutions except for few special cases. The main feature of these characteristic functions is the presence of heavytails when compared to a Gaussian distribution. The probability density functions with tails heavier than Gaussian are also denoted as impulsive. An impulsive process is characterized by the presence of large values that significantly deviates from the mean value of the distribution with nonnegligible probability. In this sense the stable distribution represents a generalization of the Gaussian distribution that allows to model impulsive processes by using only four parameters instead of an infinite number of moments. The possibility of describing the distribution of particles in anomalous diffusion phenomena by using stable distributions has numerous advantages: 1) many methods exist to perform statistical inference on stable environments nikias1995signal (); janicki1993simulation (), 2) these distributions are simple because they are completely characterized by only four parameters, 3) the use of stable distributions finds a theoretical justification in the fact that they satisfy the generalized central limit theorem which states that the limit distribution on infinitely many i.i.d. random variables, is a stable distribution, 4) they include the Gaussian distribution as a particular case for a specific set of parameters. These distributions are stable since the output of a linear system in response to stable inputs is again stable.
References
 [1] Y. Z. Povstenko. Fractional heat conduction in infinite onedimensional composite medium. Journal of Thermal Stresses, 36(4):351–363, 2013.
 [2] G. Borino, M. Di Paola, and M. Zingales. A nonlocal model of fractional heat conduction in rigid bodies. The European Physical JournalSpecial Topics, 193(1):173–184, 2011.
 [3] M. A. Ezzat. Thermoelectric mhd nonnewtonian fluid with fractional derivative heat transfer. Physica B: Condensed Matter, 405(19):4188–4194, 2010.
 [4] D. A. Benson, S. W. Wheatcraft, and M. M. Meerschaert. Application of a fractional advectiondispersion equation. Water Resources Research, 36(6):1403–1412, 2000.
 [5] D. A. Benson, R. Schumer, M. M. Meerschaert, and S. W. Wheatcraft. Fractional dispersion, lévy motion, and the made tracer tests. In Dispersion in Heterogeneous Geological Formations, pages 211–240. Springer, 2001.
 [6] J. H. Cushman and T. R. Ginn. Fractional advectiondispersion equation: a classical mass balance with convolutionfickian flux. Water resources research, 36(12):3763–3766, 2000.
 [7] S. Fomin, V. Chugunov, and T. Hashida. The effect of nonfickian diffusion into surrounding rocks on contaminant transport in a fractured porous aquifer. In Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, volume 461, pages 2923–2939. The Royal Society, 2005.
 [8] F. Mainardi. Fractional relaxationoscillation and fractional diffusionwave phenomena. Chaos, Solitons & Fractals, 7(9):1461–1477, 1996.
 [9] F. Mainardi. The fundamental solutions for the fractional diffusionwave equation. Applied Mathematics Letters, 9(6):23–28, 1996.
 [10] F. Mainardi and M. Tomirotti. On a special function arising in the time fractional diffusionwave equation. Transform methods and special functions, Sofia, 171, 1994.
 [11] F. Mainardi. Fractional calculus and waves in linear viscoelasticity: an introduction to mathematical models. World Scientific, 2010.
 [12] W. Chen and S. Holm. Modified szaboâs wave equation models for lossy media obeying frequency power law. The Journal of the Acoustical Society of America, 114(5):2570–2574, 2003.
 [13] W. Chen and S. Holm. Fractional laplacian timespace models for linear and nonlinear lossy media exhibiting arbitrary frequency powerlaw dependency. The Journal of the Acoustical Society of America, 115(4):1424–1430, 2004.
 [14] A. Yamilov, R. Sarma, B. Redding, B. Payne, H. Noh, and H. Cao. Positiondependent diffusion of light in disordered waveguides. Physical review letters, 112(2):023904, 2014.
 [15] E. Belin, F. Christnecher, F. Taillade, and M. Laurenzis. Display of an analytical model for backscattered luminance and a fullfield range gated imaging system for vision in fog. In Optical Engineering+ Applications, pages 70880O–70880O. International Society for Optics and Photonics, 2008.
 [16] Manuel E. Zevallos L., S. K. Gayen, M. Alrubaiee, and R. R. Alfano. Timegated backscattered ballistic light imaging of objects in turbid water. Applied Physics Letters, 86(1):011115, 2005.
 [17] D. Huang, E. A. Swanson, C. P. Lin, J. S. Schuman, W. G. Stinson, W. Chang, M. R. Hee, T. Flotte, K. Gregory, C. A. Puliafito, et al. Optical coherence tomography. Science (New York, NY), 254(5035):1178, 1991.
 [18] A. Ishimaru, S. Jaruwatanadilok, and Y. Kuga. Imaging through random multiple scattering media using integration of propagation and array signal processing. Waves in Random and Complex Media, 22(1):24–39, 2012.
 [19] A. P. Mosk, A.D. Lagendijk, G. Lerosey, and M. Fink. Controlling waves in space and time for imaging and focusing in complex media. Nature photonics, 6(5):283–292, 2012.
 [20] P. Sebbah. Waves and imaging through complex media. Springer Science & Business Media, 2012.
 [21] A. P. Gibson, J. C. Hebden, and S. R. Arridge. Recent advances in diffuse optical imaging. Physics in medicine and biology, 50(4):R1, 2005.
 [22] D. G. Albert and L. Liu. The effect of buildings on acoustic pulse propagation in an urban environment. The Journal of the Acoustical Society of America, 127(3):1335–1346, 2010.
 [23] M. C. Remillieux, J. M. Corcoran, T. R. Haac, R. A. Burdisso, and U. P. Svensson. Experimental and numerical study on the propagation of impulsive sound around buildings. Applied Acoustics, 73(10):1029–1044, 2012.
 [24] D. Aylor. Noise reduction by vegetation and ground. The Journal of the Acoustical Society of America, 51(1B):197–205, 1972.
 [25] A. I. Tarrero, M. A. Martín, J. González, M. Machimbarrena, and F. Jacobsen. Sound propagation in forests: A comparison of experimental results and values predicted by the nord 2000 model. Applied Acoustics, 69(7):662–671, 2008.
 [26] A. B. Baggeroer, W. A. Kuperman, and P. N. Mikhalevsky. An overview of matched field methods in ocean acoustics. IEEE Journal of Oceanic Engineering, 18(4):401–424, 1993.
 [27] D. R. Dowling and K. G. Sabra. Acoustic remote sensing. Annual Review of Fluid Mechanics, 47:221–243, 2015.
 [28] G. Casasanta and R. Garra. Fractional calculus approach to the acoustic wave propagation with spacedependent sound speed. Signal, Image and Video Processing, 6(3):389–392, 2012.
 [29] R. Schumer, D.A. Benson, M.M. Meerschaert, and S.W. Wheatcraft. Eulerian derivation of the fractional advection–dispersion equation. Journal of Contaminant Hydrology, 48(1):69–88, 2001.
 [30] Z. E. A. Fellah, S. Berger, W. Lauriks, C. Depollier, and M. Fellah. Measuring the porosity of porous materials having a rigid frame via reflected waves: A time domain analysis with fractional derivatives. Journal of Applied Physics, 93(1):296–303, 2003.
 [31] Z. E. A. Fellah and C. Depollier. Transient acoustic wave propagation in rigid porous media: a timedomain approach. The Journal of the Acoustical Society of America, 107(2):683–688, 2000.
 [32] R. Metzler and J. Klafter. The random walk’s guide to anomalous diffusion: a fractional dynamics approach. Physics reports, 339(1):1–77, 2000.
 [33] I. Goychuk. Fractionaltime random walk subdiffusion and anomalous transport with finite mean residence times: Faster, not slower. Physical Review E, 86(2):021113, 2012.
 [34] P. Barthelemy, J. Bertolotti, and D. S. Wiersma. A Lévy flight for light. Nature, 453(7194):495–498, 2008.
 [35] J. Bertolotti, K. Vynck, and D. S. Wiersma. Multiple scattering of light in superdiffusive media. Physical review letters, 105(16):163902, 2010.
 [36] J.P. Bouchaud and A. Georges. Anomalous diffusion in disordered media: statistical mechanisms, models and physical applications. Physics reports, 195(45):127–293, 1990.
 [37] A. A. Asatryan, P. A. Robinson, R. C. McPhedran, L. C. Botten, C. M. de Sterke, T. L. Langtry, and N. A. Nicorovici. Diffusion and anomalous diffusion of light in twodimensional photonic crystals. Physical Review E, 67(3):036605, 2003.
 [38] L. A. Cobus. Anderson localization and anomalous transport of ultrasound in disordered media. 2016.
 [39] J. Bertolotti. Light transport beyond diffusion. These de Doctorat, Universita degli Studi di Firenze, 2007.
 [40] M. Burresi, V. Radhalakshmi, R. Savo, J. Bertolotti, K. Vynck, and D. S. Wiersma. Weak localization of light in superdiffusive random systems. Physical review letters, 108(11):110604, 2012.
 [41] M. C. W. van Van Rossum and T. M. Nieuwenhuizen. Multiple scattering of classical waves: microscopy, mesoscopy, and diffusion. Reviews of Modern Physics, 71(1):313, 1999.
 [42] A. Ishimaru. Wave propagation and scattering in random media, volume 2. Academic press New York, 1978.
 [43] C. M. Linton and P. A. Martin. Multiple scattering by random configurations of circular cylinders: Secondorder corrections for the effective wavenumber. The Journal of the Acoustical Society of America, 117(6):3413–3423, 2005.
 [44] P. A. Martin. Multiple scattering: interaction of timeharmonic waves with N obstacles. Number 107. Cambridge University Press, 2006.
 [45] M. Kafesaki and E. N. Economou. Multiplescattering theory for threedimensional periodic acoustic composites. Physical review B, 60(17):11993, 1999.
 [46] F. Mainardi. Fractional diffusive waves in viscoelastic solids. Nonlinear waves in solids, 137:93–97, 1995.
 [47] F. Mainardi, Y. Luchko, and G. Pagnini. The fundamental solution of the spacetime fractional diffusion equation. arXiv preprint condmat/0702419, 2007.
 [48] D. Herranz, E. E. Kuruoğlu, and L. Toffolatti. An stable approach to the study of the P(D) distribution of unresolved point sources in CMB sky maps. Astronomy & Astrophysics, 424(3):1081–1096, 2004.
 [49] D.A. Benson, R. Schumer, M.M. Meerschaert, and S.W. Wheatcraft. Fractional dispersion, Lévy motion, and the MADE tracer tests. In Dispersion in Heterogeneous Geological Formations, pages 211–240. Springer, 2002.
 [50] C. L. Nikias and M. Shao. Signal processing with alphastable distributions and applications. WileyInterscience, 1995.
 [51] A. Janicki and A. Weron. Simulation and chaotic behavior of alphastable stochastic processes, volume 178. CRC Press, 1993.