On the Electromagnetic Propagation Paths and Wavefronts in the Vicinity of a Refractive Object

Maxwell's equations are transformed from a Cartesian geometry to a Riemannian geometry. A geodetic path in the Riemannian geometry is defined as the raypath on which the electromagnetic energy efficiently travels through the medium. Consistent with the spatial behavior of the Poynting vector, the metric tensor is required to be functionally dependent on the refractive index of the medium. A symmetric nonorthogonal transformation is introduced, in which the metric is a function of an electromagnetic tension. This so‐called refractional tension determines the curvature of the geodetic line. To verify the geodetic propagation paths and wavefronts, a spherical object with a refractive index not equal to one is considered. A full 3‐D numerical simulation based on a contrast‐source integral equation for the electric field vector is used. These experiments corroborate that the geodesics support the actual wavefronts. This result has consequences for the explanation of the light bending around the Sun. Next to Einstein's gravitational tension there is room for an additional refractional tension. In fact, the total potential interaction energy controls the bending of the light. It is shown that this extended model is in excellent agreement with historical electromagnetic deflection measurements.


Introduction
The transmission of electromagnetic energy along rays has been of interest in our community. The question has always been, see, for example, Cheney (2004), how accurate the ray type of approximation is in representing the actual state of affairs of sending and receiving electromagnetic signals. In fact it is a high-frequency approximation derived from Maxwell's equations, leading to the eikonal equation, see p. 111 of Born and Wolf (1959). It is assumed that the signal is well quantified by the raypath description honoring the travel times of the signal, while its amplitude is stationary. This would then be a sufficient basis to, for example, image an object. In this paper, we argue that the ray approximation does not meet up to its expectations. It fails to image the boundary of an object correctly, due to the fact that the method does not support all frequencies.
In our analysis we start with Maxwell's equations in tensor format going from a Cartesian to a Riemannian geometry with a given metric tensor which is fundamental for that geometry. We argue that the geodetic path in the Riemannian geometry determines the raypath on which the electromagnetic energy efficiently travels through the medium. From the behavior of the stationarity of the Poynting vector we conclude that the metric tensor is functionally dependent on the refractive index of the medium. We introduce a symmetric nonorthogonal transformation. The metric is a function of an electromagnetic potential, which gives rise to a tension that determines the curvature of the geodetic line. This tension is also present in vacuum outside the object. To illustrate the propagating wavefronts, we consider a ray of electromagnetic energy passing a spherical object with a refractive index not equal to 1. We carry out a full 3-D numerical simulation based on the contrast-source integral equation for the electric field vector. These experiments corroborate our wave-front analysis that the geodesics support the wavefront of the numerical simulation. This result has consequences for the explanation of light bending around the Sun. Next to the Einstein's gravitational tension, there is room for an additional refractional tension. In fact, the total potential interaction energy controls the bending of the light. We conclude this paper by showing that this model is in excellent agreement with historical electromagnetic deflection measurements.

10.1029/2019RS007021
was caused by the "heaviness" of the light in reaction to the gravitational force of the Sun. Einstein (1911) conclusion was that the gravitational force is nothing else but a curvature of space. We are not debating this conclusion, but we see, based on Maxwell's equation, that there is room, next to the gravitational tension, for an additional refractional tension. In fact, the total potential interaction energy controls the bending of the light. The gravitational tension is proportional to the inverse distance to the center of the Sun. The refractional tension is frequency dependent and proportional to the third power of the inverse distance. We conclude this paper by showing that our model is in excellent agreement with historical electromagnetic deflection measurements.

