Prevailing-frequency approximation of the coupling ray theory for electromagnetic waves or elastic S waves

The coupling–ray–theory tensor Green function for electromagnetic waves or elastic S waves is frequency dependent, and is usually calculated for many frequencies. This frequency dependence represents no problem in calculating the Green function, but may represent a great problem in storing the Green function at the nodes of dense grids, typical for applications such as the Born approximation. This paper is devoted to the approximation of the coupling–ray–theory tensor Green function, which practically eliminates this frequency dependence within a reasonably broad frequency band. In the vicinity of a given prevailing frequency, we approximate the frequency–dependent frequency–domain coupling–ray–theory tensor Green function by two dyadic Green functions corresponding to two waves described by their travel times and amplitudes calculated for the prevailing frequency. We refer to these travel times and amplitudes as the coupling–ray–theory travel times and the coupling–ray–theory amplitudes. This “prevailing–frequency approximation” of the coupling ray theory for electromagnetic waves or elastic S waves allows us to process the coupling–ray–theory wave field in the same way as the anisotropic–ray–theory wave field. This simplification may be decisive when storing the tensor Green function at the nodes of dense grids, which is typical for applications such as the Born approximation. We test the accuracy of the proposed prevailing–frequency approximation of the coupling ray theory numerically using elastic S waves in eight anisotropic velocity models. The additional inaccuracy introduced by the prevailing–frequency approximation is smaller than the inaccuracy of the standard frequency–domain coupling ray theory, and smaller than the additional inaccuracy introduced by many other approximations of the coupling ray theory.


