Presence Of Caustics Research Articles

Anisotropy affects the lattice waves and phonon distributions in GaAs

The stress-strain relationship, anisotropy, heavy- and light-hole structures, and the Phillips ionicity in GaAs are theoretically investigated by Bond matrices model as well as Zener and Every anisotropies. The results show that there is a simple correspondence between anisotropy and the extreme directions of these physical properties, i.e., [0 0 1], [1 1 0] and [1 1 1] directions. Meanwhile the lattice wave propagation, phonon focusing, and phonon distributions in GaAs are also theoretically studied in detail based on lattice-dynamical method. The slowness surface of three lattice waves are the mixing of longitudinal and transverse modes. The feature is that the lattice waves, saddle Gaussian curvatures, and caustic structures are modulated by the anisotropy. The topology in phonon distribution is discussed by the presence of caustics in the anisotropic flux of phonons emanating from a localized phonon distribution. Their topological structures show the phonon distribution and the phonon focusing. The anisotropy not only affects lattice wave propagation and phonon distribution, but influences the phonon focusing.

Elastic properties and lattice waves of TiAl

The elastic properties, lattice wave propagation and phonon focusing of TiAl are studied in detail based on Bond matrices model and lattice-dynamical method. The results contain remarkably complex three-dimensional variations in the phonon flux with angular distribution. The feature of elastic waves and phonon distribution is affected by the elastic anisotropy of the lattice. The correspondence between these features and contours of Gaussian curvature on the slowness surface is explored and the phonon topology is discussed about the presence of caustics in the anistropic flux of phonons emanating from a localized phonon distribution. The topology slowness surface gives insights into the mixing of longitudinal and transverse modes and gives the origin of the phonon focusing. The elastic anisotropy not only affects lattice wave propagation but also influences phonon focusing.

Error estimates for Gaussian beam methods applied to symmetric strictly hyperbolic systems

In this work we construct Gaussian beam approximations to solutions of linear symmetric hyperbolic systems with highly oscillatory initial data, including both strictly and non-strictly hyperbolic systems as long as they are diagonalizable. The evolution equations for each Gaussian beam component are derived. Under some regularity assumptions of the data we obtain an error estimate between the exact solution and the first order Gaussian beam superposition in terms of the high frequency parameter ε−1. The main result is that the relative local error measured in energy norm in the beam approximation decays as ε12 independent of dimension and presence of caustics, for first order beams. This result is shown to be valid when the gradient of the initial phase is bounded away from zero. Applications to some non-strictly hyperbolic systems including both the acoustic equation and the system of Maxwell equations are discussed.

A numerical study of nonlinear infrasound propagation in a windy atmosphere.

Direct numerical simulations of the two-dimensional unsteady compressible Navier-Stokes equations are performed to study the acoustic field generated by an infrasonic source in a realistic atmosphere. Some of the main phenomena affecting the propagation of infrasonic waves at large distances from the source are investigated. The effects of thermal and wind-related refraction on the signals recorded at ground level are highlighted, with particular emphasis on the phase shift induced by the presence of caustics in the acoustic field. Nonlinear waveform steepening associated with harmonic generation, and period lengthening, both of which are typical of large source amplitudes, are illustrated, and the importance of thermoviscous absorption in the upper atmosphere is clearly demonstrated. The role of diffraction in the shadow zone, around caustics and at stratospheric altitudes is also pointed out. The Navier-Stokes equations are solved using high-order finite-differences and a Runge-Kutta time integration method both originally developed for aeroacoustic applications, along with an adaptive shock-capturing algorithm which allows high-intensity acoustic fields to be examined. An improvement to the shock detection procedure is also proposed in order to meet the specificities of nonlinear propagation at long range. The modeling as well as the numerical results are reported in detail and discussed.

Babich's expansion and the fast Huygens sweeping method for the Helmholtz wave equation at high frequencies

