Radiation Field Inside the Homogeneous Moon Excited by a Vertical Magnetic Dipole

The investigation of radiation field in the moon is important for the moon subsurface detection and vertical magnetic dipoles (VMD) are especially useful as sources and probes. The present work focuses on the radiation field inside the moon excited by a VMD located on or near the moon’s surface. Analytical expressions are derived with approximate method. Contributions of the fields come from two parts: the series expansions of low-order terms using the multiplication iteration method, and the integrals of high-order terms using the asymptotic method. Meanwhile, the field of high-order terms is close to the total field. Compared with the calculation’s direct sum, our method has advantages of accuracy and is suitable for the large-size model. For a source-receiver system on or near the homogeneous moon, it is found the waves propagating inside the moon can be described in terms of ray optics. Based on the formulas, various modes of reflection are taken into account, and the intensity of the primarily reflected waves is calculated numerically. Furthermore, effects of propagation distance and angel, the operating frequency, as well as the height of the source are considered to demonstrate the properties of the field. Particularly, the proposed method and formulas have applications to the moon interior exploration.


I. INTRODUCTION
As the nearest planet to the earth, the research of the moon contributes essentially to our scientific understanding. In particular, it is imperative to develop space exploration, outline the stratigraphy structure of the lunar subsurface, and detect related parameters of the moon. Generally, in the probing of the moon, the radar detection technique has been employed, which was used for the first time during the Apollo 17 mission (Apollo Lunar Sounder Experiment, ALSE) to characterize the lunar structure [1]. Since then, experiments have continued in directions of the moon interior exploration. Chang'e-3 (CE-3), Chang'e-4 (CE-4), and Chang'e-5 (CE-5), the first attempt to survey the moon's surface using the Lunar Penetrating Radar (LPR), have landed successfully on the moon in 2013, 2018, and 2019, respectively [2]- [4]. The LPR, equipped with a low-and high-frequency antenna, has advantages of reducing the electromagnetic (EM) wave loss caused by space propagation [5]. Meanwhile, noises interference from a complex environment of the moon can be The associate editor coordinating the review of this manuscript and approving it for publication was Wei Xu . decreased. Hence, investigations of the LPR for the moon detection has significance for high-precision measurement. Notwithstanding contributions of experiments or the observation data acquired from the lunar [6]- [8], however, the theoretical treatment for an antenna near the moon's surface is still lacking, which motivates us to study this problem.
Generally, the transmitting antenna mounted on the LPR can be regarded as a dipole antenna. In the history of propagation mechanism over a sphere excited by a dipole, landmarks have been Waston's series expressions in the presence of a perfectly conducting sphere with a large radius compared to the wavelength [9]. However, Wastons's results have been left inconvenient for engineering use due to the series converges slowly. Fock [10] optimized the series of Waston to assemble faster using transform and approximate methods afterward. In the case of a sphere with finite conductivity, extension achievements were made by many investigators, especially including Van der Pol and Bremmer [11], Gray [12], North [13], and Wait [14], [15]. A more practical development was introduced by Wait [16] using the mode expansion method, of which the boundary conditions at the sphere surface were specified by surface impedance. Detailed findings were summarised in two monographs by Bremmer [17] and Wait [18], which presented asymptotic formulas in Airy function forms used to calculate the ground wave of vertical electric or magnetic dipoles (VED/VMD) over the homogeneous spherical earth. In connection of a similar method, further investigations were extended by Zhang and Pan [19] and Li et al. [20], [21] for the spherical surface of the earth with N-layered dielectric. Meanwhile, some investigators have treated the EM field over a sphere excited by a dipole under different method. Taha et al. used differential transform method to investigate the radiation field due to a dipole in the earth considering the effect of atmosphere electrical conductivity [22]. Recently, Dyadic Green function technique was applied to examine the scattering field of a perfect conducting sphere [23], [24].
Analytical approximations for the EM field have been obtained using methods mentioned above; however, not straightforwardly getting the exact closed-form solutions and overlooking the areas inside the sphere. Owing to this reason, Wu [25] invoked the concept of creeping waves based on the Poisson summation formula. Meanwhile, Nussenzveig [26], [27] applied a resemblance analysis to describe waves inside a sphere. In the same idea, Houdzoumis [28], [29] led to a new solution in terms of a series of integral that is equivalent to the exact series solution. EM fields excited by a VED over a homogeneous sphere were derived when both the source and observation point were located on the interface. Essential differences between EM waves propagation inside and outside the sphere are: 1) EM waves exponentially decreasing through the free space; 2) Ray contribution to the waves bouncing and circulating inside the sphere. A creeping-wave structure for all six components along the boundary was revealed by Margetis [30], which studied the EM field in the free space due to a horizontal electric dipole (HED) located below and tangential to the surface a homogeneous sphere. When considering the source located some distance away the surface, the results of Houdzoumis and Margetis would be invalid. Convergence problem appears when the source is off the interface. More recently, Houdzoumis [31] reexamined developing formulations for any heights of the dipole above the sphere. The emphasis was placed on expanding the integrals by dealing with poles. Unfortunately, no explicit approximate analytical formulas were presented, only computation for special cases.
Along with this research line, we reconsider same problem through a new approach. The background of this research can be divided into the lunar exploration in low frequency ranges and high frequency EM scattering analysis [32], [33]. As is mentioned above, we pay attention to the moon subsurface detection in this paper. The purposes of this paper are twofold: 1) analytical expressions are derived due to a VMD over the moon whose interior is homogeneous. Of primary interest is the field when the source is off the moon's surface; 2) propagation characteristics affecting electric field are evaluated under different conditions. In the presented treatment, a practical procedure is developed to solve the convergence problem of the series. It is accomplished through a new theoretical method, which is less restricted to the height of the source. With the propose method, the computation becomes less tedious when the sphere radius is larger than the wavelength. Firstly, the series expansions are divided into low-order terms and high-order terms. For low-order terms, we adopt the multiplication iteration method. Then, we consider the asymptotic method for transforming series into integrals for high-order terms. Finally, in view of the expressions of integrals, waves propagating inside the lunar can be described in terms of ray optics using the phase-state method.
The remainder of this paper is organized as follows: in section II, radiation fields are analyzed with the approximate method. The entire field expressions include series expansions of low-order terms and integrals of high-order terms. In section III, numerical results are presented. Some examples are considered to demonstrate the properties of the high-order terms contributions. In addition, the propagation distance, the propagation angle, and the height of the source are taken into consideration, as they affect the properties of the radiation field. In section IV, important conclusions are drawn, a discussion is presented.

