Full-wave method for the analysis of the radiation characteristics of a VLF source in the atmosphere

Journal Pre-proofs Full-wave method for the analysis of the radiation characteristics of a VLF source in the atmosphere Yin Weike, Wei Bing, Zhang Shitian PII: DOI: Reference:

S2211-3797(19)31470-6 https://doi.org/10.1016/j.rinp.2019.102682 RINP 102682

To appear in:

Results in Physics

Received Date: Revised Date: Accepted Date:

9 May 2019 28 August 2019 16 September 2019

Please cite this article as: Weike, Y., Bing, W., Shitian, Z., Full-wave method for the analysis of the radiation characteristics of a VLF source in the atmosphere, Results in Physics (2019), doi: https://doi.org/10.1016/j.rinp. 2019.102682

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2019 The Authors. Published by Elsevier B.V.

Full-wave method for the analysis of the radiation characteristics of a VLF source in the atmosphere Yin Weike1,2, Wei Bing1,2 , Zhang Shitian3 1. School of Physics and Optoelectronic Engineering, Xidian Univ., Xi’an 710071, China; 2. Collaborative Innovation Center of Information Sensing and Understanding, Xidian Univ., Xi’an 710071, China; 3. China Research Institute of Radiowave Propagation, Qingdao 266107, China

Abstract: The basis of the full-wave method is the propagation matrix method and a plane wave expansion, which is generally used for the analysis of radiation problems in layered media. The radiation field in the spatial domain is represented as an integral in the wavenumber domain. By integrating the vertical wavenumber, the radiation can be represented as a superposition of plane waves with different directions. If the atmosphere and ionosphere are regarded as layered media, the vertical wavenumber integral contains first-order singularities when the radiation source is located in the ionosphere, which are relatively easy to handle. When the radiation source is located in the air, the vertical wavenumber integral contains second-order singularities, which are more complicated. Therefore, there are few full-wave studies on radiation characteristics in the air and ionosphere when a very-low-frequency (VLF) source is located in the air. In this paper, air is considered as an isotropic medium with a microconductivity, the vertical wave number integral that has second-order singularities is solved, and the radiation fields in the air and the ionosphere are calculated. The numerical results show the correctness and effectiveness of the proposed algorithm.

Keywords: electric dipole; VLF; radiation field; full-wave method; uniform anisotropy 1 Introduction Very-low-frequency radio signals are mainly used in the fields of long-range navigation, ultralong-range strategic reliable communication, time service, low ionosphere detection and earthquake prediction. A VLF signal propagates mainly in the waveguide between the Earth surface and the lower edge of the low ionosphere. The propagation properties of a VLF signal are affected by changes in the ionosphere and the dielectric properties of the Earth's surface[1]. Sky-wave interference caused by ionosphere disturbances is the main factor for the decline in the accuracy of VLF navigation. For decades, researchers have continued to study the propagation of VLF waves. In 1975, Nagano et al. presented a propagation matrix method to calculate the field distribution of plane waves in horizontally layered anisotropic media [2]. In 1994, Yagitani, Isamu et al. calculated the up and down

VLF radiation efficiencies of low ionosphere heating in the polar region using the full-wave method [3, 4]. In 2008, Lehtinen et al. used the propagation matrix method combined with a plane wave expansion to study radiation at 700 km from the ground with a radiation source located in the D layer (stratification in the ionosphere, approximately 50 km-100 km above the ground) [5]. Wang Feng et al. used the full-wave method to calculate the wave field of a VLF/ELF wave in the horizontal layered ionosphere radiated by an electric dipole moment formed by an artificially modulated ionosphere[68]. In 2004, Bortnik et al. used the propagation matrix method to analyse the problem of a plane wave propagating from the Earth surface to low Earth orbit [9]. In 2013, Lehtinen et al. used the full-wave method to calculate the effect of the polarization, incidence angle, and azimuth on the transmission characteristics of VLF plane waves as they traverse the ionosphere [10]. The above methods are based on the propagation matrix method given by Nagano for analysing plane wave propagation in an anisotropic ionosphere. A Green’s function in the spatial domain can be obtained by a Fourier transform of the Green’s function in the wavenumber domain. The radiation field at the observation point can be expressed as an integral in the wavenumber domain. By integrating the vertical wavenumber, the radiation can be represented as a superposition of plane waves with different directions. The vertical wavenumber integral contains first-order singularities when the radiation source is located in the ionosphere, which are relatively easy to handle. The vertical wavenumber integral contains second-order singularities when the radiation source is located in the air, which are relatively complex. Therefore, there are few full-wave studies of radiation characteristics in the air and ionosphere when an electric dipole is in the air. In this paper, we use the full-wave method to calculate radiation in the Earth–ionosphere waveguide when an electric dipole is located in the air. To facilitate the calculation of the vertical wavenumber integral, we consider the air as an isotropic medium with a micro-conductivity [11]. A method for solving the vertical wavenumber integral with second-order singularities is given. The radiation in the air and the ionosphere is calculated. The numerical results show the correctness and effectiveness of the proposed algorithm. 2 Physical model and calculation method Physical model The ionosphere is regarded as an anisotropic medium in the background of the Earth's magnetic field, and the ground is regarded as an ideal conductor. The layered model of the ionosphere used in the present work has the same ionosphere parameters as in the work of Nagano [2]. The groundionosphere waveguide model is shown in Fig. 1. Z1 is the interface of the air and ionosphere, Z 0 is the layer where the electric dipole source is located, and Z 1 is the Earth surface. An electric dipole is located along the y-axis, as shown in Fig. 1.

