In building a detailed understanding of the propagation and scattering of waves, that is, their interaction with complex media and boundaries, I have followed a hierarchical description of medium variations or coefficients, starting with a smooth component, a *C*^{1,1} component, general interfaces with regularity constraints, and a class of random fluctuations.

*Multi-scale ray-wave duality*. The description of waves up to*C*^{1,1}coefficients is naturally obtained introducing the notion of wavepackets. Wavepackets were introduced in harmonic analysis and later identified as curvelets. A wavepacket can be viewed as a localized ‘fat’ plane wave, or a ‘quantum’. Wavepackets follow a scale descrip- tion (dyadic*parabolic*decomposition of phase space) and can be used to form a transform pair (generalizing in some sense the notion of beamforming). We obtained a construction of the solution to the wave equation in two steps: in step (i) we smooth the medium on the scales corresponding with the scales of the wavepackets and obtain an approximate construction using geometry (rays on each scale); in step (ii) we incorporate the scattering via a Volterra integral equation. We obtain a concentration of wave packets result; the degree of con- centration reveals the smoothness of the medium [65]. The construction also answers precisely questions about how asymptotic representations are related to full wave solutions.We generalized the construction and description from scalar waves to elastic waves, and addressed the question of the minimum regularity in Lame´ parameters for which polarized waves can be distinguished, which appears to be*C*^{1,1}[86].*Multi-scale wave-equation transmission tomography*. The above construction can be utilized to analyze the physics behind cross-correlation-based wave-equation transmission tomography, that is, the ‘finite-frequency’ analogue of the geodesic ray transform, by taking Fréchet-type derivatives of our solution construction [78], and resolves the seeming insensitivity of the transform along the unperturbed rays via a multi-scale analysis. Indeed, in the case of smooth wavespeeds, and the limit of infinite bandwdith, the geodesic ray transform is recovered [43].*Scattering series, interfaces*. Reducing the regularity of the medium further, leads us to the introduction of (general) interfaces (conormal singularities). We succeeded to track the order of scattering through the regularity of the wave constituents in the case of*C*^{1}(^{,ε}*ε**>*0) coefficients [108]. In this case, the reflected wave is more regular than the incident wave in a Sobolev sense. Via this analysis, we also obtained the first meaningful introduction of a scattering series for waves (beyond the case of planarly layered media).*Random fluctuations, moments, cross correlations*. Via random fluctuations we reduce the regularity of the medium further. Assuming a wavepacket decomposition of the source, we studied waves in random media essentially via a description in terms of moments. Through the parabolic scaling underlying the construction of wavepackets, one obtains Itô-Schrödinger diffusion processes which reproduce the moments, in the limit of fine scales. We incorporated deterministic interfaces as well, moderately tilted relative to the orientation of the wavepacket source, and obtained theoretically an extended notion of enhanced backscattering [101]. We found that the backscattering cone does not depend on the tilt.We utilized the results pertaining to the Itô-Schrödinger diffusion processes and partly coherent waves to study the retrieval of Green’s functions from ‘field-field’ cross correlations in the strongly cluttered regime. We obtain filters via Wigner transforms which satisfy a set of transport equations, and determine their asymptotic behaviors which contain information about the statistics of the random fluctuations. We found that the retrieval is partial only, contains these (blurring) filters, and depends on the directionality of the source [69,102].- We studied the ‘field-field’ correlations from data acquired on the Tibetan plateau when focussing on surface waves. We assessed the retrieval from coda and from ambient noise by comparing the results with direct surface wave observations from earthquakes [54] (including directionality [73]). We carried the comparison over to the estimation of phase velocities and enabled a fusion of the results while enlarging the bandwidth. Here, we employed a semi-classical analysis point of view inserting a small parameter in the variations in medium properties in the surface normal coordinate.
- In
*highly discontinuous and random planarly layered media*, we developed an understanding of tunneling, hidden equivalent medium averaging, and microlocal attenuation which cannot originate from standard visco-elastic models [1,2,5,8].