In some applications, it is reasonable to assume that geodesics (rays) have a consistent orientation so that the Helmholtz equation can be viewed as an evolution equation in one of the spatial directions. With such applications in mind, starting from Babich's expansion, we develop a new high-order asymptotic method, which we dub the fast Huygens sweeping method, for solving point-source Helmholtz equations in inhomogeneous media in the high-frequency regime and in the presence of caustics. The first novelty of this method is that we develop a new Eulerian approach to compute the asymptotics, i.e. the traveltime function and amplitude coefficients that arise in Babich's expansion, yielding a locally valid solution, which is accurate close enough to the source. The second novelty is that we utilize the Huygens–Kirchhoff integral to integrate many locally valid wavefields to construct globally valid wavefields. This automatically treats caustics and yields uniformly accurate solutions both near the source and remote from it. The third novelty is that the butterfly algorithm is adapted to accelerate the Huygens–Kirchhoff summation, achieving nearly optimal complexity O(Nlog⁡N), where N is the number of mesh points; the complexity prefactor depends on the desired accuracy and is independent of the frequency. To reduce the storage of the resulting tables of asymptotics in Babich's expansion, we use the multivariable Chebyshev series expansion to compress each table by encoding the information into a small number of coefficients. The new method enjoys the following desired features. First, it precomputes the asymptotics in Babich's expansion, such as traveltime and amplitudes. Second, it takes care of caustics automatically. Third, it can compute the point-source Helmholtz solution for many different sources at many frequencies simultaneously. Fourth, for a specified number of points per wavelength, it can construct the wavefield in nearly optimal complexity in terms of the total number of mesh points, where the prefactor of the complexity only depends on the specified accuracy and is independent of frequency. Both two-dimensional and three-dimensional numerical experiments have been carried out to illustrate the performance, efficiency, and accuracy of the method.

Eulerian Geometrical Optics and Fast Huygens Sweeping Methods for Three-Dimensional Time-Harmonic High-Frequency Maxwell's Equations in Inhomogeneous Media

In some applications, it is reasonable to assume that geodesics (rays) have a consistent orientation so that Maxwell's equations may be viewed as an evolution equation in one of the spatial directions. With such applications in mind, we propose a new Eulerian geometrical-optics method, dubbed the fast Huygens sweeping method, for computing Green's functions of Maxwell's equations in inhomogeneous media in the high-frequency regime and in the presence of caustics. The first novelty of the fast Huygens sweeping method is that a new dyadic-tensor-type geometrical-optics ansatz is proposed for Green's functions which is able to utilize some unique features of Maxwell's equations. The second novelty is that the Huygens--Kirchhoff secondary source principle is used to integrate many locally valid asymptotic solutions to yield a globally valid asymptotic solution so that caustics associated with the usual geometrical-optics ansatz can be treated automatically. The third novelty is that a butterfly algorithm is adapted to carry out the matrix-vector products induced by the Huygens--Kirchhoff integration in $O(N\log N)$ operations, where $N$ is the total number of mesh points, and the proportionality constant depends on the desired accuracy and is independent of the frequency parameter. To reduce the storage of the resulting traveltime and amplitude tables, we compress each table into a linear combination of tensor-product based multivariate Chebyshev polynomials so that the information of each table is encoded into a small number of Chebyshev coefficients. The new method enjoys the following desired features: (1) it precomputes a set of local traveltime and amplitude tables; (2) it automatically takes care of caustics; (3) it constructs Green's functions of Maxwell's equations for arbitrary frequencies and for many point sources; (4) for a specified number of points per wavelength it constructs each Green's function in nearly optimal complexity $O(N\log N)$ in terms of the total number of mesh points $N$, where the prefactor of the complexity depends only on the specified accuracy and is independent of the frequency parameter. Three-dimensional numerical experiments are presented to demonstrate the performance and accuracy of the new method.

Sobolev and max norm error estimates for Gaussian beam superpositions