II. METHODOLOGY
The geometry of spherical coordinates is presented in Fig.1, which consists of a spherical and electrically homogeneous moon of radius a. The center of the moon is the origin of a coordinate system (r, θ, ϕ). A VMD is located at (b = a + h b , 0, 0), and the observation point is taken at an arbitrary point P(r, θ, ϕ). The distance between source and the observation point is R = √ b 2 + r 2 − 2br cos θ. We shall assign subscripts 0 and 1 to the parameters of the free space and inside media of the lunar. Electric constants of the free space and the moon are taken as (σ 0 , ε 0 ) and (σ 1 , ε 1 ), the permeability being taken as µ 0 for both. The magnetic dipole may be realized by a small loop of wire carrying a uniform current Ida·e −iωt , where ω is the angular frequency of the current. Wave numbers of two regions are expressed as: According to the work [34], the conductivity σ 1 of the moon is small enough, and k 1 can be treated as real for the sake of simplicity. It is assumed that |k 1 | 2 |k 0 | 2 and k 0 a 1. In addition, radius of the lunar is larger than the wavelength.
In a nutshell, the physics can be described as follows: As the EM field propagates along the sphere surface, the direct field is excited when the EM reach the observation point. Simultaneously, due to the induction of the dipole source, the secondary field is generated, which exists in the entire space. The secondary field is divided into two parts: radiating and scattering field in the free space, and the field is inside the sphere in the form of rays.