On planetary scale, that is, a ball, we studied a random radially layered medium including random angular fluctuations on distinct scales as well as deterministic discontinuities. We obtained a stochastic analysis in terms of modes and partly coherent waves excited by an interior source. Through cross correlations of the trace (field at the surface) we obtained a (virtual) data set which can be utilized for the imaging of the deterministic discontinuities.

*Directional wavefield decomposition and propagation, sufficiently smooth coefficients*. In the case of sufficient medium regularity one can interchange time locally with a (curvilinear) spatial coordinate (‘pseudodepth’) to describe the propagation of waves in terms of an evolution or one-way wave equation. In certain non-trivial wavespeed profiles this can be accomplised exactly, using spectral theory (modes) and special functions [23]. Under certain conditions, we can incorporate and analyze the role of ‘evanescent’ wave constituents [24]. In general, this technique is desgned for relatively high frequency wave propagation. In appropriate coordinates, we derived a novel paraxial equation and obtained estimates for its solution in terms of an energy norm converging to the full wave solution as the frequency becomes large. This is accomplished, again, via parabolic scaling. The ‘up-down’ coupling can be incorporated through a generalized Bremmer series [11]. These consistent paraxial equations enable the development of fast algorithms [15], and are of current interest in, for example, long-range propagation of and imaging with teleseismic body waves.*Anisotropy*. We studied medium anisotropy from the point of view of polarized waves and a corresponding local slowness surface characterization (for example, by introducing anellipticity) and expansions [20], and the point of view of classification, characterization and construction of waveforms via Picard-Fuchs equations [7].

The behaviors and phenomena described above provide fundamental insight in developing multi-scale methods to estimate local continuity, anisotropy, regularity etc. of medium properties and substructures from broadband seismic observations, and can lead to a deeper understanding of the coupling of geological and geodynamical processes, (complex) phase boundaries etc. and the interaction (and, hence, probing) with seismic waves. In this context, we obtained some basic results in the case of viscous flow in the mantle wedge induced by a subducting lithospheric plate; we approximated the subduction of the Philippine Sea plate beneath Eurasia, explored the effects of different rheological and thermal constraints, and modelled the evolution of the finite strain ellipses [61].

**Algorithms**

Following the multi-scale physics of waves a family of new computational algorithms has emerged:

- We designed and implemented a
*fast discrete almost symmetric wavepacket transform pair*in three dimensions [80]. This transform pair can also be utilized for wave-based signal processing. Waves (as well as coefficient models containing interfaces) can be compressed with this transform. - Following the dyadic parabolic decomposition of phase space underlying the construction of wave packets, we developed a fast algorithm (with error estimates) for a class of so-called Fourier integral operators comprising propagators, asymptotic solution operators of the wave equations and (inverse) scattering transforms [92]. This algorithm is based on accurate low-rank separated representations,
*prolate spheroidal wave functions*, and ray- geometric computations. (From a physics of waves point of view, in the case of one-way wave propagators, these representations are related to generalizations of phase screens [34,35].) In the application to the half-wave equation, we obtained a theoretical, multi-scale description of the*Gouy phase*which affects the extraction of travel time in long-range propagation. We designed a special procedure to incorporate the (numerical) generation of general*caustics*[104]. - We estimated the regularity of the kernel of the Volterra operator in the above mentioned construction of weak, full-wave solutions in media of limited smoothness. This regularity determined the optimal choice of quadrature to discretize the integral. Moreover we estimated the approximation of solutions by finite sets or clusters of wave packets to complete the discretization of the Volterra equation and obtain a finite matrix representation [95].

We constructed *frames *of wave packets produced by parabolic dilation, rotation and translation of (a finite sum of) Gaussians and obtained asymptotics on the analogue of Daubechies’ frame criterion [110]. We also obtained a frame of wave-atom-like packets generated from pure Gaussians; these can be identified as *multi-scale Gaussian beams*. We developed an associated parametrix construction and algorithm.