INTRODUCTION
There are two different high-frequency asymptotic ray theories for electromagnetic waves or elastic S waves with frequency-independent amplitudes: the isotropic ray theory based on the assumption of equal velocities of both electromagnetic waves or both elastic S waves, and the anisotropic ray theory assuming both waves are strictly decoupled. Here the term "different" means that the isotropic ray theory is not a special case of the anisotropic ray theory for decreasing anisotropy, and that both theories yield different electromagnetic waves or elastic S waves in equal media.
This paper is equally applicable to electromagnetic waves and elastic S waves. For electromagnetic waves, we consider the ray theory developed for the magnetic vector potential with the Weyl gauge (zero electric potential) according to Klimeš (2010aKlimeš ( , 2016a. In this case, the Kelvin-Christoffel matrix is a 3×3 matrix, with two eigenvalues which may equal zero. These two eigenvalues correspond to the two polarizations of electromagnetic waves. Hereinafter we shall refer to these eigenvalues as the first two eigenvalues, and to the corresponding eigenvectors as the first two eigenvectors. The third eigenvalue of the 3 × 3 Kelvin-Christoffel matrix for electromagnetic waves cannot equal zero and does not correspond to any wave. We shall refer to this eigenvalue as the third eigenvalue, and to the corresponding eigenvector as the third eigenvector.
For elastic waves, the Kelvin-Christoffel matrix is represented by the difference of the 3×3 Christoffel matrix and the 3×3 identity matrix. We shall refer to the Swave eigenvalues as the first two eigenvalues, and to the corresponding eigenvectors as the first two eigenvectors. We shall refer to the P-wave eigenvalue as the third eigenvalue, and to the corresponding eigenvector as the third eigenvector. In the elastic case, we consider just S waves, hereinafter referred to simply as waves.
In the isotropic ray theory, the polarization vectors of waves do not rotate about the third eigenvector of the Kelvin-Christoffel matrix, whereas in the anisotropic ray theory they coincide with the first two eigenvectors of the Kelvin-Christoffel matrix which may rotate rapidly. Thomson et al. (1992) demonstrated analytically that the high-frequency asymptotic error of the anisotropic-ray-theory wave field is inversely proportional to the second or higher root of the frequency, if the ray passes through the point of the first two equal eigenvalues of the Kelvin-Christoffel matrix even in otherwise strongly anisotropic media.
In "weakly anisotropic" media, at moderate frequencies, the actual wave polarization tends to remain unrotated about the third eigenvector of the Kelvin-Christoffel matrix, but is partly attracted by the rotation of the first two eigenvectors of the Christoffel matrix. The intensity of the attraction increases with frequency. This behaviour of the actual wave polarization is described by the coupling ray theory proposed, e.g., by Kravtsov (1968), Naida (1977Naida ( , 1979 or Fuki et al. (1998) for electromagnetic waves, and by Coates and Chapman (1990) for elastic S waves. The frequency-dependent coupling ray theory is the generalization of both the zero-order isotropic and anisotropic ray theories and provides continuous transition between them. The coupling ray theory is applicable to coupled waves at all degrees of anisotropy, from isotropic to considerably anisotropic media. The numerical algorithm for calculating the frequency-dependent coupling-ray-theory tensor Green function has been designed by Bulant and Klimeš (2002).
The coupling ray theory can be derived from the asymptotic ray series referred to as the coupling ray series (Klimeš, 2013 ). This derivation demonstrates that there are many possible forms of the coupling ray theory of varying accuracy, and that the currently used coupling ray theory need not represent the most accurate option. Moreover, there are many more or less accurate approximations of this coupling ray theory (Bulant and Klimeš, 2002 ;Pšenčík, 2008, 2010 ;Pšenčík et al., 2012 ). Note that the anisotropic-common-ray approximation (Klimeš, 2006 ;Klimeš and Bulant, 2006 ) usually preserves the accuracy of the coupling ray theory and simultaneously simplifies ray tracing considerably by eliminating the slowness surface singularities. However, in approximately uniaxial (approximately transversely isotropic) media (Klimeš, 2015(Klimeš, , 2016c, we may also use the reference ordinary (SH) and extraordinary (SV) rays (Klimeš and Bulant, 2015 ).
The coupling-ray-theory tensor Green function is frequency dependent, and is usually calculated for many frequencies. This frequency dependence represents no problem in calculating the Green function, but may represent a great problem in storing the Green function at the nodes of dense grids (Klimeš and Bulant, 2013 ), typical for applications such as the Born approximation. This paper is devoted to the approximation of the coupling-ray-theory tensor Green function, which practically eliminates this frequency dependence within a limited frequency band.
In this paper, we approximate the frequency-dependent frequency-domain coupling-ray-theory tensor Green function in the vicinity of a given prevailing frequency by two dyadic Green functions corresponding to the two coupled waves described by their travel times and amplitudes calculated for the prevailing frequency. We refer to these travel times and amplitudes as the coupling-ray-theory travel times and the coupling-ray-theory amplitudes. The presented numerical examples suggest that the coupling-ray-theory travel times are usually close to the anisotropic-raytheory travel times. However, if the rays pass close to a conical or similar slowness surface singularity, the singularity acts as an interface and smoothly but very rapidly converts the actual wave polarizations from the approximately first anisotropicray-theory polarization (S1) to the approximately second anisotropic-ray-theory polarization (S2), and vice versa. In this case, the coupling-ray-theory travel times correspond approximately to the anisotropic-ray-theory travel times of the S1S2 and S2S1 waves converted at the slowness surface singularity.

COUPLING RAY THEORY
For electromagnetic waves, the time-domain tensor Green function G ij (x, x 0 , t) represents the i th component of the magnetic vector potential at point x due to the unit electric current at point x 0 oriented along the j th coordinate axis and having the Dirac-distribution time dependence (Klimeš, 2016a). For elastic waves, the time-domain tensor Green function G ij (x, x 0 , t) represents the i th displacement component at point x due to the unit force at point x 0 oriented along the j th coordinate axis and having the Dirac-distribution time dependence (Červený, 2001, Eq. 2.5.37 ).
The coupling-ray-theory approximation of its Fourier transform G ij (x, x 0 , ω) is given in Section 2.1. It is calculated along reference rays and is expressed in terms of the 2×2 complex-valued coupling-equation propagator matrix Π KL (x, x 0 , ω) which is given in Section 2.2. are orthogonal and if the numerical algorithm used to integrate coupling equation (3) is stable and accurate with respect to the singular behaviour of ϕ ′ , see Bulant and Klimeš (2002).

PREVAILING-FREQUENCY APPROXIMATION OF THE COUPLING RAY THEORY
Analogously to the anisotropic ray theory, we wish to approximate the frequencydependent coupling-ray-theory tensor Green function by two dyadic Green functions corresponding to two waves described by their travel times and amplitudes. In the vicinity of the prevailing circular frequency ω 0 , we thus wish to approximate propagator matrix Π = Π(ω) by matrix where This decomposition represents the approximation of propagator matrix Π = Π(ω) by two waves with travel times τ − D and τ + D.
In the Appendix, we demonstrate that assumption (12) does not contradict the above requirements (9) and (10).
3 . 1 . D e t e r m i n i n g t h e c o e f f i c i e n t s o f p r e v a i l i n g -f r e q u e n c y a p p r o x i m a t i o n ( 7 ) We now determine coefficients D, Π (1) and Π (2) of the prevailing-frequency approximation (7) of the coupling ray theory from conditions (9), (10) and (12).
We differentiate definition (7) and obtain Equations (9) and (10) then read and All 2×2 matrices satisfying identities (12) satisfy relation Since matrix Π(ω 0 ) in relation (14) is unimodular, The determinants of both sides of matrix relation (15) then yield The right-hand side of this equation is real-valued and non-negative, see the elements of matrix (A-8). Without loss of generality, we may choose non-negative solution Equations (14) and (15) then yield and The prevailing-frequency approximation (7) of the coupling ray theory is determined by coefficients (19)-(21). The right-hand sides of Eqs (20) and (21) are well defined even for very small D, see matrix (A-8) and definition (19). If D = 0, matrices Π (1) and Π (2) satisfying condition (14) may be chosen arbitrarily. Note that if we insert matrices (20) and (21) into the approximate propagator matrix (7), we obtain relation However, here we prefer decomposition (7) into two waves over mixed expression (22).
3 . 2 . P r e v a i l i n g -f r e q u e n c y a p p r o x i m a t i o n o f t h e t e n s o r G r e e n f u n c t i o n We now approximate propagator matrix Π KL (x, x 0 , ω) in tensor Green function (1) by prevailing-frequency approximation (7). After this insertion, the couplingray-theory tensor Green function (1) may be approximated by the superposition of two dyadic Green functions The corresponding two waves have coupling-ray-theory travel times and The complex-valued vectorial amplitudes of these two waves are and Whereas propagator matrix Π = Π(ω) in tensor Green function (1) is strongly frequency dependent at high frequencies, the frequency dependence of the couplingray-theory travel times (25)-(26) and corresponding amplitudes (27) . D i f f e r e n c e b e t w e e n t h e c o u p l i n g -r a y -t h e o r y p o l a r i z a t i o n v e c t o r s a n d t h e a n i s o t r o p i c -r a y -t h e o r y p o l a r i z a t i o n v e c t o r s The complex-valued coupling-ray-theory polarization vectors determined by complex-valued vectorial amplitudes (27) and (28) may differ considerably from both the isotropic-ray-theory polarization vectors and the anisotropic-ray-theory polarization vectors, depending on the heterogeneity and on the degree of anisotropy or bianisotropy. The possible differences in the polarization of elastic S waves are well visible in three-component synthetic seismograms.

ELASTIC S-WAVE NUMERICAL EXAMPLES
The prevailing-frequency approximation of the coupling ray theory for elastic S waves has been coded as a new option of program green.for of package CRT (Bucha and Bulant, 2012 ) using the numerical algorithm described in Section 4. Program green.for now allows synthetic seismograms of coupled elastic S-waves in heterogeneous weakly anisotropic velocity models to be calculated using either the standard frequency-domain coupling ray theory with propagator matrix (29) calculated for all frequencies, or its prevailing-frequency approximation according to Section 4. In this paper, the prevailing-frequency approximation of the coupling ray theory is numerically tested in eight weakly anisotropic velocity models. These anisotropic velocity models with related input data can be found in package DATA (Bucha and Bulant, 2012 ).

. 1 . R e f e r e n c e a n d p e r t u r b a t i o n H a m i l t o n i a n f u n c t i o n s
The reference ray is the same for the standard frequency-domain coupling ray theory and for the prevailing-frequency approximation. The reference ray for the coupling equation is calculated using the reference Hamiltonian function. We consider reference Hamiltonian function The resulting anisotropic common reference ray depends on exponent M which controls the averaging of the S-wave eigenvalues G 1 and G 2 of the Christoffel matrix. The reference ray is independent of exponent N which controls the perturbation approximation of anisotropic-ray-theory travel times τ (1) and τ (2) . Refer to Klimeš (2006) for the calculation of the reference ray and other reference quantities.
Anisotropic-ray-theory travel times τ (1) and τ (2) with corresponding S-wave eigenvectors g (1) i and g (2) i of the Christoffel matrix are usually approximated by various perturbation expansions along the reference rays (Klimeš, 2002(Klimeš, , 2010b(Klimeš, , 2016b. The perturbation approximation of the anisotropic-ray-theory travel times depends on the reference ray and on degree N of the perturbation Hamiltonian function homogeneous with respect to the slowness vector. We choose perturbation Hamiltonian function which is linear with respect to perturbation parameters f 1 and f 2 , because we believe that this linearity results in a more accurate perturbation expansion than the expansions resulting from nonlinear Hamiltonian functions, but other authors may choose other perturbation Hamiltonian functions. The same exponent N must be used in definitions (41) and (42). Different perturbation Hamiltonian functions yield different anisotropic-ray-theory travel times τ (1) , τ (2) and Christoffel matrices of varying accuracy with their eigenvectors g i . The perturbation usually yields the best results for the homogeneous perturbation Hamiltonian function of degree N = −1 with respect to the slowness vector, which was theoretically explained by Klimeš (2002, Sec. 4.4), and numerically demonstrated by Vavryčuk (2012) in examples of perturbations from real-valued reference rays to the complex-valued travel time in two isotropic attenuating media.
To calculate the reference ray, we use the reference Hamiltonian function, which is the average of the homogeneous Hamiltonian functions of degree M = −1 with respect to the slowness vector. Note that the perturbation expansions along reference rays, calculated with M = −1 and M = 2, have been compared for N = −1 by Bulant and Klimeš (2008 , Tables 7-9 and 10-12).

. 2 . N u m e r i c a l m e t h o d s b e i n g c o m p a r e d
In the numerical examples presented here, we follow the work of Pšenčík et al. (2012), who compared the synthetic seismograms calculated by their approximation of the coupling ray theory with the synthetic seismograms generated by the Fourier pseudospectral method. Bulant et al. (2011) then numerically compared three different approximations of the coupling ray theory with the Fourier pseudospectral method.
In this paper, we compare the synthetic seismograms calculated by three numerical methods: (a) the proposed prevailing-frequency approximation of the coupling ray theory; (b) the standard coupling ray theory calculated for all frequencies according to the algorithm of Bulant and Klimeš (2002); (c) the Fourier pseudospectral method which is considered here as a nearly exact reference. Refer to Pšenčík et al. (2012) for the description of the Fourier pseudospectral method used here.
Both the coupling ray theory and its prevailing-frequency approximation are calculated along anisotropic common reference rays (Klimeš, 2006 ). To calculate the reference rays, we use reference Hamiltonian function (41) which is the average of the homogeneous Hamiltonian functions of degree M = −1 with respect to the slowness vector. For perturbations, we use perturbation Hamiltonian function (42) which is a homogeneous function of degree N = −1 with respect to the slowness vector. We approximate the anisotropic-ray-theory travel times τ (1) and τ (2) in the coupling equation by the second-order perturbation expansion according to Klimeš and Bulant (2006). We approximate eigenvectors g (1) i and g (2) i of the Christoffel matrix by the unperturbed reference eigenvectors calculated using the Christoffel matrix with the unperturbed reference slowness vector.

. 3 . A n i s o t r o p i c v e l o c i t y m o d e l s
We consider eight weakly anisotropic velocity models referred to as QIH, QI, QI2, QI4, SC1 I, SC1 II, KISS and ORT. All these velocity models are laterally homogeneous. The density-reduced elastic moduli are linear functions of depth in all these velocity models. The density is constant.
A vertically heterogeneous 1-D anisotropic velocity model QI was provided by Pšenčík and Dellinger (2001, model WA rotated by 45 • ), who performed the coupling-ray-theory calculations using the programs of package ANRAY (Pšenčík, 1998 ) and compared the results with the reflectivity method. The same velocity model has been used by Bulant and Klimeš (2002) and Klimeš and Bulant (2004) to demonstrate the coupling ray theory and its quasi-isotropic approximations. For comparison with the isotropic-ray-theory and anisotropic-ray-theory seismograms in velocity model QI and for a more detailed discussion and description of this velocity model refer to Pšenčík and Dellinger (2001). Velocity model QI is approximately transversely isotropic (Klimeš, 2015(Klimeš, , 2016c) in a vertical plane which forms a 45 • angle with the source-receiver vertical plane and is situated between the positive radial and positive transverse seismogram components. The reference symmetry axis is thus horizontal and is situated right in the middle between the positive radial and negative transverse seismogram components.
Velocity models QIH, QI2 and QI4 are derived from velocity model QI and mutually differ by their degrees of anisotropy. The differences of the elastic moduli in velocity models QIH, QI, QI2 and QI4 from the elastic moduli in the reference isotropic velocity model are determined by ratio 0.5 : 1 : 2 : 4. For the elastic moduli in velocity models QI, QI2 and QI4 refer to Bulant and Klimeš (2008).
For the description of velocity models SC1 I, SC1 II, KISS and ORT and for the comparison with the anisotropic-ray-theory seismograms in these velocity models refer to Pšenčík et al. (2012).
At the depth of 0 km, velocity model SC1 I is approximately transversely isotropic and its reference symmetry axis is horizontal (Klimeš, 2015(Klimeš, , 2016c. At this depth, the slowness surface contains a split intersection singularity, whereas velocity models QIH, QI, QI2 and QI4 display no exact nor split intersection singularity (Pšenčík et al., 2012 ). At the depth of 1.5 km, velocity model SC1 I is very close to isotropic, but is slightly cubic and its symmetry axes coincide with the coordinate axes. This means that, at depths between 0 km and 1.5 km, velocity model SC1 I is very close to transversely isotropic, but is slightly tetragonal. Whereas the transversely isotropic medium contains the intersection singularity through which Stud. Geophys. Geod., 60 (2016) the rays pass without rotation of the eigenvectors of the Christoffel matrix (Vavryčuk, 2001, Sec. 4.3 ;Klimeš and Bulant, 2014a), in the approximately transversely isotropic medium, the S-wave slowness surface is split at this unstable singularity (Crampin, 1981 ) and the eigenvectors of the Christoffel matrix rapidly rotate by 90 • there. This split intersection singularity then acts as an interface and smoothly but very rapidly converts the actual S-wave polarizations from the approximately anisotropic-ray-theory polarization S1 to the approximately anisotropic-ray-theory polarization S2, and vice versa.
Note that the anisotropic-ray-theory geometrical spreading and the phase shift due to caustics cannot often be calculated numerically through the split intersection singularity (Bulant and Klimeš, 2014 ).
Velocity model SC1 II is analogous to SC1 I, but the reference symmetry axis of its approximately transversely isotropic component is tilted, refer to Pšenčík et al. (2012). The symmetry axes of its weak cubic component coincide with the coordinate axes. The split intersection singularity in velocity model SC1 II is thus positioned differently in comparison with velocity model SC1 I. In the source-receiver plane, the split intersection singularity is close to the horizontal slowness vectors. When the slowness vector of a ray passes smoothly through this split intersection singularity, the ray-velocity vector rapidly changes its direction and creates a sharp bend in the anisotropic-ray-theory ray. The deviation of the anisotropic-ray-theory rays from the anisotropic common reference rays is much greater in velocity model SC1 II than in velocity model SC1 I (Pšenčík et al., 2012 ). However, the actual wave paths in velocity model SC1 II are considerably different from both the anisotropic common reference rays and the anisotropic-ray-theory rays. The actual wave paths in velocity model SC1 II are close to the SH and SV reference rays defined by Klimeš and Bulant (2015).
Velocity model KISS represents velocity model QI described above, rotated by −44 • in order to position the reference symmetry axis, corresponding to the kiss slowness surface singularity, just 1 • from the source-receiver plane.
In the weakly orthorhombic velocity model ORT, the slowness surface contains four conical singularities. The rays leading from the source to the middle part of the receiver profile pass close to one of these singularities. This conical singularity then acts as an interface and smoothly but very rapidly converts the actual S-wave polarizations from the approximately anisotropic-ray-theory polarization S1 to the approximately anisotropic-ray-theory polarization S2, and vice versa.

. 4 . M e a s u r e m e n t c o n f i g u r a t i o n
The synthetic seismograms generated by a vertical force are calculated at the receivers located in a vertical well at a distance of 1 km from the source. The sourcereceiver configuration in velocity models QIH, QI, QI2, QI4 and KISS is displayed in Fig. 1. The source-receiver configuration in velocity models SC1 I, SC1 II and ORT is similar except for different receiver positions in the vertical well.
The source time function is the Gabor signal cos(2πf t) exp[−(2πf t/4) 2 ] with reference frequency f = 50 Hz, bandpass filtered by a cosine filter specified by For the prevailing-frequency approximation of the coupling ray theory, we naturally use the prevailing frequency f 0 = 50 Hz.

. . C o m p a r i s o n o f n a r r o w -b a n d s y n t h e t i c s e i s m o g r a m s
The synthetic seismograms in velocity models QIH, QI, QI2, QI4, SC1 I, SC1 II, KISS and ORT are displayed in Figs 2-9. In each figure, the synthetic seismograms calculated by the standard coupling ray theory and by the Fourier pseudospectral method are displayed in different colours. In each velocity model, the amplitude scaling is equal for all components and all receivers.
In all eight velocity models, the additional inaccuracy introduced by the prevailing-frequency approximation is generally smaller than the inaccuracy of the standard frequency-domain coupling ray theory. The synthetic seismograms calculated by the prevailing-frequency approximation of the coupling ray theory are plotted in red, and are often obscured by the standard-coupling-ray-theory seismograms plotted in green. The additional inaccuracy introduced by the prevailing-frequency approximation is most pronounced in velocity model QI2 (Fig. 4) and is mostly invisible or even completely invisible in the other seven velocity models.
In velocity model QIH (Fig. 2), we have only the results of the standard coupling ray theory and of its prevailing-frequency approximation. We can observe very good agreement between the prevailing-frequency approximation seismograms and the standard-coupling-ray-theory seismograms. The dependence between the strength of the anisotropy and the differences of both the coupling-ray-theory seismograms from the Fourier pseudospectral seismograms can be observed in Figs 3-5. From this dependence, we may estimate that the accuracy of both the coupling-ray-theory seismograms in velocity model QIH is very good.
In velocity model QI (Fig. 3) whose anisotropy is twice stronger than in model QIH, we can again observe the good overall fit between the prevailing-frequency approximation seismograms and the standard-coupling-ray-theory seismograms. Both these coupling-ray-theory seismograms differ from the Fourier pseudospectral seismograms only negligibly.
In velocity model QI2 (Fig. 4) whose anisotropy is twice stronger than in model QI, we can observe slight differences between the prevailing-frequency approximation seismograms and the standard-coupling-ray-theory seismograms in the second (transverse) component. This component also displays the development of S-wave splitting. We may thus conclude that the coupling of the two S waves plays an important role in this velocity model and is also related to the slight inaccuracy of the prevailing-frequency approximation. However, the additional inaccuracy, introduced by the prevailing-frequency approximation, is smaller than the inaccuracy of the standard coupling ray theory. Note that the differences of both the couplingray-theory seismograms from the Fourier pseudospectral seismograms increase with increasing anisotropy, see Figs 3-5.
In velocity model QI4 (Fig. 5) whose anisotropy is 4 times stronger than in model QI, we can see that both the prevailing-frequency approximation and the standard coupling ray theory perform equally well, even though the velocity model is considerably anisotropic, and the separation of the two S waves is about 0.06 s which corresponds approximately to three S-wave periods.
The slowness surface in velocity model SC1 I contains the split intersection singularity. In Fig. 6, the prevailing-frequency approximation seismograms are mostly obscured by the standard coupling ray theory seismograms. Both the coupling-raytheory seismograms are close to the Fourier pseudospectral seismograms despite the existence of the above-mentioned split intersection singularity in this velocity model.
The slowness surface in velocity model SC1 II contains the split intersection singularity analogous to velocity model SC1 I, but with a tilted reference symmetry axis. In the source-receiver plane, the split intersection singularity is close to the horizontal slowness vectors. When the slowness vector of a ray passes smoothly through this split intersection singularity, the ray-velocity vector rapidly changes its direction and creates a sharp bend of the anisotropic-ray-theory ray (Bulant and Klimeš, 2014 ;Klimeš and Bulant, 2014b). This sharp bend is connected with a rapid rotation of the eigenvectors of the Christoffel matrix. The sharply bent anisotropicray-theory rays thus cannot describe the correct wave propagation. The actual wave paths in velocity model SC1 II are close to the SH and SV reference rays (Klimeš, 2015(Klimeš, , 2016cKlimeš and Bulant, 2015 ) Fig. 2. Radial, transverse and vertical components of the seismograms calculated in velocity model QIH. All seismograms are scaled equally. In this velocity model, we have only the results of the standard coupling ray theory and of its prevailing-frequency approximation. The seismograms calculated by the prevailing-frequency approximation are plotted in red, and they are overlaid by the seismograms calculated by the standardcoupling-ray-theory seismograms plotted in green. The seismograms are in good agreement; the prevailing-frequency approximation seismograms are obscured by the standard-coupling-ray-theory seismograms. intersection singularity. The actual wave paths leading to the shallow receivers are significantly different from the anisotropic common reference rays. The differences between the Fourier pseudospectral method and the coupling ray theory at the shallow receivers in Fig. 7 are mostly caused by the inaccurate reference geometrical spreading and by the deviation of the actual polarization vectors of the two S waves from the reference polarization plane given by the reference S-wave eigenvectors corresponding to the reference Christoffel matrix (Klimeš and Bulant, 2015 ). The prevailing-frequency approximation is nearly identical to the standard coupling ray theory even in this situation.
The slowness surface in velocity model KISS contains a kiss singularity analogous to velocity model QI, but deviated by only 1 • from the source-receiver plane, instead of 45 • in velocity model QI. As a consequence, the vertical force generates almost no energy in the second (transverse) component. Fig. 8 shows that the differences of both the coupling-ray-theory methods from the Fourier pseudospectral method are negligible, and that the differences of the prevailing-frequency approximation from the standard coupling ray theory are invisible.
The slowness surface in velocity model ORT contains four conical singularities. The rays leading from the source to the receivers at depths from 0.72 km to 0.92 km pass close to one of these singularities. This conical singularity then acts as an interface and smoothly but very rapidly converts the actual S-wave polarizations from the approximately anisotropic-ray-theory polarization S1 to the approximately anisotropic-ray-theory polarization S2, and vice versa. In Fig. 9, we can see the very good overall fit of both the prevailing-frequency approximation seismograms and the standard-coupling-ray-theory seismograms with the Fourier pseudospectral seismograms, except for the slight differences in the singular region around the depth of 0.8 km. The differences of the prevailing-frequency approximation from the standard coupling ray theory are invisible.

. 6 . C o m p a r i s o n o f b r o a d -b a n d s y n t h e t i c s e i s m o g r a m s
In the previous section, we tested the prevailing-frequency approximation of the coupling ray theory in eight velocity models for narrow-band synthetic seismograms. The source time function was the Gabor signal cos(2πf t) exp[−(2πf t/4) 2 ] with reference frequency f = 50 Hz, bandpass filtered by a cosine filter specified by frequencies 0 Hz, 5 Hz, 60 Hz and 100 Hz. We observed the most pronounced differences between the prevailing-frequency approximation and the standard coupling ray theory in velocity model QI2 (Fig. 4), in which the frequency-dependent coupling of the two S waves is strong.
We thus test the prevailing-frequency approximation of the coupling ray theory in velocity model QI2 also for broad-band synthetic seismograms. The source time function is the Gabor signal cos(2πf t) exp[−(2πf t) 2 ] with reference frequency f = 50 Hz, bandpass filtered by a cosine filter specified by frequencies 0 Hz, 5 Hz, 200 Hz and 250 Hz.
In this case, we have the results of the standard coupling ray theory and of its prevailing-frequency approximation only. The synthetic seismograms are displayed  Fig. 3. Radial, transverse and vertical components of the seismograms calculated in velocity model QI. The results of the prevailing-frequency approximation are plotted in red, the standard-coupling-ray-theory seismograms are then plotted in green, and they are overlaid by the seismograms calculated by the Fourier pseudospectral method plotted in black. All seismograms are in good agreement, with the prevailing-frequency approximation seismograms obscured by the standard-coupling-ray-theory seismograms. The wave arriving at the shallowest receivers at about time 0.53 s is an artifact of the Fourier pseudospectral method.  Fig. 4. Radial, transverse and vertical components of the seismograms calculated in velocity model QI2. The seismograms are plotted in order red (prevailing-frequency approximation), green (standard coupling ray theory) and black (Fourier pseudospectral method). All seismograms are in overall good agreement, with the prevailing-frequency approximation seismograms mostly obscured by the standard-coupling-ray-theory seismograms. Slight differences between the prevailing-frequency approximation seismograms and the standard-coupling-ray-theory seismograms appear in the second (transverse) component. The second component also displays the development of S-wave splitting. We may thus conclude that the frequency-dependent coupling of the two S waves plays an important role in this velocity model and is also related to the slight inaccuracy of the prevailing-frequency approximation.  The seismograms are plotted in order red (prevailing-frequency approximation), green (standard coupling ray theory) and black (Fourier pseudospectral method). The prevailing-frequency approximation seismograms are mostly obscured by the standard-coupling-ray-theory seismograms. The velocity model is considerably anisotropic and the two S waves are fully separated. The coupling effects are thus less pronounced, and the prevailing-frequency approximation seismograms are in good agreement with the standard-coupling-ray-theory seismograms. The differences between both the couplingray-theory seismograms and the Fourier pseudospectral method are slightly higher than in the previous velocity models.  Fig. 6. Radial, transverse and vertical components of the seismograms calculated in velocity model SC1 I. The seismograms are plotted in order red (prevailing-frequency approximation), green (standard coupling ray theory) and black (Fourier pseudospectral method). The prevailing-frequency approximation seismograms are mostly obscured by the standard-coupling-ray-theory seismograms. Both the coupling-ray-theory seismograms are close to the Fourier pseudospectral seismograms.  Fig. 7. Radial, transverse and vertical components of the seismograms calculated in velocity model SC1 II. The seismograms are plotted in order red (prevailing-frequency approximation), green (standard coupling ray theory) and black (Fourier pseudospectral method). The prevailing-frequency approximation seismograms are mostly obscured by the standard-coupling-ray-theory seismograms. The seismograms from the shallow receivers indicate problems with the inaccurate reference polarization vectors and the inaccurate reference geometrical spreading. However, the prevailing-frequency approximation represents a very good approximation to the standard coupling ray theory even in this situation.  The seismograms are plotted in order red (prevailing-frequency approximation), green (standard coupling ray theory) and black (Fourier pseudospectral method). The prevailing-frequency approximation seismograms are obscured by the standard-coupling-ray-theory seismograms. All seismograms are in good agreement, similarly as in Fig. 3 Fig. 9. Radial, transverse and vertical components of the seismograms calculated in velocity model ORT. The seismograms are plotted in order red (prevailing-frequency approximation), green (standard coupling ray theory) and black (Fourier pseudospectral method). The prevailing-frequency approximation seismograms are obscured by the standard-coupling-ray-theory seismograms. Both the coupling-ray-theory seismograms are close to the Fourier pseudospectral seismograms despite the existence of conical singularities in this velocity model.  Fig. 10.Radial, transverse and vertical components of the seismograms calculated in velocity model QI2 for a broad-band signal. In this case, we have only the results of the standard coupling ray theory and of its prevailing-frequency approximation. The seismograms are plotted in order red (prevailing-frequency approximation) and green (standard coupling ray theory). The differences in the polarization due to the strong frequency-dependent coupling of the two S waves are visible in the second (transverse) component.
in Fig. 10. The differences visible in the second (transverse) component are mostly caused by the inaccurate approximation of the coupling-ray-theory seismograms by the prevailing frequency approximation at low frequencies, for which the couplingray-theory polarization approaches the isotropic-ray-theory polarization. There are also differences in the second (transverse) component at the shallow receivers caused by inaccurate approximation of the coupling-ray-theory seismograms by the prevailing frequency approximation at high frequencies. These high-frequency differences vanish at the deep receivers, where the two S waves are split at the prevailing frequency (Fig. 4).

. 7 . C o u p l i n g -r a y -t h e o r y t r a v e l t i m e s
The coupling-ray-theory travel times (25) and (26) are often close to the anisotropic-ray-theory travel times τ (1) and τ (2) . The most important exception is observed in velocity model ORT, at the receivers at depths from 0.72 km to 0.92 km. The rays leading from the source to these receivers pass close to the conical singularity at the slowness surface, and the eigenvectors of the Christoffel matrix smoothly but very rapidly rotate by 90 • (Vavryčuk, 2005, Fig. 1 ). This conical singularity then acts as an interface and smoothly but very rapidly converts the actual S-wave polarizations from the approximately anisotropic-ray-theory polarization S1 to the approximately anisotropic-ray-theory polarization S2, and vice versa. The coupling-ray-theory travel times then correspond approximately to the S1S2 and S2S1 anisotropic-ray-theory travel times of the converted waves, rather than to the S1 and S2 anisotropic-ray-theory travel times. This behaviour of the coupling-ray-theory travel times is evident in Table 1, where the differences of the Table 1. Receivers in velocity model ORT influenced by the conical singularity at the slowness surface. Differences D of the coupling-ray-theory travel times from the mean travel time in seconds at the frequency of 50 Hz are compared with the differences 1 2 |τ (2) −τ (1) | of the anisotropic-ray-theory travel times from the mean travel time in seconds. The anisotropic-ray-theory travel times are calculated by the second-order perturbation expansion along the reference rays traced in steps of 0.01 s. The unknown error of the perturbation expansion cannot influence the comparison because differences D are calculated using differences 1 2 |τ (2) −τ (1) |. The rounding errors should not influence the displayed differences. coupling-ray-theory travel times (25) and (26) from the mean travel time (2) at the frequency of 50 Hz are compared with the differences of the anisotropic-ray-theory travel times τ (1) and τ (2) from the mean travel time (2). The coupling-ray-theory travel times at these receivers are practically frequency independent between 10 Hz and 100 Hz. Analogous behaviour of the coupling-ray-theory travel times is also observed in velocity model SC1 I, at the deepest receivers at depths from 1.12 km to 1.32 km. The rays leading from the source to these receivers pass through the split intersection singularity at the slowness surface, and the eigenvectors of the Christoffel matrix smoothly but very rapidly rotate by 90 • . This split intersection singularity then acts as an interface and smoothly but very rapidly converts the actual S-wave polarizations from the approximately anisotropic-ray-theory polarization S1 to the approximately anisotropic-ray-theory polarization S2, and vice versa. The coupling-ray-theory travel times then correspond approximately to the S1S2 and S2S1 anisotropic-ray-theory travel times of the converted waves, rather than to the S1 and S2 anisotropic-ray-theory travel times. This behaviour of the coupling-raytheory travel times is evident in Table 2, where the differences of the coupling-raytheory travel times (25) and (26) from the mean travel time (2) at the frequency of 50 Hz are compared with the differences of the anisotropic-ray-theory travel times τ (1) and τ (2) from the mean travel time (2). The coupling-ray-theory travel times at these receivers are again practically frequency independent between 10 Hz and 100 Hz.
Apart from the above mentioned singularities at the slowness surface, there are also less pronounced differences between the coupling-ray-theory travel times and the anisotropic-ray-theory travel times. As a rule, these "soft" differences approach Table 2. Receivers in velocity model SC1 I influenced by the split intersection singularity at the slowness surface. Differences D of the coupling-ray-theory travel times from the mean travel time in seconds at the frequency of 50 Hz are compared with the differences 1 2 |τ (2) −τ (1) | of the anisotropic-ray-theory travel times from the mean travel time in seconds. The anisotropic-ray-theory travel times are calculated by the second-order perturbation expansion along the reference rays traced in steps of 0.001 s. The unknown error of the perturbation expansion cannot influence the comparison because differences D are calculated using differences 1 2 |τ (2) −τ (1) |. The rounding errors should not influence the displayed differences. zero with increasing frequency ω. The importance of these soft differences can be evaluated in terms of expression The most important soft differences between the coupling-ray-theory travel times and the anisotropic-ray-theory travel times are observed in velocity model SC1 II, around the receiver at the depth of 0.96 km and around the frequency of 80 Hz, where the difference of the coupling-ray-theory travel time from the mean travel time is whereas the difference of the anisotropic-ray-theory travel time from the mean travel time is 1 2 τ (2) − τ (1) = 0.003767 s .
The rays leading to these receivers do not pass close to any slowness surface singularity (Pšenčík et al., 2012 ;Klimeš and Bulant, 2014b). The coupling-raytheory travel times in velocity model SC1 II approach the anisotropic-ray-theory travel times with increasing frequency. The next most important soft differences between the coupling-ray-theory travel times and the anisotropic-ray-theory travel times are observed in velocity model SC1 I, around the receiver at the depth of 0.04 km and around the frequency of 40 Hz, where the difference of the coupling-ray-theory travel time from the mean travel time is D = 0.007136 s , whereas the difference of the anisotropic-ray-theory travel time from the mean travel time is 1 2 τ (2) − τ (1) = 0.007488 s .
The coupling-ray-theory travel times at the shallow receivers in velocity model SC1 I approach the anisotropic-ray-theory travel times with increasing frequency. Within frequency band 10 Hz to 100 Hz, the differences between the coupling-raytheory travel times and the anisotropic-ray-theory travel times in velocity models QIH, QI, QI2 and QI4 are an order of magnitude less important than those in velocity model SC1 II, or at the shallow receivers in velocity model SC1 I. The coupling-ray-theory travel times approach the anisotropic-ray-theory travel times with increasing frequency.
The differences between the coupling-ray-theory travel times and the anisotropic-ray-theory travel times in velocity model ORT at the receivers outside the singular region, which extends from the depth of 0.72 km to 0.92 km, are even much smaller within the 10 Hz to 100 Hz frequency band. These coupling-ray-theory travel times approach the anisotropic-ray-theory travel times with increasing frequency.
The coupling-ray-theory travel times in velocity model KISS are equal to the anisotropic-ray-theory travel times up to rounding errors, independently of frequency within the 10 Hz to 100 Hz band.

CONCLUSIONS
Within a limited frequency band, we may efficiently approximate the frequencydependent frequency-domain coupling-ray-theory tensor Green function by two dyadic Green functions corresponding to two waves, described by their coupling-raytheory travel times and the corresponding vectorial amplitudes calculated for the prevailing frequency. The additional inaccuracy introduced by this prevailing-frequency approximation is usually smaller than the inaccuracy of the standard frequencydomain coupling ray theory, and smaller than the additional inaccuracy introduced by many other approximations of the coupling ray theory. The prevailing-frequency approximation of the coupling ray theory is accurately applicable in all situations we have tested numerically. The prevailing-frequency approximation is very accurate at both very weak anisotropy (Fig. 2) and strong anisotropy (Fig. 5), but it should further be tested with respect to stronger heterogeneities and broader frequency bands.
We have noticed no limitation of the applicability of the prevailing-frequency approximation of the coupling ray theory for electromagnetic waves or elastic S waves. On the other hand, before the prevailing-frequency approximation is applied in a particular medium, its accuracy should be tested there. If the results of the tests are satisfactory, we can apply the prevailing-frequency approximation.
The prevailing-frequency approximation of the coupling ray theory for electromagnetic waves or elastic S waves allows us to process the coupling-ray-theory wave field in the same way as the anisotropic-ray-theory wave field. This simplification may be decisive when storing the tensor Green function at the nodes of dense grids (Klimeš and Bulant, 2013 ), which is typical for applications such as the Born approximation. The prevailing-frequency approximation with its couplingray-theory travel times also offers a new way of understanding the results of the coupling ray theory. (12) In this section, we propose a convenient numerical representation of matrices Π, D, Π (1) or Π (2) . We also employ this representation in demonstrating that the 2×2 matrices (20) and (21) satisfy our assumption (12).