Comparison of geometric optics and «aperture» methods for calculation of the electromagnetic radiation caused by charged particles flying by dielectric objects

Full Text

Introduction

A problem in the physics of wave processes is determining the electromagnetic field excited by a charged particle or a collection of particles flying past or through a certain dielectric object. As a rule, exact solutions to such problems cannot be constructed because they are cumbersome, poorly interpreted, and require labor-intensive calculations. Simultaneously, numerical calculations require large computer resources. For this reason, most problems are solved using an approximation.

This article discusses in detail the application of the stationary phase method to aperture integrals, which provide a solution in the case of large (on the wavelength-scale) objects, including those in the geometric–optical region. This approach has some difficulties; namely, systems of equations for finding the saddle point can often be solved only numerically, the second derivatives of the integrands for calculating the determinant turn out to be too cumbersome, etc. In contrast, when the saddle point is found for a single integral, the stationary phase method greatly simplifies the problems [1].

An alternative to this approach is geometric–optical calculation [2] based on the method for reflected rays described by V.A. Fok [3]. The essence of this method can be described in several steps. First, the field excited inside an object without external boundaries is calculated. Second, the field at the external boundary of the object is determined using Fresnel coefficients. Finally, the path of the rays and the expansion of the ray tubes in the area outside the object are considered.

Inverse problems exist, where the properties of the medium are determined based on the change in the field after reflection or transmission [4–6].

This work compares two approaches for determining the field in the geometric–optical region using the example of a problem with a dielectric prism, as described in Section 1, and calculates caustics in Section 2 for some of the most common dielectric objects based on the use of equations for ray tube expansion. All studies were performed for the Fourier transform of the electric field (the calculation of the field itself, which reduces to the calculation of the corresponding inverse Fourier integral, is beyond the scope of this work).

1. Advantages of the geometric–optical method

Let us consider the advantages of the geometric–optical method using the following problem as an example. A particle with charge q and speed $\overrightarrow{V}=\beta \overrightarrow{c}$ (where c is the speed of light) at a distance a from its lower boundary travels along the boundary of a dielectric prism with dielectric permeability ε and magnetic permeability $\mu =1$. The size of the prism is considered large compared to the wavelength under consideration. The prism is located in a medium in which differences from vacuum are unsubstantial (Fig. 1). We assume that radiation is emitted only from a rectangular region centered at a point $\left(x_c, y_c, z_c\right)$ with sides d and b parallel to the axes $\xi$ and $\eta$ respectively; accordingly, the remainder of the prism is closed to radiation.

 

Fig. 1. Central cross-section of the prism with main notations

Рис. 1. Центральное сечение призмы с основными обозначениями

 

A geometric–optical solution is known for this problem [7]. Let us further consider the aspects of solving this problem using the “aperture” method (i.e., the method based on the Stratton–Chu equations). This solution is much more general because it enables the determination of the field in areas where geometric optics is invalid. However, this method also has certain disadvantages if the task consists of analyzing the field only in the region where geometric optics is valid.

The field components obtained using the aperture methods are seen in [8]. They have the following form:

$\left\{\begin{array}{l}{E}_{\xi }^{\left(h\right)}\left(R\right)\\ E_{\eta }^{\left(h\right)}\left(R\right)\\ E_{\zeta }^{\left(h\right)}\left(R\right)\end{array}\right\}=\frac{ik}{4\pi }\underset{-d/2}{\overset{d/2}{\int }}d\xi \text{'}\underset{-b/2}{\overset{b/2}{\int }}d\eta \text{'}\frac{G\left(\left|\Delta \overrightarrow{R}\text{'}\right|\right)}{{\left(\left|\Delta \overrightarrow{R}\text{'}\right|\right)}^2}\times$

$\times \left\{\begin{array}{l}-\left[{\left(\Delta \eta \text{'}\right)}^2+{\zeta }^2\right]H_{\eta }\left(\overrightarrow{R}\text{'}\right)-\Delta \xi \text{'}\Delta \eta \text{'}{H}_{\xi }\left(\overrightarrow{R}\text{'}\right)\\ \left[{\left(\Delta \xi \text{'}\right)}^2+{\zeta }^2\right]H_{\xi }\left(\overrightarrow{R}\text{'}\right)+\Delta \xi \text{'}\Delta \eta \text{'}{H}_{\eta }\left(\overrightarrow{R}\text{'}\right)\\ -\zeta \Delta \eta \text{'}{H}_{\xi }\left(\overrightarrow{R}\text{'}\right)+\zeta \Delta \xi \text{'}{H}_{\eta }\left(\overrightarrow{R}\text{'}\right)\end{array}\right\},$ (1)