Scaled Maxwell's Equations
We consider electromagnetic wave propagation with complex time factor exp(−i t), where i is the imaginary unit, is the radial frequency, and t is the time. In a source-free domain, with Cartesian coordinates x ∈  3 , we write Maxwell's equations in the frequency domain as where E = E(x, ) is the electric field vector, H = H(x, ) is the magnetic field vector, = (x, ) is the electric permittivity, and = (x, ) is the magnetic permeability. We neglect absorption, so that all material parameters are real valued.
Next we introduce the scaled electromagnetic field vectors as in which Z = √ ∕ is the electromagnetic wave impedance. Using this scaling in Maxwell's equations, we where the refractive index n = n(x, ) follows from in which c = 1∕ √ and c 0 = 1∕ √ 0 0 are the electromagnetic wave speeds in a material medium and in vacuum, respectively. Note that Equation 3 represents the Maxwell equations for the scaled electromagnetic field vectors in a "background" medium with permittivity 0 and permeability 0 . At this point, we may not assume that the waves in this background medium travel with the wave speed c 0 . The curl operators in Maxwell's equations are replaced by the medium-dependent curl operators, which determine the spatial dependency of the electromagnetic wavefield. Let us denote the fastest path of the wave as the geodetic line. For vacuum in the whole  3 we have constant electric permittivity = 0 and magnetic permeability = 0 (hence n = 1). Then, the curl operators are medium independent and the electromagnetic waves travel with wave speed c 0 along straight geodetic lines. In a vacuum subdomain in the vicinity of an object  with refractive index n ≠ 1, we are not allowed to conclude that the geodetic lines in that subdomain are straight. The geodetic lines are not equivalent to the raypaths in optics. The latter paths follow from a high-frequency approximation of Maxwell's equations. In the neighborhood of domain , these optical rays in vacuum remain straight when they pass , because within the ray approximation the interaction with matter in  is neglected. However, the presence of the object  leads to diffraction of the incident wave and this may influence the path of propagation. In fact, the geodetic line may become curved. Although, with the help of present-day computer codes a more or less complete solution of Maxwell's equations is possible, the structure of the geodetic lines is hard to observe from the numerical solution. We therefore investigate the nature of Maxwell's equations in a different coordinate system.

Maxwell's Equations in Tensor Notation
We introduce a Riemannian geometry with position vector x and symmetric metric tensor and its conjugate as, for example, Synge and Schild (1978):

10.1029/2019RS007021
Note that the Einstein summation convention with repeated indices is employed. In tensor notation, the scaled Maxwell's equations of (3) are written as the contravariant equations: where E i , H i , and i are the electric field vector, the magnetic field vector, and the partial derivative in the Riemannian geometry, respectively. These vectors are defined as The permutation tensor ijk is related to the Levi-Civita symbol as i k = e i k ∕ √ g, where g is the determinant of the metric tensor g ij . Using this relation, we obtain the covariant equations It is obvious that any solution of the Maxwell equations in a Riemannian geometry needs the specification of the refraction index and the impedance in whole space. Both material parameters occur in the curl operators.
To investigate the energy transport in the Riemannian space in more detail, we introduce the complex in which the asterisk denotes complex conjugate. Next we contract the complex conjugate of the left equation of (8) with E i and the right equation of (8) with H i * . Adding the two results and combining various terms, while using Taking the real part of (10) and multiplying the result with n √ g, we arrive at This equation represents the conservation law of energy transport in the Riemannian geometry. In the Cartesian space ( √ g = 1), this conservation law is written as Within the accuracy of geometrical optics, the curvature of the ray is determined by the refractive index only, see p. 114 of Born and Wolf (1959). In view of (11) and (12), we choose the metric tensor g ij to be a function of the refraction index n only.

Specification of the Metric Tensor
In this section, we make two specific choices for the metric tensor and investigate the consequences for the wave propagation.

Orthogonal Transformation
Choosing for the simple orthogonal transformation yields the diagonal metric tensor This choice of transformation just scales the local spatial behavior of the refraction index. It does not take into account the global refraction dependency. Later in this paper, we show that the geodetic line evaluation based on this metric leads to the well-known ray theory. In a vacuum domain outside the object , it leads to propagation along a straight line in the Cartesian space. Inside the object , this transformation leads to curved paths according to the standard ray theory based on a high-frequency approximation of the wavefield.