This work is concerned with the accuracy of Gaussian beam superpositions, which are asymptotically valid high frequency solutions to linear hyperbolic partial differential equations and the Schr\"odinger equation. We derive Sobolev and max norms estimates for the difference between an exact solution and the corresponding Gaussian beam approximation, in terms of the short wavelength $\varepsilon$. The estimates are performed for the scalar wave equation and the Schr\"odinger equation. Our result demonstrates that a Gaussian beam superposition with $k$-th order beams converges to the exact solution as $O(\varepsilon^{k/2-s})$ in order $s$ Sobolev norms. This result is valid in any number of spatial dimensions and it is unaffected by the presence of caustics in the solution. In max norm, we show that away from caustics the convergence rate is $O(\varepsilon^{\lceil k/2\rceil})$ and away from the essential support of the solution, the convergence is spectral in $\varepsilon$. However, in the neighborhood of a caustic point we are only able to show the slower, and dimensional dependent, rate $O(\varepsilon^{(k-n)/2})$ in $n$ spatial dimensions.

Refractometry of optically inhomogeneous media based on recording the position of caustics using structured laser radiation

This paper deals with the refractometry of optically inhomogeneous media under conditions of significant refraction of probing beams due to substantial gradients of the refractive index and the length of the test medium. For these conditions, wave and beam refraction models of laser beams in the presence of caustics are developed. The advantages of using structured beams to record caustics for the purpose of quantitative diagnostics of inhomogeneities of the refractive index are discussed. Based on the developed models, possible methods are proposed to solve the inverse refraction problem and reconstruct the refractive index values in inhomogeneity cross-sections.

Fast Huygens sweeping methods for Helmholtz equations in inhomogeneous media in the high frequency regime

In some applications, it is reasonable to assume that geodesics (rays) have a consistent orientation so that the Helmholtz equation may be viewed as an evolution equation in one of the spatial directions. With such applications in mind, we propose a new Eulerian computational geometrical-optics method, dubbed the fast Huygens sweeping method, for computing Green functions of Helmholtz equations in inhomogeneous media in the high-frequency regime and in the presence of caustics. The first novelty of the new method is that the Huygens–Kirchhoff secondary source principle is used to integrate many locally valid asymptotic solutions to yield a globally valid asymptotic solution so that caustics associated with the usual geometrical-optics ansatz can be treated automatically. The second novelty is that a butterfly algorithm is adapted to carry out the matrix–vector products induced by the Huygens–Kirchhoff integration in O(NlogN) operations, where N is the total number of mesh points, and the proportionality constant depends on the desired accuracy and is independent of the frequency parameter. To reduce the storage of the resulting traveltime and amplitude tables, we compress each table into a linear combination of tensor-product based multivariate Chebyshev polynomials so that the information of each table is encoded into a small number of Chebyshev coefficients.The new method enjoys the following desired features: (1) it precomputes a set of local traveltime and amplitude tables; (2) it automatically takes care of caustics; (3) it constructs Green functions of the Helmholtz equation for arbitrary frequencies and for many point sources; (4) for a specified number of points per wavelength it constructs each Green function in nearly optimal complexity in terms of the total number of mesh points, where the prefactor of the complexity only depends on the specified accuracy and is independent of the frequency parameter.Both two-dimensional (2-D) and three-dimensional (3-D) numerical experiments are presented to demonstrate the performance and accuracy of the new method.

Optimal Error Estimates for First-Order Gaussian Beam Approximations to the Schrödinger Equation

Gaussian beams are generally local asymptotic solutions to the linear wave equations in the high-frequency regime. Each Gaussian beam is concentrated around a specific ray path determined by the underlying Hamiltonian system. Expressed as some superposition of Gaussian beams, Gaussian beam approximation is expected to be a high-frequency asymptotic solution which remains globally valid even around caustics. We derive optimal first-order error estimates for first-order Gaussian beam approximations, in both continuous and discrete levels, to the Schrodinger equation equipped with a WKB initial data. Our error estimates are valid for any spatial dimension and unaffected by the presence of caustics. Some numerical tests are presented to validate the theoretical results.