$\left\{\begin{array}{c}{E}_{\xi }^{\left(e\right)}\left(R\right)\\ E_{\eta }^{\left(e\right)}\left(R\right)\\ E_{\zeta }^{\left(e\right)}\left(R\right)\end{array}\right\}=\frac{ik}{4\pi }\underset{-d/2}{\overset{d/2}{\int }}d\xi \text{'}\underset{-b/2}{\overset{b/2}{\int }}d\eta \text{'}\frac{G\left(\left|\Delta \overrightarrow{R}\text{'}\right|\right)}{\left|\Delta \overrightarrow{R}\text{'}\right|}\times$

$\times \left\{\begin{array}{c}-\zeta E_{\xi }\left(\overrightarrow{R}\text{'}\right)\\ -\zeta E_{\eta }\left(\overrightarrow{R}\text{'}\right)\\ \Delta \xi \text{'}{E}_{\xi }\left(\overrightarrow{R}\text{'}\right)+\Delta \eta \text{'}{E}_{\eta }\left(\overrightarrow{R}\text{'}\right)\end{array}\right\}.$ (2)

The dashed line marks the coordinates of the prism surface (integration variables). For clarity, the electric field strength is divided into two components, as indicated by the superscripts (e) and (h), and the integrand functions depend on the electric and magnetic field strength on the prism surface. The projections $\overrightarrow{E}$ and $\overrightarrow{H}$ of the electromagnetic field in the integrands denote the field on the upper outer surface of the prism (they are given in [8]). In Eqs. (1) and (2), the following notations are used: $k=\omega /c$ is the wave number; i is the imaginary unit, $\Delta \xi \text{'}=\xi -\xi \text{'}$, $\Delta \eta \text{'}=\eta -\eta \text{'}$, $\overrightarrow{R}\text{'}=(\xi \text{'},\eta \text{'},0)$, and $G\left(∆\overrightarrow{R}\text{'}\right)$ is Green’s function:

$G\left(\left|\Delta \overrightarrow{R}\text{'}\right|\right)=\frac{e^{ik\left|\Delta \overrightarrow{R}\text{'}\right|}}{\left|\Delta \overrightarrow{R}\text{'}\right|}=\frac{e^{ik\sqrt{{\left(\xi -\xi \text{'}\right)}^2+{\left(\eta -\eta \text{'}\right)}^2+{\zeta }^2}}}{\sqrt{{\left(\xi -\xi \text{'}\right)}^2+{\left(\eta -\eta \text{'}\right)}^2+{\zeta }^2}}.$ (3)

To calculate the integrals of Eqs. (1) and (2), the stationary phase method can be used. N-fold integrals with a large positive parameter of the form Ω

${\widehat{I}}_n\left(\Omega \right)=\int \underset{-\infty }{\overset{+\infty }{}}\int f\left(x\right)e^{i\Omega q\left(x\right)}d{x}_1\mathrm{...}d{x}_n,\Omega >0$ (4)

are equal to [9; 10]

${\widehat{I}}_n\left(\Omega \right)={\left(\frac{2\pi }{\Omega }\right)}^{n/2}\times$

$\times \exp \left[i\Omega q\left(x_s\right)+\frac{i\pi }{4}\sigma \right]\frac{f\left(x_s\right)+O\left({\Omega }^{-1}\right)}{{\left|\det \frac{{\partial }^{2}q}{\partial x_{is}\partial x_{js}}\right|}^{1/2}}.$ (5)

Here,

$\sigma =\sum _{i=1}^n\mathrm{sgn}\left(d_i\right),$

where d_i is the eigenvalue of the matrix

$\frac{{\partial }^{2}q}{\partial x_{is}\partial x_{js}}$, $i,j=\mathrm{1,...,}n$,

and x_js are the coordinates of the stationary point. To determine these coordinates, we must determine the points of simultaneous zeroing of the derivatives of the phase of the integrand with respect to $\xi \text{'}$, $\eta \text{'}$. 

In Eqs. (1) and (2), considering the field on the prism surface, equations for which are given in [8], a phase is identified that is identical for all components, which are sixfold:

${\Phi }_1=ik\sqrt{{\left(\xi -\xi \text{'}\right)}^2+{\left(\eta -\eta \text{'}\right)}^2+{\zeta }^2}$, (6)

${\Phi }_2=-a\sqrt{{\kappa }^2+k_y^2},$

${\Phi }_3=i\left(x-a\right)g_m,$

${\Phi }_4=i{k}_{y}y,$

${\Phi }_5=ik{\beta }^{-1}z,$

${\Phi }_6=-i\pi /4$

where

$\kappa =\frac{k}{\beta }\sqrt{1-{\beta }^2}$, $k_y=\frac{k\sqrt{\epsilon {\beta }^2-1}\eta \text{'}}{\beta \sqrt{{\left(x_c+\xi \text{'}\sin \alpha -a\right)}^2+\eta {\text{'}}^2}}$, $g_m=\sqrt{\frac{k^2\left(\epsilon {\beta }^2-1\right)}{{\beta }^2}-k_y^2}$, $x=x_c+\xi \text{'}\sin \alpha$, $y=\eta \text{'}$, $z=z_c-\xi \text{'}\cos \alpha$.

The results of differentiating the imaginary part of the phase $\Phi$ are presented below:

$\frac{\partial \Phi }{\partial \xi \text{'}}=g_m\sin \alpha \left[1+\frac{k_y^2}{g_m^2}\frac{{\left(x_c+\xi \text{'}\sin \alpha -a\right)}^2}{{\left(x_c+\xi \text{'}\sin \alpha -a\right)}^2+\eta {\text{'}}^2}\right]-$

$-\frac{k_y\eta \text{'}\left(x_c+\xi \text{'}\sin \alpha -a\right)\sin \alpha }{{\left(x_c+\xi \text{'}\sin \alpha -a\right)}^2+\eta {\text{'}}^2}-\frac{k\cos \alpha }{\beta }-$

$-\frac{k\left(\xi -\xi \text{'}\right)}{\sqrt{{\left(\xi -\xi \text{'}\right)}^2+{\left(\eta -\eta \text{'}\right)}^2+{\zeta }^2}},$ (7)

$\frac{\partial \Phi }{\partial \eta \text{'}}=k_y\left[1+\frac{{\left(x_c+\xi \text{'}\sin \alpha -a\right)}^2}{{\left(x_c+\xi \text{'}\sin \alpha -a\right)}^2+\eta {\text{'}}^2}\right]-$

$-\frac{k_y^2}{g_m\eta \text{'}}\frac{{\left(x_c+\xi \text{'}\sin \alpha -a\right)}^3}{{\left(x_c+\xi \text{'}\sin \alpha -a\right)}^2+\eta {\text{'}}^2}-$

$-\frac{k\left(\eta -\eta \text{'}\right)}{\sqrt{{\left(\xi -\xi \text{'}\right)}^2+{\left(\eta -\eta \text{'}\right)}^2+{\zeta }^2}}.$

The stationary point is the solution of the following system of equations with respect to $\xi \text{'}, \eta \text{'}$:

$\frac{\partial \Phi }{\partial \xi \text{'}}=f_1=\mathrm{0,}$

$\frac{\partial \Phi }{\partial \eta \text{'}}=f_2=0.$ (8)

This nonlinear system of equations can be solved numerically. For the problem under consideration, Newton’s method was chosen for the greatest convenience.

Newton’s method involves searching for a solution over a certain number of iterations until the absolute difference between the new and previous values is less than a predetermined threshold.

To solve the system $\left(\begin{array}{c}{f}_1\\ f_2\end{array}\right)=0$, we use the following recursion:

$\frac{\partial f_1}{\partial \xi \text{'}}\left({\xi }_n^{\text{'}},{\eta }_n^{\text{'}}\right)\delta {\xi }_n^{\text{'}}+\frac{\partial f_1}{\partial \eta \text{'}}\left({\xi }_n^{\text{'}},{\eta }_n^{\text{'}}\right)\delta {\eta }_n^{\text{'}}=-f_1\left({\xi }_n^{\text{'}},{\eta }_n^{\text{'}}\right),$

$\frac{\partial f_2}{\partial \xi \text{'}}\left({\xi }_n^{\text{'}},{\eta }_n^{\text{'}}\right)\delta {\xi }_n^{\text{'}}+\frac{\partial f_2}{\partial \eta \text{'}}\left({\xi }_n^{\text{'}},{\eta }_n^{\text{'}}\right)\delta {\eta }_n^{\text{'}}=-f_2\left({\xi }_n^{\text{'}},{\eta }_n^{\text{'}}\right),$ (9)

where $\delta {\xi }_n^{\text{'}}={\xi }_{n+1}^{\text{'}}-{\xi }_n^{\text{'}}$ and $\delta {\eta }_n^{\text{'}}={\eta }_{n+1}^{\text{'}}-{\eta }_n^{\text{'}}$ are changes in the values of the required points at the n-th step, respectively. These changes are easy to calculate using Cramer’s method, and the result can then be compared with a predetermined value. To stop recurrence, both increments in absolute value must be less than a predetermined threshold. It is also necessary to set the starting point $({\xi }_0^{\text{'}},{\eta }_0^{\text{'}}),$ for example, at the origin with a slight offset.

The results of differentiating $f_1, f_2$ are given below:

$\frac{\partial f_1}{\partial \xi \text{'}}=\frac{\partial g_m}{\partial \xi \text{'}}\sin \alpha \left[1-\frac{{\Xi }^2}{{\Sigma }^2}\right]+\frac{\partial k_y}{\partial \xi \text{'}}\sin \alpha \frac{X}{{\Sigma }^2}\left[2\Xi -\eta \text{'}\right]+$

$+\frac{k_y{\sin }^2\alpha }{{\Sigma }^4}\left[2\Xi \eta {\text{'}}^2-\eta \text{'}\left(\eta {\text{'}}^2-X^2\right)\right]+$

$+\frac{k{\left(\eta -\eta \text{'}\right)}^2+k{\zeta }^2}{{\left|\Delta \overrightarrow{R}\text{'}\right|}^3},$ (10)

$\frac{\partial f_1}{\partial \eta \text{'}}=\frac{\partial g_m}{\partial \eta \text{'}}\sin \alpha \left[1-\frac{{\Xi }^2}{{\Sigma }^2}\right]+\frac{\partial k_y}{\partial \eta \text{'}}\sin \alpha \frac{X}{{\Sigma }^2}\left[2\Xi -\eta \text{'}\right]+$

$+\frac{k_y\sin \alpha X}{{\Sigma }^2}\left[\frac{2\eta {\text{'}}^2-2\Xi \eta \text{'}}{{\Sigma }^2}-1\right]-\frac{k\left(\xi -\xi \text{'}\right)\left(\eta -\eta \text{'}\right)}{{\left|\Delta \overrightarrow{R}\text{'}\right|}^3},$

$\frac{\partial f_2}{\partial \eta \text{'}}=\frac{\partial g_m}{\partial \eta \text{'}}\frac{1}{\eta \text{'}}\frac{{\Xi }^{2}X}{{\Sigma }^2}+\frac{\partial k_y}{\partial \eta \text{'}}\left[1+\frac{X^2}{{\Sigma }^2}-\frac{2\Xi X^2}{\eta \text{'}{\Sigma }^2}\right]+$

$+\frac{k_y{X}^2}{{\Sigma }^2}\left[\frac{\Xi }{\eta {\text{'}}^2}-\frac{2\eta \text{'}}{{\Sigma }^2}+\frac{2\Xi }{{\Sigma }^2}\right]+\frac{k{\left(\xi -\xi \text{'}\right)}^2+k{\zeta }^2}{{\left|\Delta \overrightarrow{R}\text{'}\right|}^3},$

$\frac{\partial f_2}{\partial \xi \text{'}}=\frac{\partial f_1}{\partial \eta \text{'}}.$

Here,

$X=x_c+\xi \text{'}\sin \alpha -a,$

$\Sigma =\sqrt{{\left(x_c+\xi \text{'}\sin \alpha -a\right)}^2+\eta {\text{'}}^2},$ and

$\Xi =\frac{k_y}{g_m}\left(x_c+\xi \text{'}\sin \alpha -a\right).$

All results obtained were programed and compared with those described in [7]. However, first, the stationary point within the framework of this problem cannot be calculated analytically, which reduces the accuracy of the calculations. Second, when searching for stationary points, “blind zones” for calculations appear, i.e., observation points for which the point of the stationary phase cannot be determined using the described methods. For example, in the plane ${\zeta }_{obs}=\text{const},$ where ${\zeta }_{obs}$ is the distance from the top surface of the prism (Fig. 1), areas (with a peak field value) can be seen where the stationary point is not located. The further the plane is from the prism surface, the larger the area. Such a blind zone is presented in Fig. 2.

 

Fig. 2. Example of «blind spot» for stationary phase points calculations. Areas of undefined stationary points are shown in black. All the distances are in $c/\omega$. Parameters for calculations: $\epsilon =2, a=1, \alpha =\pi /6, \beta =0,9, \zeta =30, b=d=50$

Рис. 2. Пример «слепой зоны» точек для вычисления стационарной фазы. Черным обозначены области, в которых не вычисляется точка стационарной фазы. Для наглядности расстояния измеряются в единицах $c/\omega$. При этом использованы следующие параметры: $\epsilon =2, a=1, \alpha =\pi /6, \beta =0,9, \zeta =30, b=d=50$

 

Let us compare the electric field obtained using the stationary phase method and by considering the ray tube, for lines $\zeta =\text{const}$ and $\xi =\text{const}$. The distribution pattern of the field amplitude in space is identical, and the field magnitude is of the same order. However, with increasing distance from the prism within the framework of the geometric optics approximation and increasing the dielectric permeability at peak values, the discrepancy in field strength can be up to 15 %.

Both approaches correctly reflect the structure and direction of the field, but the stationary phase method provides a smoother solution.

Figure 3 shows a comparison of the field at various charged particle velocities and distances from the prism surface.

 

Fig. 3. Comparison of two methods for Fourier transform electric field magnitude measured in $V·s/m$ for prism where: $q=1$ nC, $\epsilon =4, a=c/\omega , \alpha =\pi /6, b=d=50a$

Рис. 3. Сравнение результатов, полученных методом стационарной фазы и методом геометрической оптики, для Фурье-образа электрического поля в $V·s/m$ для призмы с параметрами: $q=1$ нКл, $\epsilon =4, a=c/\omega , \alpha =\pi /6, b=d=50a$

 

2. Calculation of caustics

According to [11], caustics can be calculated (or their absence can be shown) using equations for the cross section of the ray tube D(l). For this purpose, the equation $D\left(l\right)=0$ must be solved together with the system of ray path equations. Let us consider three dielectric objects: a cone, a prism, and a sphere (Fig. 4). The source of radiation, as before, is a moving point charge. Geometric–optical solutions to these problems, as well as equations for the cross section of the beam tube, are given in [2; 7; 12], respectively.

 

Fig. 4. Considered dielectic objects in central plane cross section with main notations

Рис. 4. Центральное сечение рассматриваемых диэлектрических объектов

 

Calculating caustics for a cone with a channel (Fig. 4, left) is not difficult. Because $D\left(l\right)=\frac{\rho }{\rho \text{'}}\cos {\theta }_t$ [2], the solution is as follows:

${\theta }_t=\frac{\pi }{2}+\pi n,n\in \mathbb{Z},$

$\rho =0$. (11)

Thus, caustics can form on the surface of the cone at the limiting angle of total internal reflection or on the axis of the cone. Both situations are of no practical interest.

For a prism (Fig. 4, center), the situation is also quite simple. The equation for the cross section of the ray tube is as follows [7]:

$D\left(l\right)=\cos {\theta }_t+$

$+\frac{l}{\rho \text{'}}\sqrt{\epsilon }\sin {\theta }_p\left[\cos {\theta }_t+\frac{{\cos }^2\alpha {\sin }^2\phi \text{'}}{\cos {\theta }_t}\left(\epsilon -1\right)\right].$ (12)

The equation $D\left(l\right)=0$ can be reduced to the form

${\cos }^2{\theta }_t\left(1+\frac{l}{\rho \text{'}}\sqrt{\epsilon }\sin {\theta }_p\right)=$

$=-\frac{l}{\rho \text{'}}\sqrt{\epsilon }{\cos }^2\alpha {\sin }^2\phi \text{'}\sin {\theta }_p\left(\epsilon -1\right).$ (13)

During the generation of Cherenkov radiation, $\sin {\theta }_p$ takes real positive values, so Eq. (13) cannot be satisfied that indicates the absence of caustics.

For a sphere (Fig. 4, right), the equation for the ray tube cross section has the following form [12]:

$D\left(l\right)=\cos {\theta }_t+$

$+\frac{l}{R\text{'}}\left[\frac{\sin \left({\theta }_t-{\theta }_i\right)}{\sin {\theta }_i\cos {\theta }_t}-\frac{\sin \left(\theta \text{'}-{\theta }_t\right)\cos {\theta }_t}{\sin \theta \text{'}}\right]+$

$+{\left(\frac{l}{R\text{'}}\right)}^2\frac{\sin \left({\theta }_t-\theta \text{'}\right)\sin \left({\theta }_t-{\theta }_i\right)}{\cos {\theta }_t\sin {\theta }_i\sin \theta \text{'}}.$ (14)

The quadratic equation $D\left(l\right)=0$ has two solutions:

$\frac{l}{R\text{'}}=\frac{\sin \theta \text{'}}{\sin \left({\theta }_t-\theta \text{'}\right)},$

$\frac{l}{R\text{'}}=\frac{{\cos }^2{\theta }_t\sin {\theta }_i}{\sin \left({\theta }_t-{\theta }_i\right)},$ (15)

where $\left(R, \theta , \phi \right)$ and $\left(R\text{'}, \theta \text{'}, \phi \text{'}\right)$ are the spherical coordinates at the observation point and the point on the sphere surface, respectively.

To determine the relationship between the coordinates of the caustic points, we use the following equation for the ray length:

$l=\sqrt{R^2+R{\text{'}}^2-2RR\text{'}\cos \left(\theta -\theta \text{'}\right)}.$ (16)

Because of the cylindrical symmetry of the problem, it is sufficient to consider a section through the channel axis, and without loss of generality, we can proceed from the surfaces of the caustics to the lines they formed in the section under consideration. Hereinafter, such lines will be called caustic lines. Let us first determine whether such lines are located on the rays. For this purpose, we use the initial condition that at the point of the ray exiting the sphere, we have $l=0, R=R\text{'}$, and $\theta =\theta \text{'}$. In this case, we obtain two solutions of Eq. (15), 1) a nonphysical solution $\theta \text{'}=\pi n, n=0, \pm 1, \pm 2, ...$ and 2) a solution suggested by the geometry of the problem, namely, the condition $\theta \text{'}={\theta }_i+{\theta }_p$ where the sign of ${\theta }_i$ is determined by the position relative to the normal (positive on the left and negative on the right if viewed on the sphere in Fig. 4). Considering Snell’s law $\sin {\theta }_t=\sqrt{\epsilon }\sin {\theta }_i$, we obtain the following final equation for the ray exit points from the sphere:

$\theta \text{'}={\theta }_p\pm \arcsin \frac{1}{\sqrt{\epsilon }}+\pi n,n=\mathrm{0,}\pm \mathrm{1,}\pm \mathrm{2,}\mathrm{...}$ (17)

The rays emitted from the resulting points cut off the areas with the intersections of the rays from the side of the sphere, as shown in Fig. 5. Now, we need to determine the curved lines of the caustics, which are clearly visible in Fig. 5 with simulated rays.

 

Fig. 5. Dielectric ball caustics for case of $\epsilon =4, \beta =0,999$. Thin (red in online version of the article) lines are rays and thick (blue in online version of the article) lines show caustics

Рис. 5. Каустики шара. Для примера рассмотрен случай с $\epsilon =4, \beta =0,999$. Тонкими (красными в онлайн-версии статьи) линиями изображены лучи, жирными (синими в онлайн-версии статьи) – вычисленные линии каустик

 

If the caustics are not located on the rays, the initial conditions cannot be applied to determine them, and a different approach is required. To begin with, we connect $\theta , \theta \text{'}$ and $R$ from the three Eqs. (15)–(16) above, containing l:

$R=R\text{'}\left[\cos \left(\theta -\theta \text{'}\right)\pm \sqrt{\frac{{\cos }^4{\theta }_t{\sin }^2{\theta }_i}{{\sin }^2\left({\theta }_t-{\theta }_i\right)}-{\sin }^2\left(\theta -\theta \text{'}\right)}\right],$ (18)

$R=R\text{'}\left[\cos \left(\theta -\theta \text{'}\right)\pm \sqrt{\frac{{\sin }^2\theta \text{'}}{{\sin }^2\left({\theta }_t-{\theta }_i\right)}-{\sin }^2\left(\theta -\theta \text{'}\right)}\right].$

However, to construct caustic lines, a connection of R and $\theta$ is required, and the angles ${\theta }_t, {\theta }_i$ are expressed through $\theta \text{'}$ and the constants of the problem. Therefore, we will use the ray equation as an additional equation:

$\frac{x-x\text{'}}{{\kappa }_x^{*}}=\frac{z-z\text{'}}{{\kappa }_z^{*}},$ (19)

where [12]:

${\overrightarrow{\kappa }}^{*}=\frac{\sin {\theta }_t}{\sin {\theta }_i}\overrightarrow{\kappa }-\frac{\sin \left({\theta }_t-{\theta }_i\right)}{\sin {\theta }_i}\overrightarrow{N},$ (20)

$\overrightarrow{N}=\left(\begin{array}{l}{N}_x\\ N_y\\ N_z\end{array}\right)=\left(\begin{array}{l}\cos \phi \text{'}\sin \theta \text{'}\\ \sin \phi \text{'}\sin \theta \text{'}\\ \cos \theta \text{'}\end{array}\right),$ (21)

$\overrightarrow{\kappa }=\left(\begin{array}{l}{\kappa }_x\\ {\kappa }_y\\ {\kappa }_z\end{array}\right)=\left(\begin{array}{l}\cos \phi \text{'}\sin {\theta }_p\\ \sin \phi \text{'}\sin {\theta }_p\\ \cos {\theta }_p\end{array}\right).$ (22)

Here, $\overrightarrow{\kappa }, {\overrightarrow{\kappa }}^{*}$ and $\overrightarrow{N}$ are the vectors directed along the incident and refracted rays and the normal restored at the point of incidence, respectively.

Because of the transformations and considering the cylindrical symmetry of the problem, we obtain the following equation:

$R=\frac{R\text{'}\sin {\theta }_t}{\sin \left({\theta }_t+\theta -\theta \text{'}\right)}.$ (23)

The equation for connecting points is a complex transcendental equation, so the problem of determining points on the caustic line was solved numerically when only the points that obeyed Eqs. (18) and (23) were selected simultaneously. In Fig. 5, thick blue lines indicate caustics (located on the rays and curves), and thin red lines indicate the rays emitted from the sphere, for clarity.

Conclusion

The use of geometric optics, which involves determining the path of rays and the expansion of beam tubes, allows for improved visualization of the field for complex cases compared with the stationary phase method. In addition, after obtaining the cross section of the beam tube, the wave field can be further examined for the presence of caustics, the analysis of which enables us to determine the weak points of the geometric–optical approach. The use of geometric optics is often less labor-intensive than calculating aperture integrals using the stationary phase method.

This work was partially supported by the Russian Science Foundation grant No. 18-72-10137.

The author expresses gratitude to A.V. Tyukhtin for useful discussions.

About the authors

Ekaterina S. Belonogaya

Saint Petersburg University

Email: ekaterinabelonogaya@yandex.ru
ORCID iD: 0000-0003-0557-1710

master of physics, research engineer of the Department of Radiophysics of Faculty of Physics, Saint Petersburg University, Saint Petersburg, Russia.

Research interests: Cherenkov radiation, geometric optics of heterogeneous media, wave processes.

Russian Federation, 7-9, Universitetskaya Embankment, Saint Petersburg, 199034, Russia

Dmitriy S. Klyuev

Povolzhskiy State University of Telecommunications and Informatics

Author for correspondence.
Email: klyuevd@yandex.ru
ORCID iD: 0000-0002-9125-7076

Doctor of Physical and Mathematical Sciences, head of the Department of Radioelectronic Systems, Povolzhskiy State University of Telecommunications and Informatics, Samara, Russia. Author of over 250 scientific papers.

Research interests: electrodynamics, microwave devices, antennas, metamaterials.

Russian Federation, 23, L. Tolstoy Street, Samara, 443010, Russia

References

Supplementary files