Nonorthogonal Transformation
If the refractive index is equal to 1 in the whole space (vacuum), the transformation is trivial (g = 1). As a consequence, the raypaths are straight lines, both in the Cartesian and Riemannian space. But we surmise that the presence of object  changes the geodetic structure in the vicinity of the object. Let us assume that g is a twice differentiable function not equal to one inside , then g ≡ 1 does not hold at all points in vacuum. This implies that g is a harmonic function determined by its values at the interface of the object. In other words, the presence of object  changes the direction of energy transport in each finite domain, because g ≢ 1. Physically, it means that a wave approaching the object  is "feeling" the object before it has reached it. It will propagate along a path where the variation of √ g Re{S k } is stationary. The direction will change at locations where √ g ≢ 1.
In this way the object shows its emergence to the incident wavefield, see also Feynman (1964), Chapter 26.5 A more precise statement of Fermat's principle.
In order to develop a transformation and hence a geodetic formulation that adheres the global dependency of the refraction index, we proceed as follows. The divergence of x k (x) (trace of the operator) is taken from the contraction of the second relation of (13) as Subsequently, we introduce the difference vector between the spatial points x k and x k as in which f k is a continuous and differentiable vector field. Next taking the divergence of f k and using (15), According the Helmholtz decomposition theorem, see p. 38 of Helmholtz (1858), a twice continuously differentiable vector field is uniquely determined by its curl-free component and its divergence-free component. Note that curl-free component of f is uniquely determined by (17). However, the divergence-free part of f is obtained as curl f = curlx, because curl x = 0. However,x is still unknown. Hence, a coordinate transformation can only be constructed if we require that the tension f is curl free or in other words that This means that the transformation matrix is required to be symmetric. After a differentiation of (16) with respect to x k , it follows that the new nonorthogonal transformation is given by Using the property that the tension f is curl free, together with the expression for the divergence of f, see (17), the Helmholtz decomposition theorem for a curl-free vector provides us the unique expression:

10.1029/2019RS007021
Obviously, f k is the tension due to the difference in refractive index with respect to vacuum. We define Φ as the refractional potential and we denote f k as the refractional tension. This representation is valid under the condition that n − 1 vanishes at the boundary surface of the domain . Equations 18 and 19 define our transformation matrix. They hold for any distribution of the refractive index inside domain . Note that the expression of the refractive potential yields a nonzero value outside  and this confirms that the refractive index distribution inside the object  determines the spatial coordinate transformation not only inside this object but also outside. Hence, the geodetic lines in the vacuum domain around  are influenced by the inner refractive index of the object.
Substitution the expression for the tension f into (18) yields the transformation matrix as In order to compare this transformation matrix with the one of (13), we analyze the domain integral on the right-hand side of (20) in more detail. When x ∈ , the evaluation of the domain integral has to be interpreted as the Cauchy principal value, where the contribution around the singular is excluded symmetrically and calculated analytically (Fokkema & Van den Berg, 1993). To this end, we consider the contribution of the integration over a spherical domain  with vanishing radius and center point x.
we observe that (20) becomes Here, f denotes that the integral has to be interpreted as its Cauchy principle value. We now assume that the refraction index is a slowly varying function in space. For decreasing distance |x−x ′ |, we observe that in the integral on the right-hand side of (22) the value of n(x) − n(x ′ ) vanishes. In addition, for increasing distance, the value of 1/|x−x ′ | vanishes. If we neglect the integral completely, we have established that the transformation matrix becomes identical to the one of (13). Hence, we confirm the well-known fact that the standard ray theory is only applicable for slowly varying refraction index and for locations far away from significant changes in the refraction index. In conclusion we remark that on the right-hand side of (22) the first term represents the local and orthogonal part of the transformation, while the second term stands for the global part.

Construction of the Geodetic Line
In the construction of the geodetic line we consider the scalar arclength ds along the curved geodetic line, given by At this point, we switch to the matrix representation of the tensors and introduce the curvature matrix  i as a representation of the symmetric transformation tensor x ∕ x i of (18), namely, Since the matrix  k is real and symmetric, an eigenvalue decomposition with positive eigenvalues exists and the sum of the eigenvalues is equal to the trace. By inspection of (18), we learn that in the transform geometry the coordinate axes are spanned by the components of the tension f. The latter is directed in the normal direction to the surface Φ = constant. This normal is an eigenvector belonging to the matrix  i , because