The Principle of Stationary Action in Biophysics: Stability in Protein Folding

We conceptualize protein folding as motion in a large dimensional dihedral angle space. We use Lagrangian mechanics and introduce an unspecified Lagrangian to study the motion. The fact that we have reliable folding leads us to conjecture the totality of paths forms caustics that can be recognized by the vanishing of the second variation of the action. There are two types of folding processes: stable against modest perturbations and unstable. We also conjecture that natural selection has picked out stable folds. More importantly, the presence of caustics leads naturally to the application of ideas from catastrophe theory and allows us to consider the question of stability for the folding process from that perspective. Powerful stability theorems from mathematics are then applicable to impose more order on the totality of motions. This leads to an immediate explanation for both the insensitivity of folding to solution perturbations and the fact that folding occurs using very little free energy. The theory of folding, based on the above conjectures, can also be used to explain the behavior of energy landscapes, the speed of folding similar to transition state theory, and the fact that random proteins do not fold.

Multiscale Discrete Approximations of Fourier Integral Operators Associated with Canonical Transformations and Caustics

We develop an algorithm for the computation of general Fourier integral operators associated with canonical graphs. The algorithm is based on dyadic parabolic decomposition using wave packets and enables the discrete approximate evaluation of the action of such operators on data in the presence of caustics. The procedure consists in the construction of a universal operator representation through the introduction of locally singularity-resolving diffeomorphisms, enabling the application of wave packet driven computation, and in the construction of the associated pseudo-differential joint-partition of unity on the canonical graphs. We apply the method to a parametrix of the wave equation in the vicinity of a cusp singularity.

Error estimates for Gaussian beam superpositions

Gaussian beams are asymptotically valid high frequency solutions to hyperbolic partial differential equations, concentrated on a single curve through the physical domain. They can also be extended to some dispersive wave equations, such as the Schrodinger equation. Superpositions of Gaussian beams provide a powerful tool to generate more general high frequency solutions that are not necessarily concentrated on a single curve. This work is concerned with the accuracy of Gaussian beam superpositions in terms of the wavelength epsilon. We present a systematic construction of Gaussian beam superpositions for all strictly hyperbolic and Schrodinger equations subject to highly oscillatory initial data of the form Ae(i Phi/) (epsilon). Through a careful estimate of an oscillatory integral operator, we prove that the k-th order Gaussian beam superposition converges to the original wave field at a rate proportional to epsilon(k/2) in the appropriate norm dictated by the well-posedness estimate. In particular, we prove that the Gaussian beam superposition converges at this rate for the acoustic wave equation in the standard, epsilon-scaled, energy norm and for the Schrodinger equation in the L-2 norm. The obtained results are valid for any number of spatial dimensions and are unaffected by the presence of caustics. We present a numerical study of convergence for the constant coefficient acoustic wave equation in R-2 to analyze the sharpness of the theoretical results.

Differential optics for illumination design in the presence of caustics

Alternatively to using time-consuming Monte-Carlo simulations the irradiation at a target plane can also be calculated by differential optics methods. In the case of caustics, these methods yield to infinite irradiance and its results are not directly comparable to those of Monte-Carlo simulations. In this paper, a differential based algorithm for an on axis point source and a rotationally symmetric optical system is presented, which yields to the same results as a Monte-Carlo simulation. However, the differential optical methods are about three orders of magnitude faster than the latter one, thus allowing fast trial and error design of such kind of illumination systems. An applet is presented that uses sliders to change the shape of the lens and other properties of the illumination system whereas the irradiance profile is nearly immediately perceived. For beginners in the field, this does not only accelerate the design process itself but also the learning process is improved considerably. Some extensions and special cases are shortly discussed.

Extended isochron rays in prestack depth (map) migration