z

... ...

Zm-1 Zm-2

ionosphere

Zj Z2 Z1

air dipole ground



Z0 y Z-1

Figure 1 The ground-ionosphere stratification model Calculation method Plane wave expansion of the radiation of an electric dipole located in air For a plane electromagnetic wave with an angular frequency of  , the electric field satisfies the wave equation [12, 13]:     E  k 02   E   j0 J ,

(2)

where k0 denotes the wavenumber in free space,  0 is the permittivity of free space, and  0 is the magnetic permeability of free space. Assuming that air has a micro-conductivity and   is the socalled dielectric tensor of air, n2     0 0 

where n 2  1  i

0 n2 0

0  0, n 2 

 ,  is the conductivity of the air, i denotes the imaginary unit, and J denotes a  0

current source. The electric field in (2) can be integrated by the Green's function G  r, r   and the current source in the region where the source is located: E  r    j0  G  r, r J  r  dr  ,

(3)

where r and r  are vectors for the observation point and source point, respectively. and dr is equal to dx dy dz  , G  r, r   

1

3  2  

Ω1 exp k 02

 jk r  rdk ,

(4)

Ω 

1 kk   k 2 I +k02   k 02

 q2  s y 2  n 2   s y sx  qsx 

sx s y  q 2  sx 2  n 2 qs y

， sx q  sy q   sx 2  sy 2  n 2 

where k is the wavevector, k  k , k x  k0 sx , k y  k0 s y , q  n 2  sx2  s 2y , and  I  is a unit dyadic tensor. When the electric dipole is located r =  0, 0, z0  , the integral in the source region in (2) then becomes

 J r exp  jkr dr=Il exp  jk qz  m , 0

0

where m is the direction of the source dipole, Il is the dipole moment of the electric dipole, and then (3) becomes E r  

j0 k0 Il

 2 

3

 Ω

1





m exp ik0  sx x  s y y  q  z  z0   dsx ds y dq .



(5)

By integrating q , we can represent the radiation of (4) as a superposition of a series of plane waves with fixed sx and s y . The field of the plane wave is expressed as E  sx , s y , z  

j0 k0 Il Ω1 exp  ik0 q  z  z0   mdq ,  2

Ω1  adj Ω / det Ω ,

(6) (7)

where adj Ω denotes the adjoint of the matrix Ω and det Ω denotes the determinant of Ω : adj  Ω  Ω0    q  qu  q  qd  ,

where    n 2  sx2   sx s y sx q   ,   n 2  s 2y  sy q Ω0    s x s y   2 2   s q s q  n  q   x y  

and det Ω  a  q  q u   q  q d  . 2

2

a is coefficient of q 4 , qu and qd are eigenvalues of Ω denoting the up-going and down-going waves, respectively, qu  n 2  sx2  s 2y , qd   n 2  sx2  s 2y . Simplifying (7), we obtain

Ω1 

Ω0  . a  q  qu  q  qd 

(8)

Obviously, the wavenumber domain solution of the electric field can be obtained by completing the integral of (6). However, when the radiation source is located in the air, if the air is regarded as a lossless medium, the singularities in (8) are along the real axis; when sx2  s 2y  n 2 , and qu  qd  0 , the

integral (6) is uncertain. If the air is considered as an isotropic medium with a micro-conductivity, there are no singularities along the real axis in (8), and then (6) can be calculated using the residue theorem. Completing the integral calculation, the electric fields of the up-going and down-going waves are expressed as Eup  sx , s y , z   0 k0 Il