10.1029/2019RS007021
In our further analysis we choose as the principal unit vector of the transformation, complemented with the other two eigenvectors 1 and 2 with corresponding eigenvalues 1 and 2 . Since the trace of a symmetric operator is equal to the sum of the eigenvalues, from (15) it follows that + 1 + 2 = 3n. The location of the local frame is such that the vector = { 1 , 2 } is tangential to the surface Φ = constant, while is perpendicular to it. Next we consider the scalar arclength ds of (24). Using the eigenvalue decomposition, we may write Introducing the unit vectorŝ To investigate the dynamic behavior, see p. 114 of Born and Wolf (1959), we consider the optical length of the geodetic path, which is defined by the actual length of the path times the index of refraction. Hence, the left-hand side of (27) represents the optical length of the path. Therefore, we introduce the virtual refractive index n g along the geodetic path as In general, the eigenvalue decomposition has to be determined numerically, except for a rotationally symmetric medium or a horizontally layered medium. In the latter case, the eigenvalues are identical and equal to n. Then, the ray theory applies, because n g = n and the horizontally layered medium has no curvature.
The virtual refractive index n g (x,ŝ) controls the path of the geodetic line in a similar way as the refractive index n(x) controls the path of optical rays. Note that the virtual refractive index is not only determined by the local position of the geodetic line but it also depends on the direction of the geodetic line at this position. We construct this geodetic line by considering the classic differential equation for the evolution of an optical raypath, see p. 121 of Born and Wolf (1959), but we replace the physical refractive index n by the virtual counterpart n g , namely, where x j = x j (s) is the trajectory of the geodetic line and s is the parametric distance along this trajectory, whileŝ is the tangential unit vector along the geodetic line. We note that this differential equation applies to refractive indices, which are invariant for the direction of the geodetic path. However, the explicit Euler integration of this differential equation updates the ray position and ray direction, so that only the previous information of position and direction about the associated path segment is used. The path directions are taken to be constant during each integration. This is consistent in keeping the refractive index constant during the update step.

Radially Inhomogeneous Medium
At this point we use spherical coordinates. For a rotationally symmetric configuration, the eigenvalues can be determined analytically. The refractional potential is given by The tension depends on R only and is directed in the radial direction. Hence, = = 0 and the radial component R = R Φ is given by

10.1029/2019RS007021
From (25) it follows that the eigenvalue in the radial direction is given by while the eigenvalues in the tangential directions follow from the trace, In view of the axial symmetry of our configuration, the tangential eigenvalues are the same. We therefore confine our analysis to the plane in which the geodetic path is defined. Hence, we suffice with the computation of , namely, The eigenvalues depend only on R f R . From (17) we observe that From (32) it follows that and The virtual refractive index is then obtained as, compare (28), whereŝ R = cos( − ) andŝ = sin( − ) are the components of the unit vectorŝ. Here, is the angle betweenr and the x 1 direction, while is the angle betweenŝ and the x 1 direction.

Numerical Reconstruction of the Geodetic Paths for a Sphere
We consider a sphere with a radius a = 60 m. In the inner domain, the refractive index is n = 1.5. To avoid numerical problems, we require that the refractive index is a continuously differentiable function as function of R. We apply a cosinus type or tapering of the refractive index at the boundary region. As width of this region we take Δa = 0.01 m. Within this region, the refractive index varies from 1 to 1.5, for decreasing R, namely, Differentiation with respect to R yields R n(R) = −0.25 Δa sin Note that this function is continuous for all R; it vanishes everywhere, except in the small boundary region of 0.01 m width. We remark that, for the present refractive index, the integral of the tension f R is calculated analytically. The numerical construction of the geodetic path is performed by an Euler integration of (29). In order to have enough integration steps for a geodetic line, passing the boundary region of the sphere, we take an integration step of 0.1Δa. The starting point of all geodetic paths is x S = (−120,0, 0) m. To facilitate the computation of the wavefront curvature propagating along a set of geodetic lines in the plane x 3 = 0, we compute a large set of geodesics starting with an angle of 0.1 • between each other. For each path in the plane x 3 = 0, we store the locations and the travel times to reach each location. The wavefront curvatures are computed by connecting the points of equal travel time of the various paths. In Figures 1 and 2, we show some geodetic lines in the plane x 3 = 0 of the Cartesian space and some wavefronts for a travel time of t = 0.25, 0.50, and 0.75 μs, respectively. In Figure 1, we present the raypaths using the orthogonal metric tensor of (14). In other words, in our computer code we have replaced the virtual refractive index n g with the actual index n and we get the results of the standard ray theory. There are two remarkable observations. First, there is always a "shadow zone," where there is no wavefront. In this domain the wavefront of the arrival time seems to disappear abruptly. Second, at the points where these shadow zones arise, there is a duplication point. This is a point where a very small displacement of the raypath results in large refraction. Note that this phenomenon occurs despite the fact that the refractive index is continuously differentiable.
These types of artifacts are due to the local character of ray theory, which is based on a high-frequency approximation of the wave equations, in which the object does not "feel" the arrival of a passing wave.  In Figure 2, we present the geodesics using the virtual refractive index n g . First, we note that there are no shadow zones and that all geodesics are bent outside the sphere. This bending of the geodesic line is stronger for a path closer to the sphere. It is obvious that the wave inside the sphere travels slower than outside. Given the continuous property of n g , the wavefront cannot break and travels in a rather peculiar way along the interface of the sphere. One may wonder to what extent this theoretical model of wave propagation along geodesics describes the actual situation of wave propagation. In the next section, we discuss a validation against the numerical solution of a full-wave integral-equation model as solution of Maxwell's equations.

