Seismic simulation methods can be broadly categorized into two groups: (a) direct solving the wave equation using numerical algorithms, such as finite-differences (Hestholm & Ruud, 1998;Pitarka, 1999), finite-elements (Basabe & Sen, 2009), pseudo-spectral (Huang, 1992 and spectral-element approaches (Komatitsch & Tromp, 1999), and (b) the high-frequency asymptotic solution of the wave equation (Bullen, 1961;Červený, 2005). All these methods have their own advantages and drawbacks, and all are useful for specific seismic problems.Ray theory is representative of the high-frequency asymptotic methods for modeling seismic wavefields by approximating the harmonic solution of the wave equation using an infinite asymptotic series (Červený, 2005;Popov, 2002). Thus, it is also known as asymptotic ray theory. The high-frequency approximation simplifies the wave equation into an eikonal equation and a transport equation, which determine the traveltimes and amplitudes of seismic wavefields, respectively. The characteristics of the eikonal equation define seismic rays; the traveltime Abstract Ray theory gives a high-frequency asymptotic solution of the wave equation, which has been widely used in the seismic study because of its high computational efficiency and good physical insights for elementary waves. However, the fundamental high-frequency assumption makes it difficult to accurately simulate the finite-frequency effect of wave propagations. On the other hand, the Gaussian beam, as one of the advanced ray-based methods, overcomes the amplitude singularity problem near the caustics by introducing a complex-valued paraxial approximation. But Gaussian beams only use the velocity and up to second-order velocity derivative of central rays and thus does not accurately simulate the response of heterogeneity away from the rays. To bridge the gap between the ray theory and wave-equation methods, we present a two-way beam wave method and apply it to seismic imaging. A fan of central rays are first shot to extract local ray-centered model parameters in the global Cartesian coordinates. Then, a finite-difference method is used to solve the acoustic wave equation in the ray-centered coordinates to simulate wave propagations near the central rays. The beam wave method can accurately simulate the response of the velocity heterogeneities in the ray tubes and therefore produces more accurate wavefields. In addition, we discretize the ray-centered wave equation onto a non-orthogonal grid, which considerably reduces the computational cost. We apply the two-way beam wave method in seismic imaging and develop a beam wave reverse-time migration. It inherits the advantages of both ray-based and wave equation migrations, and produces clear images for complicated structures with fewer artifacts than traditional imaging methods.Plain Language Summary Seismic ray theory has evolved from classical asymptotic ray theory, through paraxial ray theory and Maslov's method, and then to the Gaussian beam method. Although the computational accuracy of the ...