Many processes in seismic data analysis and seismic imaging can be identified with solution operators of evolution equations. These include data downward continuation and velocity continuation. We have addressed the question of whether isochrons defined by imaging operators can be identified with wavefronts of solutions to an evolution equation. Rays associated with this equation then would provide a natural way of implementing prestack map migration. Assuming absence of caustics, we have developed constructive proof of the existence of a Hamiltonian describing propagation of isochrons in the context of common-offset depth migration. In the presence of caustics, one should recast to a sinking-survey migration framework. By manipulating the double-square-root operator, we obtain an evolution equation that describes sinking-survey migration as a propagation in two-way time with surface data being a source function. This formulation can be viewed as an extension of the exploding reflector concept from zero-offset to sinking-survey migration. The corresponding Hamiltonian describes propagation of extended isochrons (fronts with constant two-way time) connected by extended isochron rays. The term extended reflects the fact that two-way time propagation now takes place in high-dimensional space with the following coordinates: subsurface midpoint, subsurface offset, and depth. Extended isochron rays can be used in a natural manner for implementing sinking-survey migration in a map-migration fashion.

Seismic imaging with the generalized Radon transform: a curvelet transform perspective**This research was supported in part under NSF CMG grant DMS-0724808. MVdH was also funded in part by START-grant Y237-N13 of the Austrian Science Fund. HS was supported by NSF grant DMS-0654415. GU was also funded in part by a Walker Family Endowed Professorship.

A key challenge in the seismic imaging of reflectors using surface reflection data is the subsurface illumination produced by a given data set and for a given complexity of the background model (of wave speeds). The imaging is described here by the generalized Radon transform. To address the illumination challenge and enable (accurate) local parameter estimation, we develop a method for partial reconstruction. We make use of the curvelet transform, the structure of the associated matrix representation of the generalized Radon transform, which needs to be extended in the presence of caustics and phase linearization. We pair an image target with partial waveform reflection data, and develop a way to solve the matrix normal equations that connect their curvelet coefficients via diagonal approximation. Moreover, we develop an approximation, reminiscent of Gaussian beams, for the computation of the generalized Radon transform matrix elements only making use of multiplications and convolutions, given the underlying ray geometry; this leads to computational efficiency. Throughout, we exploit the (wave number) multi-scale features of the dyadic parabolic decomposition underlying the curvelet transform and establish approximations that are accurate for sufficiently fine scales. The analysis we develop here has its roots in and represents a unified framework for (double) beamforming and beam-stack imaging, parsimonious pre-stack Kirchhoff migration, pre-stack plane-wave (Kirchhoff) migration and delayed-shot pre-stack migration.

Numerical method for calculating the surface wave field in the presence of caustics

Due to the lateral heterogeneity of the upper layers of the Earth, paths of surface waves deviate from arcs of great circles. Because of the sphericity of the Earth, the paths intersect on a hemisphere opposite to the epicenter and form caustics consisting of two branches, with their tangent point being a cusp. For this reason, the field of surface waves cannot be analyzed in terms of the ray theory at distances larger than 90°. The asymptotic approach to the analysis of the field in the vicinity of such caustics is very ill-suited for numerical implementation. The difficulties of such an approach to the field calculation are aggravated by the fact that such caustics are superimposed in some regions. Therefore, it is suggested to use the theorem of representation, according to which the field within a certain contour is expressed as an integral whose integrand contains values of the function itself, its derivative along the normal to the contour, and Green’s function. The field on the contour (the circle bounding a hemisphere centered at the epicenter) is calculated by the ray method because rays do not intersect on this hemisphere. These data are used for the construction of the field on the opposite hemisphere assumed to be homogeneous, which enables the construction of Green’s function for this hemisphere. This limitation is not very stringent because the configuration of rays and caustics on this hemisphere is mainly determined by the field on the circle. The integral in the representation theorem is calculated numerically. Numerical examples are presented for models in which one caustic or two superimposed caustics form. These calculations yield constraints on variations in the amplitude and phase of the wave. Rayleigh wave fields are also calculated for a model of the real Earth. It is shown that, at some points, the Rayleigh wave spectrum can be strongly distorted because caustics corresponding to different periods differ in shape.