Verification Using Contrast-Source Integral Equations
In the case that the magnetic permeability is equal to the vacuum permeability, the unknown electric-field vector E(x, ), for x ∈ , may be obtained from a contrast-source type of integral equation. For a frequency independent permittivity, we define a contrast function and a contrast source w E as respectively. We then have the integral equation for the 3-D unknown contrast source distribution w E ; see, for example, Zwamborn and Van den Berg (1992) and Abubakar and Van den Berg (2004): Here, E inc (x, ) represents the known electric-field vector of the incident wave in the whole space, in absence of the contrasting domain . In the present case we take the electric-field distribution of an electric dipole oriented in the x 1 direction. Furthermore, the Green function G is given by For the numerical solution of the integral equation we define a regular grid of square subdomains. This grid includes the scattering object . At grid points outside , we enforce the contrast function to 0. The integral equation is solved iteratively, using the BiGSTAB method developed by Van der Vorst (1979). In view of the convolution structure of the integral operator, the operator is computed efficiently by using fast Fourier transform ( In Figure 3, for t = 0.5 μs, we show this power distribution, both for the incident wavefield E inc and for the total wavefield E. In the left picture, the center of the spherical wavefront is the blue curve. This is the reference of the position of the wavefront. Note that the wavefront vanishes at x 1 = 0, because the dipole is oriented in the horizontal direction. The right picture of the total wavefront is now used as benchmark of the wavefronts predicted by the ray theory and the geodesic theory. We make similar images for t = 0.25 μs and for t = 0.75 μs and use the three images as overlays over those of Figures 1 and 2 to arrive at Figures 4 and 5. Comparing the middle pictures of Figures 4 and 5, we observe that the inner wavefront is delayed and an extra bending of the wavefront occurs to bridge the difference in wave speed along the interface of the sphere. The standard ray method leads to a shadow zone at the boundary of the object, while the geodetic theory indeed predicts this extra bending. In the right picture of Figure 5, the extra curvature of the actual wavefront, consisting of two wavefronts with different wave speeds, is predicted by the geodetic theory. The transition between outer and inner wavefront is bridged by a wavefront propagating along the interface. Physically, we interpret it as a surface wave which is launched along the interface of the sphere; see Figure 6. Note that the wavefronts in ray theory are orthogonal to the raypaths, but in the latter figure it is obvious that the wavefronts in the vicinity of the object are not orthogonal to the geodetic lines. In the three pictures, we observe that the wavefronts at the boundary of the object propagate along this boundary, with an amplitude decaying in the perpendicular direction: acting as a surface wave.

Additional Bending of Radio Waves by Sun's Refractive Index
We now turn to the consequences of our findings. They shed another light on the bending of electromagnetic waves while closely passing an object with a contrasting index of refraction. Fokkema and Van den Berg (2019) discussed the phenomenon that a light ray passing the Sun has an additional deflection outside the object, which is controlled by the refractional potential in a similar way as the mass-density potential changes the path in the theory of gravity which was shown by Einstein (1911). This means that the total potential energy is a superposition of a gravitational and an electromagnetic constituent of Sun's interior. The electric permittivity and magnetic permeability determine the velocity of light c 0 in vacuum. In material media the electromagnetic wave speed c is less. As we have shown this permits us to characterize the medium by its index of refraction n = c 0 ∕c. The deflection intensity is constituted by the refractional deflection which is frequency dependent and proportional to R −3 (Fokkema & Van den Berg, 2019), while the gravitational deflection is frequency independent and proportional to R −1 . In the remainder of this article we look at historical data for evidence for our claim.

Validation on Historical Data
An overview of optical deflection measurements are given by Von Klüber (1960), Will (2015), and Shapiro (1999). Based on the gravitational model, the deflection angle is given by d GR = R  ∕R, where = 1.75 (in arcsec) is the Einstein value and R  is the effective radius of the Sun. Mikhailov (1959) analyzed in detail the six eclipses during the period 1919-1952. These historical optical deflection measurements are tabulated by Pathria (2003), and he concluded that the spatial dependence is correct, but the spreading of around the Einstein value is significant in the near region of the Sun. Shapiro (1964) and Shapiro (1971) suggested that a more accurate deflection measurement follows from radio interferometry. In radio experiments, Sun's corona effects the (frequency-dependent) deflection to a larger extent than in the optical experiments. Seielstadet al. (1970) showed discrepancies up to 20% in neglecting the coronal effects. Muhleman et al. (1970) incorporated the coronal plasma effect and observed a spreading from 10% to 15%. Later radio experiments confirm this frequency dependence, while the spatial variation differs from the inverse-distance relation. This was explained by extending the GR model with the local bending due to the frequency-dependent coronal medium. However, satisfactory fitting to the measurements was only possible in a restricted range of R; see, for example, Figure 1 of Merat et al. (1974). Fokkema and Van den Berg (2019) investigated the optical-deflection data collected by Merat et al. (1974) and showed that the fitting to the measurements over the whole radial range of observations improved substantially, once the additional bending by Sun's interior refractive index is taken into account. In the model the frequency dependence of the data has not taken into account. In this paper, we consider the radio deflection data, which are certainly frequency dependent. At this point, we return to the work of Merat et al. (1974). They conclude on basis of radio deflection observations made by Muhleman et al. (1970), that for R < 5 R  deviations from the Einstein prediction become statistically significant. They have collected the whole set of radio deflection measurements into four samples; see the fifth column of Table 3 of their paper. The weighted mean of the distance R/R S has been given, together with the range of deviations of every measurement. The deviations from the GR effect are significant for R/R S < 5. But that is mainly due to the frequency dependency of the data. After subtraction of the gravitation term, 1.75 R S /R, we obtain the electromagnetic constituents of the radio deflections and denote them as d EM . Since we surmise that the upper and lower bounds are related to different frequency ranges, we consider the upper and lower bounds of measured deviations separately. They are denoted by their superscripts. Hence, we have two sets of four data points, namely, d EM ∶= d upper and d EM ∶= d lower , respectively.

Influence of the Naked Sun
Let us first consider the additional bending of an electromagnetic wave passing the Sun, and we neglect the presence of the corona. Following the pure gravity light-bending theory of Maccone (2009), we also denote this as the naked-Sun situation. For small values of the refraction index of the Sun, Fokkema and Van den Berg (2019) have shown that the electromagnetic deflection is asymptotically given by To find the unknown factor B from the four data points d EM ∶= d upper , we carry out a least squares fit, which minimizes the residuals where R 0 is the smallest value of R on the geodetic line. The value of R 0 ∕R  is often denoted as the impact parameter. The minimum residuals are given in the third column of Table 1. The mean of these residuals   (46), the deflection function d EM (R) is presented as the solid blue line in the left picture of Figure 7. A similar procedure for the four data points d EM ∶= d lower is carried out. The minimum residuals are given in the third column of Table 2. The mean square of these residuals amounts to 7.3%. Substituting the resulting value of B in Equation 46, the deflection function d EM (R) is presented as the solid blue line in the right picture of Figure 7. We now observe that the deflection curve is negative. This is typically the effect of plasma of the outer region of the Sun. The large discrepancies of the two curves with the measured data may be explained by the "coronal mantle" outside the Sun.

Influence of the Coronal Mantle
In the corona, we only take into account the local effect of the refractive index of the corona. In order to include the plasma effects of the corona, we start with the refractive index described as a superposition of powers of R  ∕R, with constant factors p , namely, The data under consideration are obtained for R > 3R S , and we employ the refractive index described in Muhleman et al. (1970), namely, where p 1 = 6 and p 2 = 2.33. We conclude that the electromagnetic deflection may be written as For the range of R > 3R  we determine the coefficients C p 1 and C p 2 by a least squares fitting of (50) to the four data points. For the upper bounds we define the residual error as The minimum residuals are given in the fourth column of Table 1. The mean of these residuals amounts to 6.3%. Substituting the resulting value of the coefficients C p 1 and C p 2 in (50), the deflection function d EM (R) is presented as the dashed blue line in the left picture of Figure 7. A similar procedure for the lower bounds of the data yields the dashed blue line in the right picture. The discrepancies with these data points are presented in the fourth column of Table 2, with a mean error of 3.3%.

Influence of the Naked Sun and the Coronal Mantle
For small deflections, we take a linear superposition of the naked-Sun part and the mantle part (the corona). We conclude that the total electromagnetic deflection may be written as When we apply a least squares fitting procedure of this function with three unknown coefficients to four data points, we observed that the system matrix is heavily ill posed and impossible to invert numerically. A stable result is obtained by preconditioning. We rewrite (52) as 10.1029/2019RS007021 where C 1 = C p 1 ∕B and C 2 = C p 2 ∕B. This nonlinear equation is solved with an iterative Gauss-Newton method. As starting values we take zero values for C 1 and C 2 and determine B by a direct least squares minimization. After carrying out a few Gauss-Newton iterations, a stable result is obtained. The minimum residuals are given in the fifth column of Table 1. The mean of these residuals amounts to 1.0%. Substituting the resulting value of the coefficients C p 1 and C p 2 in (52), the deflection function d EM (R) is presented as the red line in the left picture of Figure 7. A similar procedure for the lower bounds of the data yields the red line in the right picture. The discrepancies with the data points are presented in the fifth column of Table 2, with a mean error of 0.2%.
Under condition that we keep the GR prediction unchanged, we claim that the near-field correction due to the tension of the Sun's interior refractive index is a prerequisite to obtain an accurate model in solar gravitational lensing; see, for example, Eshleman (1979) and Maccone (2009).

Conclusions
The propagation of electromagnetic energy over the fastest paths is investigated: (1) using the standard ray theory and (2) a novel approach based on the theory of geodesics. The analysis of raypaths showed that there is always a "shadow zone." Moreover, where this zone arises there is a duplication point. These types of artifacts are due to the local character of ray theory and are not present in the geodetic description, which has a global character. In addition, the geodesics bend in the vicinity of the object and the wavefronts are nonorthogonal to the geodetic lines. At a curved boundary of the object it predicts the propagation of surface waves, where both the wavefront and geodetic propagate parallel to the boundary surface. For a spherical object, the geodetic wavefronts are verified using a full 3-D numerical simulation based on a contrast-source integral equation. The conclusion is that the present theory of geodesics offers a reliable physical insight in the actual propagation of electromagnetic waves.
The theory of geodesics has consequences for the explanation of the light bending around the Sun; namely, next to Einstein's gravitational tension there is room for an additional refractional tension. With this extended model historical "radio light" deflection measurements has been investigated. The conclusion is that it explains these measurements very well. It adds a significant correction to solar gravitational lensing and interstellar radio communication. On 12 August 2018, the Parker Solar Probe mission (http:// parkersolarprobe.jhuapl.edu) has been launched with dedicated instruments for measuring the electromagnetic fields and two-way radio transmissions with the Earth station at different frequencies (Sokol, 2018). This mission will create excellent conditions for collecting the electromagnetic properties of the Sun.

Data Availability Statement
All numerical results can be reproduced by using data and information available in the listed references, tables and figures included in the paper; in particular the geodetic lines are constructed numerically via a predictor-corrector version of the recursive scheme of (38) of Fokkema and Van den Berg (2019).