Interpolation of the coupling-ray-theory Green function within ray cells

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 pose a significant challenge in storing the Green function at the nodes of dense grids, typical for applications such as the Born approximation or non–linear source determination. Storing the Green function at the nodes of dense grids for too many frequencies may be impractical or even unrealistic. We have already proposed the approximation of the coupling–ray–theory tensor Green function, in the vicinity of a given prevailing frequency, by two coupling–ray–theory dyadic Green functions described by their coupling–ray–theory travel times and their coupling–ray–theory amplitudes. The above mentioned prevailing–frequency approximation of the coupling ray theory enables us to interpolate the coupling–ray–theory dyadic Green functions within ray cells, and to calculate them at the nodes of dense grids. For the interpolation within ray cells, we need to separate the pairs of prevailing–frequency coupling–ray–theory dyadic Green functions so that both the first Green function and the second Green function are continuous along rays and within ray cells. We describe the current progress in this field and outline the basic algorithms. The proposed method is equally applicable to both electromagnetic waves and elastic S waves. We demonstrate the preliminary numerical results using the coupling–ray–theory travel times of elastic S waves.


INTRODUCTION
In sufficiently smooth media, the anisotropic ray theory (Babich, 1956 ;Keller et al., 1956 ;Červený, 2001 ) can be used to calculate elastic P waves at all degrees of anisotropy, including isotropic media. However, neither the isotropic ray theory (Luneburg, 1944 ;Babich, 1961 ;Červený, 2001 ) nor the anisotropic ray theory is applicable to calculating electromagnetic waves or elastic S waves in many cases of heterogeneous anisotropic or bianisotropic media. We should use 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.
Unfortunately, the coupling-ray-theory tensor Green function is frequency-dependent, and is calculated separately for each frequency. This frequency dependence represents the main obstacle in interpolating the coupling-ray-theory tensor Green function within ray cells. Bulant (2012, 2016) found the approximation of the coupling-raytheory tensor Green function in the vicinity of a given prevailing frequency by two prevailing-frequency coupling-ray-theory dyadic Green functions corresponding to two elementary coupling-ray-theory waves. Each prevailing-frequency couplingray-theory dyadic Green function is described by its coupling-ray-theory travel time and its complex-valued coupling-ray-theory dyadic amplitude tensor, both calculated for the given prevailing frequency. Bulant (2012, 2016) numerically demonstrated that the prevailing-frequency coupling-ray-theory dyadic Green functions are well applicable in a reasonably broad frequency band around the given prevailing frequency. The prevailing-frequency coupling-ray-theory dyadic Green functions are calculated along the reference rays, which may be represented, e.g., by isotropic common reference rays, anisotropic common reference rays, anisotropic-ray-theory rays, or ordinary (SH) and extraordinary (SV) reference rays (Klimeš and Bulant, 2014, 2015Bulant and Klimeš, 2017 ). For the sake of conciseness, we shall refer hereinafter to the reference rays as the rays.
Each ray corresponds to two ray parameters and is represented by a point in the ray-parameter domain. The ray-parameter domain is triangulated according to Bulant (1996Bulant ( , 1999, see Fig. 9. The vertices of the triangles correspond to rays. Each triangle corresponds to a ray tube, and its three vertices correspond to three rays which form the edges of the ray tube. Each ray tube is sliced into ray cells according to Bulant and Klimeš (1999).
At each point of each ray, we have two prevailing-frequency coupling-ray-theory dyadic Green functions described by their coupling-ray-theory travel times and their coupling-ray-theory dyadic amplitude tensors. We wish to interpolate them within ray cells using the algorithm designed by Bulant and Klimeš (1999). To interpolate within ray cells, we need to separate the pairs of prevailing-frequency coupling-raytheory dyadic Green functions so that both the first Green function and the second Green function are continuous along rays and within ray cells. This paper is devoted to emerging challenges related to this separation.
In this paper, we implement the separation in two steps: In the first step, we copy each old ray to a pair of identical new rays and match the pair of prevailingfrequency coupling-ray-theory dyadic Green functions with the pair of new rays so that each Green function is continuous along the corresponding new ray. As a result, each of the three edges of each ray tube is represented by two new rays instead of one ray. In the second step, we double each ray tube and match the three pairs of new edge rays with the pair of ray tubes so that the Green function is continuous within both the ray tubes.
In this paper, we assume that both prevailing-frequency coupling-ray-theory tensor Green functions are calculated along a single anisotropic common reference ray (Klimeš, 2006 ;Klimeš and Bulant, 2006 ;Bulant and Klimeš, 2008 ). The anisotropic common reference rays then determine the ray cells for interpolation according to Bulant and Klimeš (1999).
The proposed method is equally applicable to both electromagnetic waves and elastic S waves. We demonstrate the preliminary numerical results using the coupling-ray-theory travel times of elastic S waves in Section 5.

PREVAILING-FREQUENCY APPROXIMATION OF THE COUPLING-RAY-THEORY GREEN FUNCTION
In the prevailing-frequency approximation, the coupling-ray-theory tensor Green function from point x 0 to point x is composed of two prevailing-frequency dyadic Green functions corresponding to two elementary coupling-ray-theory waves (Klimeš and Bulant, 2012, Eq. 39 ;2016, Eq. 23 ), (1) Each of the prevailing-frequency coupling-ray-theory dyadic Green functions on the right-hand side of relation (1) is described by its coupling-ray-theory travel time T (J) (x, x 0 , ω 0 ) and its complex-valued coupling-ray-theory dyadic amplitude tensor A (J) ij (x, x 0 , ω 0 ) (Klimeš and Bulant, 2012, Eq. 40 ;2016, Eq. 24 ), G both calculated for the given prevailing frequency ω 0 . The waves have couplingray-theory travel times (Klimeš and Bulant, 2012, Eqs 41-42 ;2016, Eqs 25-26 ) and Since we know the reference slowness vectors at the vertices of ray cells in addition to reference travel times τ , we interpolate reference travel time τ using tricubic interpolation within ray cells according to Bulant and Klimeš (1999). Since we do not know the spatial gradient of travel-time difference D, we separately interpolate travel-time difference D by trilinear interpolation within ray cells, and compose the coupling-ray-theory travel times according to the above definitions (3) and (4).
Before proceeding to the determination of the sign of travel-time difference D in Section 3, let us mention some properties of matrices Π and D useful for derivations in the next sections of this paper.
The evolution of matrix D with respect to the chosen monotonic parameter is given by ordinary differential equation which is obtained by differentiation of the coupling equation (9) with respect to ω. The algorithm for the numerical calculation of matrix D (Klimeš and Bulant, 2012, Sec. 4 ;2016, Sec. 4 ) is based directly on definition (11). Ordinary differential equation (13) is not used for numerical calculation of matrix D, but we shall use it in Section 3.
Unitary and unimodular complex-valued matrix can be expressed as (Klimeš and Bulant, 2012, Eq. 28 ; 2016, Eq. A-5 ) with real-valued coefficients Π 0 , Π 1 , Π 2 and Π 3 . Analogously, complex-valued matrix can be expressed as (Klimeš and Bulant, 2012, Eq. 31 ; 2016, Eq. A-8 ) with real-valued coefficients D 0 , D 1 , D 2 and D 3 . Its inverse matrix then reads Stud. Geophys. Geod., 61 (2017) 3. CONTINUITY OF THE COUPLING-RAY-THEORY GREEN FUNCTION ALONG RAYS Bulant (2012, 2016) proposed an algorithm for calculating the prevailing-frequency coupling-ray-theory dyadic Green functions at the receiver situated at the end of a two-point ray. The same algorithm can be used to calculate the prevailing-frequency coupling-ray-theory dyadic Green functions along any ray. The main problem that remains unsolved is the correct succession of the prevailingfrequency coupling-ray-theory dyadic Green functions. As we have seen above, this succession is determined by the sign of travel-time difference D. It is obvious from relations (7)-(8) that we should choose the sign of D so that matrix D/D is continuous along the ray.
For the sake of conciseness, we shall omit arguments x 0 and ω 0 of such functions as D(x, x 0 , ω 0 ) or D(x, x 0 , ω 0 ) in this section. We shall thus write just D(x) or D(x).
We calculate numerically all quantities at discrete points x k along the ray. Determinant used in definition (12) can be calculated as see (17). In order to control the change of sign from D(x k−1 ) to D(x k ), we define quantity Analogously to relations (19) and (20), quantity can be calculated as see (17) and (18). In order to make matrix D/D continuous along the ray, we put If the whole matrix D(x k ) vanishes, see Eq. (20). To calculate matrix D(x k )/D(x k ) used in relations (7) and (8), we apply l'Hospital's rule with respect to the derivative along the ray in this case, and substitute D ′ for D. Omitting multiplication factor ε ′ /ω, derivative (13) with matrix (14) yield substitution In the case of (25), we apply substitution (27) to calculate S(x k , x k−1 ), S(x k+1 , x k ) and D(x k )/D(x k ), but we naturally keep D(x k ) = 0 unchanged for calculating matrix D along the ray.

CONTINUITY OF THE COUPLING-RAY-THEORY GREEN FUNCTION WITHIN RAY TUBES
In place of each old ray, we now have two new rays corresponding to two prevailing-frequency coupling-ray-theory dyadic Green functions. The Green function is continuous along the corresponding new ray. As a result, each of the three edges of each old ray tube is represented by two new rays instead of one old ray. We need to double each old ray tube and match the three pairs of new edge rays with the pair of new ray tubes so that the Green function is continuous within both the new ray tubes.
At the points of each new ray, we calculate positive-definite Hermitian matrix Otherwise, if A 11 < A 22 , we define this vector as Complex-valued polarization vector e i at the point of the new ray then reads The real part of the polarization vector is Stud. Geophys. Geod., 61 (2017) and the imaginary part of the polarization vector is Each of the two edges of each side of an old ray tube is formed by two new rays corresponding to two prevailing-frequency coupling-ray-theory dyadic Green functions. For each side of an old ray tube, we thus have to compare the polarization vectors between four possible combinations of new edge rays. For each of these four combinations, we match the corresponding points of the two new edge rays of the combination, and compare polarization vectors e i andẽ i at each pair of the corresponding points. For each pair of the corresponding points, we calculate the squared cosine of angle Φ between the complex-valued polarization vectors e i andẽ i : We calculate minimum value [cos(Φ)] 2 min and maximum value [cos(Φ)] 2 max of (35) over all pairs of the corresponding points along the two rays. If the rays correspond to the same prevailing-frequency coupling-ray-theory dyadic Green function.
the rays correspond to different prevailing-frequency coupling-ray-theory dyadic Green functions. If both conditions (36) and (37) are violated, the rays cannot be compared. For example, if a conical singularity is situated close to the source within a sufficiently long old ray tube, for at least one side of the old ray tube we may expect [cos(Φ)] 2 max ≥ 3/4 with [cos(Φ)] 2 min ≪ 1/2 or [cos(Φ)] 2 min ≤ 1/4 with [cos(Φ)] 2 max ≫ 1/2 for all four possible pairs of corresponding new rays. In this case, we can create neither the corresponding side of a new ray tube, nor a new ray tube.
If two new edge rays of the side of an old ray tube correspond to the equal prevailing-frequency coupling-ray-theory dyadic Green function, and the other two new edge rays also correspond to the equal prevailing-frequency coupling-raytheory dyadic Green function, while the other two combinations of the new edge rays correspond to different prevailing-frequency coupling-ray-theory dyadic Green functions, the former two new edge rays form a possible side of a new ray tube, and the latter two new edge rays form a possible side of another new ray tube.
If the four new edge rays of each side of an old ray tube form two possible sides of two new ray tubes, and if three possible sides form a single new ray tube, the other three possible sides form another new ray tube. In this case, we can replace the old ray tube by the two new ray tubes.
The old ray tubes, which cannot be split into pairs of new ray tubes using the above described algorithm, should be studied further.

ELASTIC NUMERICAL EXAMPLES
The proposed method is equally applicable to both electromagnetic waves and elastic S waves. Here we numerically test the interpolation of the coupling-raytheory dyadic Green functions within ray cells using the elastic S waves. Following Pšenčík et al. (2012) and Bulant (2012, 2016), we consider eight anisotropic velocity models referred to as QIH, QI, QI2, QI4, KISS, SC1 I, SC1 II 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. For a sketch of the source-receiver configuration refer to Klimeš and Bulant (2016, Fig. 1).
We plot the relative coupling-ray-theory travel-time difference (half the relative coupling-ray-theory travel-time splitting) d = |D/τ | in the vertical rectangular section bounded by the point source from the left, and by the vertical well from the right (Klimeš and Bulant, 2012, Fig. 1 ; 2016, Fig. 1 ). The distance of the vertical well with receivers from the source is 1 km. The vertical extent of the rectangular section corresponds to the length of the vertical receiver profile considered by Klimeš and Bulant (2012, Figs 2-9;2016, Figs 2-9) for the calculation of the coupling-raytheory seismograms: QIH, QI, QI2, QI4 and KISS 0.6 km; SC1 I and SC1 II 1.4 km; ORT 1.6 km.
A vertically heterogeneous 1-D anisotropic velocity model QI was provided by Pšenčík and Dellinger (2001, model WA rotated by 45 • about the vertical axis). The same velocity model was used by Bulant and Klimeš (2002) and Klimeš and Bulant (2004) to demonstrate the coupling ray theory and its quasi-isotropic approximations. For a more detailed discussion and description of this velocity model refer to Pšenčík and Dellinger (2001). Velocity model QI is approximately uniaxial (approximately transversely isotropic). Its reference symmetry axis is horizontal and forms a 45 • angle with the vertical source-receiver plane (Klimeš, 2015(Klimeš, , 2016a. 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).
The relative coupling-ray-theory travel-time difference d = |D/τ | in the vertical source-receiver section in velocity models QIH, QI, QI2 and QI4 is displayed in Figs 1-4. Since the colour scale in each figure corresponds to the degree of anisotropy, the changes between Figs 1-4 illustrate the different development of S-wave coupling and splitting in dependence on anisotropy, which was discussed byČervený et al. (2007, Fig. 21). Figures 1-4 would be practically identical and close to Fig. 4 for anisotropic-ray-theory travel times.
For the description of velocity models KISS, SC1 I, SC1 II and ORT and the corresponding elastic moduli refer to Pšenčík et al. (2012).
Velocity model KISS represents velocity model QI described above, rotated by −44 • about the vertical axis in order to position the reference symmetry axis, corresponding to the kiss S-wave singularity, just 1 • from the source-receiver plane.  The relative coupling-ray-theory travel-time difference d = |D/τ | in the vertical source-receiver section is displayed in Fig. 5. Figures 2 and 5 thus represent two different vertical sections in the same velocity model.
Velocity model SC1 I is approximately uniaxial (approximately transversely isotropic) and its reference symmetry axis is horizontal (Klimeš, 2015(Klimeš, , 2016a. Its slowness surface contains a split intersection singularity, whereas velocity models QIH, QI, QI2 and QI4 display no exact or split intersection singularity (Pšenčík et al. 2012 ). Refer to Klimeš and Bulant (2016) for a more detailed discussion. The relative coupling-ray-theory travel-time difference d = |D/τ | in the vertical source-receiver section is displayed in Fig. 6.
Velocity model SC1 II is analogous to SC1 I, but its reference axis of symmetry is tilted. 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. The relative coupling-ray-theory travel-time difference d = |D/τ | in the vertical source-receiver section is displayed in Fig. 7.
In the 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 elastic S-wave polarization from the approximately anisotropic-ray-theory polarization S1 to the approximately anisotropic-ray-theory polarization S2, and vice versa. That  Stud. Geophys. Geod., 61 (2017) is probably why the ray tubes in the vicinity of the conical singularity in velocity model ORT cannot be split into the pairs of tubes with continuous prevailingfrequency coupling-ray-theory dyadic Green functions. The relative coupling-raytheory travel-time difference d = |D/τ | in the vertical source-receiver section is displayed in Fig. 8. The diagonal grey zone corresponds to the ray tubes which cannot be split into the pairs of tubes with continuous prevailing-frequency couplingray-theory dyadic Green functions according to Section 4. The ray tubes in the rayparameter domain in velocity model ORT are displayed in Fig. 9. We can observe the regions around three of the four conical singularities where the ray tubes cannot be split into the pairs of tubes with continuous prevailing-frequency coupling-raytheory dyadic Green functions according to Section 4. The bottom-right region corresponds to the diagonal grey zone in Fig. 8.

CONCLUSIONS
The prevailing-frequency approximation of the coupling ray theory allows us to process the coupling-ray-theory wave field in terms of travel times and amplitudes, i.e. in the same way as the anisotropic-ray-theory wave field. The prevailingfrequency approximation of the coupling ray theory can thus be included in wavefront tracing and in the interpolation within ray cells in anisotropic or bianisotropic media, provided that we can separate the pairs of prevailing-frequency coupling-ray-theory dyadic Green functions so that both the first Green function and the second Green function are continuous within the ray cells. In this paper, we have proposed the preliminary version of the separation algorithm and tested it using eight numerical examples.
The interpolation of the prevailing-frequency coupling-ray-theory dyadic Green functions within ray cells is essential for the ray-based Born approximation, nonlinear source determination and other applications in anisotropic or bianisotropic media.
In the proposed preliminary version of the separation algorithm, the interpolation is possible only in ray tubes in which the prevailing-frequency coupling-ray-theory dyadic Green functions are continuous. This requirement may result in many whole ray tubes where we cannot interpolate, as in Fig. 8 in velocity model ORT. In order to refine the algorithm and extend the region where we can interpolate, we should study the continuity of the prevailing-frequency coupling-ray-theory dyadic Green functions in the individual ray cells rather than in whole ray tubes.
First Normalized Ray Parameter Second Normalized Ray Parameter Fig. 9. Ray tubes in velocity model ORT for the interpolation of coupling-ray-theory travel times. The rays are represented by points (here crosses) and the ray tubes by triangles in the ray-parameter domain. The three regions with missing triangles (ray tubes) are situated around the conical singularities where the ray tubes cannot be split into the pairs of tubes with continuous prevailing-frequency coupling-ray-theory dyadic Green functions according to Section 4. The bottom-right region corresponds to the diagonal grey zone in Fig. 8. Stud. Geophys. Geod., 61 (2017)