Seismic inverse scattering in the downward continuation approach

Seismic data are commonly modeled by a linearization around a smooth background medium in combination with a high frequency approximation. The perturbation of the medium coefficient is assumed to contain the discontinuities. This leads to two inverse problems, first the linearized inverse problem for the perturbation, and second the estimation of the background, which is a priori unknown (velocity estimation). Here we give a reconstruction formula for the linearized problem using the downward continuation approach. The reconstruction is done microlocally, up to an explicitly given pseudodifferential factor that depends on the aperture. Our main result is a characterization of the wave-equation angle transform, derived from downward continuation, that generates the common image point gathers as an invertible Fourier integral operator, microlocally. We show that the common image point gathers obtained with this particular angle transform are free of so called kinematic artifacts, even in the presence of caustics. The assumption is that the rays in the background that are associated with the reflections due to the medium perturbation are nowhere horizontal. Finally, making use of the mentioned angle transform, pseudodifferential annihilators of the data are constructed. These annihilators detect whether the data are contained in the range of the modeling operator, which is the precise criterion in migration velocity analysis to determine whether a background medium is acceptable, even in the presence of caustics.

Linearized 2.5-dimensional parameter imaging inversion in anisotropic elastic media

SUMMARY In this paper we derive 2.5-D short-period seismic modelling and imaging-inversion formulae in the Born approximation for anisotropic elastic media. The 2.5-D approach encompasses 3-D wave scattering measured in a common-azimuth acquisition geometry subject to 2-D computations under appropriate assumptions. The lowest possible symmetry of the medium in this approach, in principle, is monoclinic, while the medium must be translationally invariant in the normal direction to the associated symmetry plane. In the presence of caustics, artefacts may be generated by the imaging-inversion procedures. We show that in the 2.5-D approach the analysis of artefacts in the 2-D symmetry plane implies the corresponding analysis in 3-D in the framework of the common-azimuth acquisition geometry. An interesting aspect of our results is the occurrence of out-of-plane geometrical spreading in the least-squares removal of the contrast source radiation patterns on the data. We finally introduce the 2.5-D generalized Radon transform that generates common-image-point gathers. The reconsideration of 2.5-D scattering theory is motivated by the increase in use of ocean-bottom acquisition technology. It is not uncommon for ocean bottom cable (OBC) seismic data to be collected along a single line, and the question arises of how to make optimal use of these data. We show the effect of the factors making up the amplitude in the 2.5-D generalized Radon transform on OBC field data from the North Sea.

The applicability of dip moveout/azimuth moveout in the presence of caustics

Reflection seismic data continuation is the computation of data at source and receiver locations that differ from those in the original data, using whatever data are available. We develop a general theory of data continuation in the presence of caustics and illustrate it with three examples: dip moveout (DMO), azimuth moveout (AMO), and offset continuation. This theory does not require knowledge of the reflector positions. We construct the output data set from the input through the composition of three operators: an imaging operator, a modeling operator, and a restriction operator. This results in a single operator that maps directly from the input data to the desired output data. We use the calculus of Fourier integral operators to develop this theory in the presence of caustics. For both DMO and AMO, we compute impulse responses in a constant-velocity model and in a more complicated model in which caustics arise. This analysis reveals errors that can be introduced by assuming, for example, a model with a constant vertical velocity gradient when the true model is laterally heterogeneous. Data continuation uses as input a subset (common offset, common angle) of the available data, which may introduce artifacts in the continued data. One could suppress these artifacts by stacking over a neighborhood of input data (using a small range of offsets or angles, for example). We test data continuation on synthetic data from a model known to generate imaging artifacts. We show that stacking over input scattering angles suppresses artifacts in the continued data.