Ω0  |q q m  exp  ik0 qu  z  z0   , a  qu  qd 

Edown  sx , s y , z   0 k0 Il

u

Ω0  |q q m  exp  ik0 qd  z  z0   . a  qd  qu  d

The electric field in the spatial domain can be calculated by the integration of a large number of elementary plane waves： E r  

E  sx , s y , z  exp  ik0  sx x  s y y   dsx ds y . 2   2   1

(9)

When the amplitude of E  sx , s y , z  decreases to -30 dB of the maximum value in the wavenumber domain, the upper and lower limits of the integral are selected, and (9) can be written in a discrete form that can facilitate numerical calculations [3]. Propagation of VLF waves in the low ionosphere The dielectric tensor of the lower ionosphere has the same form in each layer and varies with





height: ε j   0 I   M  j . I is the unit matrix, and  M  j is the plasma dielectric tensor for layer j [14]  U 2  l 2Y 2 X  inYU  lmY 2 M j = 2 2  U U  Y   2  imYU  ln Y X

inYU  lmY 2 U 2  m 2Y 2 ilYU  mnY 2

imYU  ln Y 2   ilYU  mnY 2  , U 2  n 2Y 2 

(10)

 p2   , Y  h , Z  e , U  1  iZ ,  p is the plasma frequency, h is the magnetic cyclotron 2   

frequency,  is the angular frequency of the wave, and  e is the electronic collision frequency. l , m ,and n are the projections of the magnetic field directions on the x-axis, y-axis, and z-axis,

respectively. The plasma dielectric tensor for each layer is determined by the electron density, collision frequency, and magnitude and direction of the geomagnetic field in this layer. Suppressing

the

time

variation

eit ,

the

incident

wave

can

be

expressed

as

exp it  ik0  sx x  s y y  qz   , and electromagnetic waves in a medium satisfy the following Maxwell

equations[15-17]:   H  i 0  I   M  E   E  i0 H

Eliminating E z and H z from (11), we can obtain

.

(11)

Ex Ey de  ik0 Te , e  , Z dz 0H x Z0 H y

(12)

where T is a 4  4 matrix [14, 18, 19]: 

sx M zx 1  M zz s y M zx

T

1  M zz 1  M xx 

sx s y

1  M zz

1  M zz



1  M zz M yz M zx

sx M zy

 M yx  s x s y

s y M zy 1  M zz

1  s  M yy  2 x

M xz M zx 2  sy 1  M zz

M xz M zy 1  M zz

s2y

1

M yz M zy

sy

1  M zz

 M xy  sx s y

sx s y

1  M zz 

1  M zz 

 M yz

sx M yz

1  M zz 

sy

1  sx2  M zz 1  M zz

M xz 1  M zz

,

(13)

1  M zz 

sx M xz 1  M zz

and Z 0  0 /  0 is the wave impedance in free space. From (12)，we can find that the field at z  z j can be depicted as e  z j   B jA j ,

(14)

where B j are the eigenvectors of the matrix T , Aj  is the amplitude of the  th characteristic wave in layer j,   1, 2 are up-going waves, and   3, 4 are down-going waves. We can obtain e  z j 1   B j Δ j A j , e Δj 

 

 jk0  j 1

0

0 e

 

 jk0  j 2

0

0

0

0

e

(15)

0

0

0

0

 

 jk0  j 3

0

,

0 e

 

 jk0  j 4

 j   q j   z j 1  z j  ,

(16)

where q j  is the  th eigenvalue of the characteristic wave in layer j. Eliminating A j from(14) and (15), we have e  z j 1   K j e  z j  ,

(17)

where K j  B j Δ j B j 1 , B j 1 is the inverse of B j , and the medium above zm 1 is considered to be uniform. The wave at the boundary of zm 1 is x (1) 0 e1  zm 1   B m , 0 0

0 x (2) e 2  zm 1   B m . 0 0

(18)

The amplitudes of the two up-going waves x 1 and x  2  remain unknown. On the basis of (17), the field at z  z1 can be represented by the unfolding of the field at z  zm 1 : m 1

m 1

s 2

s 2

e  z1   x (1) K s e1  zm 1   x (2) K s e2  zm 1  .

(19)

The field at z  z0 can be represented by the field at z  z1 [5] e  z0   Δ 0  1e  z1  ,  nz   nz    cos   sin    cos  n n  nz   nz  cos    cos    sin    sin  , n n nz cos  n sin   n z cos  n sin  n cos  nz sin   n z sin  n cos   sin 