A. SERIES EXPANSIONS OF RADIATION FIELD EXCITED BY A VMD
The EM field in two regions (j = 0, 1) satisfy the Maxwell equations where Since the source is a VMD, it is apparent from symmetry that only B jr , B jθ , E jϕ are to be considered. Electric fields radiated by a dipole on the sphere have been investigated by Houdzoumis [28] and Margetis [30] in detail. According to their works, with no diffraction conditions, direct waves E in 0ϕ due to a VMD in the free space are expressed as follows: Here, j m is the spherical Bessel function. h (1) m refers to the spherical Hankel function, and P 1 m are the Legendre function of the first kind. It is natural to set where E sc 0ϕ represents the scattering field in the free space and the E 1ϕ is the total field inside the lunar. T m is the transmission coefficient and S m is the scattering coefficient.
It is then deduced from the boundary that at the interface of the surface, the tangential component of the inside and the outside electric fields are equal, i.e., According to the boundary conditions, the coefficients are obtained The scattering field in the free space and the total field in the moon are obtained by subscribing (6)-(8) into (4) and (5).
We have

B. MULTIPLICATION ITERATION METHOD FOR LOW-ORDER TERMS
From equations above, it is straightforward to obtain exact series solutions for the fields. However, the series expressions converge slowly when the radius of the moon is much larger than the wavelength in the free space. In this section, a new solution has been proposed in response to this problem. In the solution process, terms m are divided into low-order terms and high-order terms for calculation. Due to the large size of the moon, even in the low-order terms, the Bessel function and the Hankel function will turn to be large or small enough that over the range for the computer calculation. Therefore, for low-order terms, the multiplication iteration method of the Bessel and Hankel function is presented to solve this problem. The solution is where j m (x) is the Bessel function, j m (x) refers to the derivation of the Bessel function. α, β, γ and η are iteration coefficients, respectively. The specific derivations are shown in Appendix. We have Eqs.(11)- (14) are also applicable for the Hankel function.
For the sake of generality, j m (x) (j m (x)) is used for the Bessel function iteration matrix and h m (y) (h m (y)) is for the Hankel function iteration matrix. We write and Multiplying (15) and (16), recurrence formulas are expressed in matrix form by transposing. This gives where W refers to the recurrence coefficient matrix, When m = 0, the initial values of the multiplication iteration functions are in the following form Grouping Eqs. (19)- (22) into (15)-(18), the entire function sequence of low-order terms is recursively deduced. We have where, m c is the critical term of low-order terms and highorder terms. For low-order terms of the electric field, the terms vary from 0 to m c . The variations of the electric field E 1ϕ of low-order terms along the number of terms are shown in Fig.2, where the solid line refers to the results calculated by the proposed method, and the dashed line represents the results obtained by the calculation directly the sum of the series algorithm. The calculation is both under the condition that shown in Section III. Taking the operating frequency f = 60 kHz, k 0 a = 2184, and m c is 6000. From Fig.2, it is found that two methods agree well. The difference between the method proposed in this paper and the traditional recursion method is 2 dB. With the presented method, the resonance peak appear later than the direct sum method. Therefore, the multiplication iteration method can calculate more terms m. After calculation and verification, the method proposed in this paper has higher stability and calculation accuracy than traditional methods and is suitable for the large-scale model.

C. ASYMPTOTIC METHOD FOR TRANSFORMING SERIES INTO INTEGRALS FOR HIGH-ORDER TERMS
In the case of high-order terms (i.e., m = [m c + 1, ∞]), asymptotic integral methods are applied. Asymptotic expressions of the Bessel and Hankel function are derived as follows [35]: Combining (24) and (25), the electric field of high-order terms is obtained In the remarkable paper by Wu [25], the Poisson summation formula was invoked to analyze radiation waves inside a sphere. We have As a result, Eq.(26) is replaced by integral sums: In particular, the last term of Eq.(28) is series expressions, which can be calculated through direct iteration.
Applying relation between the Bessel (Hankel) function and the half-order Bessel (Hankel) function, we have Therefore, C(λ) is defined by

1) GEOMETRICAL RAY WAVES INSIDE THE MOON
According to the characteristic of Bessel functions, we have Then, Eq.(31) can be transformed as From above formulas, [1 − E(λ)] −1 can be expanded into series and (34) is rewritten as It is obvious that E 1ϕ of high-order terms has two parts, the first part is the wave propagating along the surface of the moon depending with the poles, the second part represents behavior of the wave propagating within the moon.
For the first term, the integral oscillates due to the P 1 aλ−1/2 (cos θ) and the exponential e i2πnaλ . Concerning the second group of terms, the factor F(λ)E M (λ) is of oscillatory nature, as well. Moreover, the oscillation of this factor interferes with P 1 aλ−1/2 (cos θ) and e i2πnaλ in a way that can render the phase stationary.