where   a tan  s y / sx  , and the field at z  z0 can be represented by the field at z  z0 : e  z0   e  z0   e down  e up .

(20)

The field at z  z1 can be represented by the field at z  z0 : e  z1   Δ 1 1e  z0  .

(21)

In terms of the VLF band, the Earth surface can be regarded as a perfect electric conductor. e  z1   0, 0, x (3), x(4) .

(22)

The unknown coefficients x 1 , x  2  , x  3 ,and x  4  can be determined by applying the boundary conditions at z  z1 using (19)-(21); then, e(z m1 ) can be determined, and (17) can be used to obtain the horizontal electromagnetic fields in the remaining layers. 3 Simulation results The radiation field of an electric dipole in infinite air with a micro-conductivity is calculated using the full-wave method. The results obtained using the full-wave method are compared with the far-field analytical solution of an electric dipole in vacuum. The electric field along the x-axis of the plane 45 km from the electric dipole plane is shown in Fig. 2. The operating frequency is 10 kHz, the electric dipole direction is along the y-axis, and the electric dipole moment is Il  8.4  103 A  m .

1.5

real(Ey)/mv

1.0

Full-wave method =4x10-10 S/m Analytical solution =0 S/m

0.5 0.0 -0.5 -1.0 -1.5 -500-400-300-200-100 0 100 200 300 400 500 x/km

Figure 2 Simulation result in infinite air The ground-ionosphere model is reduced to a ground-air half-space model. The full-wave method is used to calculate the electric field along the x-axis of the plane 10 km above the ground. Compared with the analytical solution, as shown in Fig. 3, the electric dipole is oriented along the y-axis, 45 km above the ground, and the operating frequency is 10 kHz. The assumption of considering the Earth's surface as a perfect planar conductor is examined. 2.5 2.0

Full-wave method =4x10-10 S/m Analytical solution =0 S/m

real(Ey)/mv

1.5 1.0 0.5 0.0 -0.5 -1.0 -1.5 -2.0 -2.5 -500-400-300-200-100 0 100 200 300 400 500 x/km