2) PHASE STATE METHOD
According to the work [36], the optics-ray integral can be evaluated by the method of stationary phase. The solutions by the method of stationary phase agree with those obtained using numerical integration method. In the following, wave propagating through the lunar will be mainly considered and analyse their contributions. We have Using the approximations of Bessel and Hankel functions for aλ 1 [37], that is where To examine the roots of stationary phase, we pay attention to the phase part of the field inside the moon. This gives E nM ± = k 1 a 2(M +1) 1−(λ/k 1 ) 2 −(λ/k 1 ) cos −1 (λ/k 1 ) where the ''±'' sign originates form the Legendre functions whose expression involves a cosine.
The above functions represent the part of the phase that depends on λ, we have Therefore, the phase becomes stationary when For λ = k 1 cos ϕ nM ± , we have here ϕ nM ± is the angle between the ray and the corresponding local tangent. Subscribing (43)-(49) into (39), the preceding consideration lead to the field E ray 1ϕ inside the lunar, which is addressed as follows: ×exp{[i2(M + 1)k 1 a sin ϕ nMs + k 1 a cos ϕ nMs ln where S M + = n : 0 ≤ 2n ≤ M , S M − = n : 0 ≤ 2n ≤ M + 1. G(k 1 cos ϕ nMs ) is the distance covered by rays which start from the north pole forming an angel. R M (k 1 cos ϕ nMs ) equals the reflection coefficient. The rays described by (50) circulate around the origin while they are multiply reflected at the spherical boundary, as depicted in Fig.3. Waves are determined by all possible pairs (n, M , s), of which n is the number of circulations, M presents the reflection number, s refers to the clockwise (+) and counterclockwise (−). Primary waves from source and various reflected few times at the spherical surface are superposed together. For simplicity the source is assumed to oscillate with no diminution. The diminution of the radiation field will become remarkable only after once or twice reflected.
So far, the total field excited by a VMD inside the moon are obtained. It can be readily seen that different terms of series are described in series and integrals, respectively. ×exp{[i2(M + 1)k 1 a sin ϕ nMs + k 1 a cos ϕ nMs ln In order to verify corrections of the proposed method in this paper, computations are carried out using the traditional method (i.e., direct sum method) and the fast multipole method (FMM) which is addressed in the work [38]. The operating frequency, critical term, and the dielectric constant inside the sphere are taking as f = 3 GHz, 1 = 2, and m c = 300, respectively. Assuming both the dipole and the observation point are on the surface, the radius of the sphere equals 1. Comparison results are shown in Fig.4 and Table 1.
For sake of clarity, the normalized propagation distance η = ρ/ρ c is introduced by Houdzoumis [29], and ρ c = a( k 0 a 2 ) −1/3 . From different methods, it is found numerical results by the method proposed in this paper are agreement to the direct sum method and the FMM. Particularly, the method proposed in this paper has a faster convergence rate. It is noted that from Table 1 that most of the relative error of results between the direct sum method and the proposed method is less than 6%, whereas those between the direct sum method and FMM is much larger than 11 %.

III. COMPUTATIONS AND DISCUSSIONS
According to above derivations and analyses, radiation field components inside the moon excited by a VMD are studied in detail. In this section, radiation fields will be treated numerically via MATLAB. In this paper, we adopt the cold   homogeneous moon model which has been proposed by Ward [39]. Cold homogenous model is characterized by frequency dependent effective values of dielectric constant and conductivity which are uniform throughout the lunar. Parameters used in calculating are indicted in Table 2. According to the work [34], the conductivity of the interior moon close to the surface is about σ 1 = 10 −12 . Therefore, the imaginary part of k 1 is smaller than the real part, which can be ignored in the low frequency (LF) ranges. Moreover, the observation point is located inside the moon and close to the surface (i.e., r → a).
Firstly, Fig.5 shows the comparison between the results of this paper and those by Pan [21]. The operating frequency is f = 150 kHz and other parameters are taken from Table 2. Apparently, the analytical curve of this paper and Pan's curve variation tendencies are approximate and unanimous, completely confirming the high accuracy of our method. In addition, the phenomenon of interference arises in the curve of this paper, which shows our method has advantages for computation.
To illustrate the attenuation of the electric field depending on the propagation distance, single bounce rays are taken as examples (i.e., n = 1, M = 2 and n = 2, M = 2). The radiation component E 1ϕ varying with normalized propagation distance is shown in Fig.6. Different cases are chosen with f = 60 kHz, f = 80 kHz, and f = 150 kHz to gain insight into the influence of frequency. According to Fig.6, there are three sets of curves in which a comparison is made among the low-order terms, the high-order terms, and the total field. For a large-size model, the high-order terms of electric field play an important role in the total field. It is evident that as the normalized propagation distance increases, the electric field intensity decreases. A phenomenon of interference appears when the normalized propagation distance exceeds η ≥ 0.3. When η ≥ 3, the interference will strengthen, which reaches the peak at the observation point. It is caused by the superposition of waves of different paths and various reflections at the observation point. Comparing the results of different frequencies, electric fields increase as frequencies advance.
For same parameters values, typical computations of electric fields E 1ϕ are carried out versus the propagation angle with the operating frequencies f = 80 kHz and f = 150 kHz, respectively. From Fig.7, the electric field trend varying with the propagation angle is similar with that of the propagation distance. For a large-size model, the propagation distance and the propagation angle have a following relation: ρ = aθ. As is shown in Fig.7, the amplitude of E 1ϕ decreases when the propagation angle increases and electric fields increase as the frequencies increase. It is noted that the first reflection angle appears early as the (n, M,s) enhances. Fig.8 demonstrates relations between the height of the source and the electric field. The source is located at h b = 30 m, h b = 10 m, and h b = 0 m, respectively. Theoretical results are shown at the operating frequency f = 75 kHz. Under different heights of the source, the electric field highly oscillates with the normalized propagation distance increasing. In the three cases, the regular patterns of electric field  varying with normalized propagation distance are similar. It is due to the height of the source is much smaller than the radius of the moon, so the difference of electric field at different heights is slight. Meanwhile, a distinction arises in the electric field generated by a VMD in the moon. When the source is located on the boundary, the amplitude of electric field is larger than that off the surface.
Spatial distributions of electric field E 1ϕ at the operating frequencies f = 60 kHz, f = 100 kHz, and f = 150 kHz are presented in Fig.9, respectively. According to Fig.9, the electric field in the moon show a divergence distribution in space. It demonstrates the amplitude of the electric field is relatively large near the source, and the field intensity decreases as the propagation distance increases. In particular, the minimum of the E 1ϕ aggregates at the location near the source. This is probably because the minimum of the E 1ϕ corresponds to the change point of the phase. By comparing the results at two frequencies, the radiation pattern is similar. It is noted that the frequency influences the strength of electric fields, where the amplitude of E 1ϕ increases as frequency increases.

IV. CONCLUSION
In this paper, the EM fields inside the moon, which are radiated by a VMD located on or near the surface, have been proposed in series form of low-order terms and integral form of high-order terms. To solve the converge problem of a largesize model, the multiplication iteration method for low-order terms and the asymptotic method for high-order terms are employed to formulate the expressions of the field. When the source is on or near the boundary, the integral expansions of high-order terms dominate the total field. Comparison results show that the radiation fields against the direct sum method with the relative error of 6%. The proposed method can give a better performance in accuracy for a large-size model. It is found that field attenuates strongly with the propagation distance, the propagation angle, and heights of the source increasing, but phenomenon to frequency is electric fields increase as frequencies advance. Meanwhile, the interference appears when the normalized propagation distance exceeds η ≥ 0.3. When η ≥ 3, the interference will strengthen, which reaches the peak at the observation point. The analytical method presented in this paper provides theoretical support for the moon subsurface exploration.
The present work admits several extensions and also points to pending issues. Results presented in this paper are only for analytical and numerical computations. Experiments are needed to demonstrate the theory. In actual situations, the environment inside the moon is complex, which would require more complicated integrals by imposing some restrictions. These issues will be investigated in our future work.

APPENDIX
According to the recurrence relation of the Bessel function, Therefore, iteration coefficients α, β, γ and η are obtained According the work of [40], when m = 0, the initial values of the spherical Bessel function and the Hankel function are