Figure 3 Simulation result in the air-ground half space In the ground-ionosphere waveguide model and air-ground half-space model, we calculate the radiation field of the VLF horizontal electric dipole. The electric field along the x-axis of the plane 1 km above the ground is shown in Fig. 4. Fig. 4 shows that the ionosphere plays a major role in VLF propagation at long distances. Snapshots of the electric field in the xoz plane in both models are shown in Figs. 5 and 6. The electric dipole is oriented along the y-axis, 1 km above the ground. Fig. 5 shows that the wavelength becomes shorter in the ionosphere than in the air. In the ionosphere, the energy of the electromagnetic field mainly follows the vertical direction (the direction of the Earth's magnetic field), and a small part of the electromagnetic field can penetrate the ionosphere.

0.10 0.05

Half space =4x10-10 S/m Ground-ionosphere waveguide =4x10-10 S/m

real(Ey)/mv

0.00 -0.05 -0.10 -0.15 -0.20 -0.25 -0.30 -300

-200

-100

0

100

200

300

x/km

Figure 4 Electric field in the ground-ionosphere waveguide and half space

Figure 5 Electric field in the xoz plane in the ground-ionosphere waveguide

Figure 6 Electric field in the xoz plane in the air-ground half space 4 Conclusion In this paper, a propagation model of the ground, the air, and a layered anisotropic ionospheric space is established. When an electric dipole is located in the air, the vertical wavenumber integral that has a second-order singularity is solved. The radiation in the air and the ionosphere is calculated. The calculation results of the full-wave method are compared with the analytical solutions in infinite

air and the air-ground half space. The numerical results show the correctness and effectiveness of the proposed algorithm. References：

[1] W. M. McRae and N. R. Thomson, "Solar flare induced ionospheric D-region enhancements from VLF phase and amplitude observations," Journal of Atmospheric and Solar-Terrestrial Physics, vol. 66, pp. 77-87, 2004/01/01/ 2004. [2] I. Nagano, M. Mambo, and G. Hutatsuishi, "Numerical calculation of electromagnetic waves in an anisotropic multilayered medium," Radio Science, vol. 10, pp. 611-617, 1975. [3] N. Isamu, M. Kazutoshi, Y. Satoshi, K. Iwane, O. Toshimi, H. Kozo, et al., "Full wave calculation method of VLF wave radiated from a dipole antenna in the ionosphere-analysis of joint experiment by HIPAS and akebono satellite," Electronics and Communications in Japan (Part I: Communications), vol. 77, pp. 59-71, 1994. [4] S. Yagitani, I. Nagano, K. Miyamura, and I. Kimura, "Full wave calculation of ELF/VLF propagation from a dipole source located in the lower ionosphere," Radio Science, vol. 29, pp. 39-54, 1994. [5] L. N. G. and I. U. S., "Radiation of ELF/VLF waves by harmonically varying currents into a stratified ionosphere with application to radiation by a modulated electrojet," Journal of Geophysical Research: Space Physics, vol. 113, 2008. [6] Z. Shu-Fan, L. Li, and Z. Xue-Min, "Trans-ionospheric VLF wave power absorption of terrestrial VLF signal," Chinese Journal of Geophysics (in Chinese), vol. 60, pp. 3004-3014, 2017. [7] F. Wang, Z.-Y. Zhao, S.-S. Chang, G.-X. Zhao, and H.-Y. Qing, "Radiation of ELF waves by ionospheric artificial modulation into a stratified ionosphere," Diqiu Wuli Xuebao, vol. 55, pp. 2167-2176, 2012. [8] C. Shan-Shan, Z. Zheng-Yu, and N. I. Bin-Bin, "Propagation and application of VLF waves generated by HF heating of ionosphere in the magnetosphere," Chinese Journal of Geophysics (in Chinese), vol. 56, pp. 3960-3968, 2013. [9] B. T. Bortnik J, "Full Wave Calculation of the Source Characteristics of Seismogenic Electromagnetic Signals as Observed at LEO Satellite Altitudes," AGU Fall Meeting Abstracts, pp. T51B-0453, 2004. [10] K. L. Graf, N. G. Lehtinen, M. Spasojevic, M. B. Cohen, R. A. Marshall, and U. S. Inan, "Analysis of experimentally validated trans-ionospheric attenuation estimates of VLF signals," Journal of Geophysical Research: Space Physics, vol. 118, pp. 2708-2720, 2013. [11] K. V. Kumar, S. Sampath, S. M. Das, and V. S. Kumar, "Atmospheric electrical conductivity variations over different environments," Geophysical Journal International, vol. 122, pp. 89-96, 1995. [12] J. A. Kong, "Theory of electromagnetic waves," New York, Wiley-Interscience, 1975. 348 p., 1975. [13] L. Huang and Y. Chen, "Localized excitations and interactional solutions for the reduced Maxwell-Bloch equations," Communications in Nonlinear Science and Numerical Simulation, vol. 67, pp. 237-252, 2019/02/01/ 2019. [14] K. G. Budden, The Propagation of Radio Waves: The Theory of Radio Waves of Low Power in the Ionosphere and Magnetosphere: Cambridge University Press, 1988.

[15] L. Huang and Y. Chen, "Nonlocal symmetry and similarity reductions for the Drinfeld–Sokolov–Satsuma–Hirota system," Applied Mathematics Letters, vol. 64, pp. 177-184, 2017/02/01/ 2017. [16] L. Huang, Y. Yue, and Y. Chen, "Localized waves and interaction solutions to a (3+1)-dimensional generalized KP equation," Computers & Mathematics with Applications, vol. 76, pp. 831-844, 2018/08/15/ 2018. [17] Y. Yue, L. Huang, and Y. Chen, "Localized waves and interaction solutions to an extended (3+1)-dimensional Jimbo–Miwa equation," Applied Mathematics Letters, vol. 89, pp. 70-77, 2019/03/01/ 2019. [18] P. Bespalov and V. Mizonova, "Propagation of a whistler wave incident from above on the lower nighttime ionosphere," in Annales Geophysicae, 2017, p. 671. [19] H. Zheng and D. Ge, "Electromagnetic wave reflection and transmission of anisotropic layered media by generalized propagation matrix method," Acta Phys. Sin., vol. 49, p. 1702, 2000-09-20 2000.

This paper established the propagation model of ground, air and layered anisotropic ionospheric space. Considering the electric dipole located in the air, the vertical wave number integral that has second-order singularity was solved and calculated the radiation in the air and the ionosphere. The dipole radiation in the half-space and in the Earth's atmosphere is compared. The graphs of the tangential component of the electric field of a short dipole at different distances from it are given. The field of the horizontal electric dipole in the layer of the earth's atmosphere with the ionosphere and in the half-space is shown and compared.