FDTD analysis of ELF radio wave propagation in the spherical Earth-ionosphere waveguide and its validation based on analytical solutions

. The FDTD model of electromagnetic wave propagation in the Earth-ionosphere cavity was developed under assumption of axisymmetric system, solving the reduced Maxwell’s equations in a 2D spherical coordinate system. The model was validated on different conductivity profiles for the electric and magnetic field components for various locations on Earth along the meridian. The characteristic electric and magnetic altitudes, the phase velocity and attenuation rate were calculated. 10 We compared the results of numerical and analytical calculations and found good agreement between them. The undertaken FDTD modeling enables us to analyze the Schumann resonances and the propagation of individual lightning discharges occurring at various distances from the receiver. The developed model is particularly useful when analyzing ELF measurements.


Introduction
The finite-difference time-domain (FDTD) is a numerical analysis technique based on time-dependent differential Maxwell's equations.It was originally developed for Cartesian coordinate system, but after elaboration of the code for spherical coordinates, it found applications in studies of ELF and VLF radio wave propagation in the Earth-ionosphere waveguide (Holland, 1983;Hayakawa and Otsuyama, 2002;Simpson and Taflove, 2002;Otsuyama et al., 2003;Yang and Pasko, 2005;Yu et al., 2012;Samimi et al., 2015).Similar analysis can be performed for Mars and other planets (Soriano et al., 2007;Navarro et al., 2007;Yang, et al., 2006;Navarro et al., 2008).
When a small part of the Earth-ionosphere cavity needs to be analyzed, a local volume can be divided into FDTD grid in 2D cylindrical (Cummer, 2000;Hu and Cummer, 2006;Marshall, 2012, Qin et al., 2019), or 3D Cartesian coordinate systems (Araki et al., 2018;Suzuki et al., 2016).It facilitates taking into account some complex inhomogeneities and ionospheric anisotropy in the analysis of ELF/VLF radio waves.
FDTD modeling in 3D Cartesian coordinates system was used for verification of Wait's and Cooray-Rubinstein analytical formulas, describing lightning-radiated electric and magnetic fields for a mixed propagation path (vertically stratified conductivity) and for the fractal rough ground surface (Zhang et al., 2012a;Zhang et al., 2012b;Li et al., 2013;Li et al., 2014).
The analysis of propagation of over a mountainous terrain was performed in 2D axial symmetric model (Li et al., 2016), which was further developed into FDTD model in 2D spherical axisymmetric coordinate system (Li et al., 2019).In such modeling the authors have investigated the effect of the Earth-ionosphere waveguide structure and medium parameters, including the effect of the ionospheric cold plasma characteristics, the effect of the Earth curvature, and the propagation effects over a mountainous terrain.
A full-wave finite element method was used to calculate electromagnetic fields in a horizontally stratified ionosphere treated as a magnetized plasma (Lehtinen and Inan, 2008).The authors have implemented the source currents with arbitrary vertical and horizontal distributions.The electromagnetic field was calculated both in the Earth-ionosphere waveguide and in the ionosphere as an upward propagating whistler mode wave.This model was used to simulate trans-ionospheric propagation of VLF electromagnetic waves from ground-based transmitters up to satellite altitudes (Lehtinen and Inan, 2009).
In a recent paper of Nickolaenko et al. (2021) the propagation of natural wide-band pulsed radio signals in the spherical Earth-ionosphere cavity is numerically simulated.Applying the realistic vertical profiles of ionospheric conductivity the authors have obtained ELF-VLF propagation parameters by the full wave solution in form of Riccati differential equation.
Using the day and night ionosphere models the correspondent parameters were calculated for different source-observer distances.
The influence of ionospheric disturbances on the Schumann resonances was analyzed using 3D FDTD model in (Navarro et al., 2008).The authors estimated the role of day-night asymmetry, polar non-uniformities associated to solar proton events, and X-ray bursts.
Few other numerical approaches were used for the estimation the parameters of Earth-ionosphere resonator, like the twodimensional telegraph equation (TDTE) in a two-dimensional spherical transmission line model of the Earth-ionosphere cavity (Kulak et al., 2003, Prácser et al., 2019, Bozóki et al., 2019), Transmission Line Matrix (TLM) numerical method in Cartesian and spherical coordinate systems (Morente et al., 2003;Toledo-Redondo et al., 2016).
Besides the Schumann resonances, another important aspect of ELF studies concerns propagation of ELF waves generated by individual lightning discharges, including Q-burst phenomena (Ogawa et al., 1967).Their waveforms observed at large distance from the source are significantly influenced by the dispersive properties of the Earth-ionosphere waveguide and the over-the-world propagation.Inverse solutions in such case should take into account full non-uniform solutions of the Maxwell equations.This is particularly important for lighting discharges that have a long continuing current phase or are associated with Transient Luminous Events (TLE) when they occur at large distance from the receiver.
The development of an FDTD model was motivated by our ELF systems: World ELF Radiolocation Array (WERA) (Kulak et al., 2014) and a new European ELF radiolocation system (EERS) (Mlynarczyk et al., 2018), which are operating in the Extremely Low Frequency (ELF) range, as well as for possible stations on Mars.
In this paper, we present an FDTD uniform model in 2D spherical coordinates and introduce new approach for validation the FDTD models by introducing computation of complex altitudes of the Earth-ionosphere waveguide, which allows us to compare numerical results with two-scale-height analytical solutions.We also infer the resonance frequencies of Earthionosphere waveguide for various conductivity profiles using the decomposition method (Kulak et al., 2006).

FDTD model
We have created the FDTD model following the ideas, which were originally proposed in (Holland, 1983) and developed in further studies (Hayakawa and Otsuyama, 2002;Otsuyama et al., 2003;Yang and Pasko, 2005;Navarro et al., 2007;Yang, et al., 2006;Navarro et al., 2008).We chose a spherical coordinate system to be able to study the wave's propagation in the Earth-ionosphere cavity and analyze the Schumann resonances.
Since in the present work we did not intended to study the azimuthal dependence of propagation parameters, we reduced a 3D system of Maxwell's equations to a 2D axisymmetric system.It allowed us to significantly decrease the required computational power and enabled longer simulation times, which leads to a better frequency resolution (up to df = 0.01 Hz).

Update equations
In the case of spherical system, assuming no dependence on  coordinate, the system of six Maxwell's equations in 3D spherical coordinates ( r ,  ,  ) can be reduced to 2D spherical system ( r ,  ) with tree equations for r E ,  E , and  H field components (Holland, 1983;Inan and Marshall, 2011) is the conductivity profile of Earth-ionosphere cavity, 0  and 0  are the vacuum permittivity and permeability respectively.
These equations were discretized using central-difference approximations to the space and time partial derivatives.The resulting finite-difference equations are solved in a leapfrog manner: the electric field vector components in a volume of space are solved at a given instant in time; then the magnetic field vector components in the same spatial volume are solved at the next instant in time and the process is repeated over and over again until the desired transient or steady-state electromagnetic field behaviour is fully evolved (Inan and Marshall, 2011).So the update equations for r E ,  E and  H are the following (Holland, 1983) where t  is the time step, r  and   are the sizes of grid cell in r and  coordinates respectively, R is the mean Earth's radius.Superscript n signifies that the quantities are to be evaluated at time , and i , j represent the point ( ) in the spherical grid.The half time steps indicate that the electric and magnetic fields are calculated alternately.
The update equation for r E cannot be applied for poles (where 0 =  and   = ) because of  sin in denominator.To solve this singularity problem the Holland's approach was used in Holland (1983).Namely, for the poles the integral form of the curl equation with a small contour around the poles was applied, which leads to new update equations for poles without singularities.After that, the update equation for  0 =  takes the form (Holland, 1983): And for N is the number of grid cells in  direction.

Mesh
The size of grid cell (in r and  directions) is defined by the desired range of frequencies, which are going to be analyzed.
For the present analysis we used , which lets us analyze frequencies up to 1 kHz.The ground was modeled as a perfect electric conductor.Also, the upper boundary max R for the model is a perfect conductor but it was placed at high enough altitude to make sure there is no reflection of the waves from this boundary.
The total simulation time is defined by the desired frequency resolution of FFT ( , where c is the speed of light and a is the Earth radius (Inan     and Marshall, 2011).In order to get sufficient frequency resolution (in our case it's 0.01 Hz, to be able to study the differences between different models) we have conducted the FDTD simulation up to 100 s.

Source
In order to analyze the propagation of electromagnetic waves originating from lightning discharges, one has to use the source model with similar characteristic to real lightning discharge.We used the time profile of lightning discharge proposed in (Kulak et al., 2010;Rakov, 2007), which has a flat spectrum in a wide frequency range.But taking into account the restrictions connected with FDTD computational grid, we had to modify the sourcethe highest frequency (or smallest wavelength) is defined by the cell size, and according to the generally accepted rule the smallest wavelength should be at least 10-20 times larger than the cell size.Therefore, the source has to be filtered with a loss-pass filter in order to remove frequency components that does not fit the size of computational cell.Also, it should be noted, that if the source contains the direct current (DC) in its spectrum, it can introduce artifacts, which are not physical (Li et al., 2013).Therefore, additional filtering has to be used for the lowest frequencies by a high-pass filter.
For such filtering we have used a 5 th order Butterworth lowpass filter with cutoff frequency of 300 Hz, and a 3th order Butterworth highpass filter with cutoff frequency of 2 Hz.The final spectrum of this modified source is almost flat in the range 5-100 Hz but decreases rapidly for 3 < f Hz and 1000 > f Hz.
We placed the source in the FDTD grid at the pole (  0 =  ) in radial direction in nodes with i = 0 -6, which corresponds to the source length L = 6 km.Also, for the implementation of the source into the FDTD grid we took into account the lightning stroke speed, and for our analysis we used v=10 8 m/s (Rakov, 2007).We assumed a cloud-to-ground (CG) discharge and implemented it by the time delay between adjacent nodes in the source.

Validation of the model
The main purpose of the current work is the analysis of resonance phenomenon in the Earth-ionosphere cavity.In order to validate the developed FDTD model, we analyzed several configurations of spherical waveguide with known analytical solutions and compared it with our FDTD results.We compared the resonance frequencies for them in relative and absolute units.

Lossless spherical waveguide
The first model for testing was a lossless spherical waveguide with perfect electric conductors at the ground and the upper boundary.The distance between the conductors was set to h = 74 km.The precise theoretical resonance frequencies for such waveguide for a given height were presented by the equation , where c is the speed of light, a is the Earth radius, n is the mode of resonance frequency (Bliokh et al., 1977).A comparison of analytical and numerical results for such model is presented in Table 1.

Two-layered spherical waveguide
The third model for testing was a two-layered model with a perfect conductor at the ground, a constant conductivity numerical results for this model is presented in Table 3.

Application of a realistic atmospheric conductivity profile
The next validation of our FDTD model was done for a realistic atmospheric conductivity profile (Kudintseva et al., 2016;Nickolaenko et al., 2016).We did not take into account the influence of Earth magnetic field.The resulting anisotropy of the conductivity has little influence in the analyzed frequency range (i.e., 0-100 Hz).As it was shown in Yu et al. (2012), the 185 anisotropy has a more significant influence at higher frequencies.Since the upper boundary of this profile has the conductivity much below 1 S/m, which is not enough to attenuate the ELF waves, reflections in FDTD grid would occur from the PEC at its upper boundary.To avoid reflections and let the waves attenuate gradually at high altitudes we extended this profile up to 200 km using the IRI model.We chose it in such a way that 195 the combined profile is smooth (Figure 2).The required IRI profile was found on January 22, 2006, at 6:00 UT for the location typical for mid-latitudes (49N, 23E).

Characteristic electric and magnetic altitudes
To be able to compare the numerical results with the analytical solutions, we have extracted the complex characteristic electric and magnetic altitudes from the FDTD model.The corresponding altitudes can be extracted using their definition (Mushtak and Williams, 2002;Kirillov, 1993): where N is number of radial nodes, i r E and i H  are the complex field values at the radial node i for a given frequency, and 0 r E , 0


H are the values of these fields at the ground level.The complex electric altitude can also be expressed using the conductivity profile by the analytical equation in a normalized form (Mushtak and Williams, 2002) , Assuming a similar dependence for complex magnetic altitude and taking into account that the characteristic conductivity value for magnetic altitude is (Greifinger and Greifinger, 1978) , ) ( 4 where the scale height , we can write a similar equation for the complex magnetic altitude .We compared the characteristic electric and magnetic altitudes calculated analytically from equations ( 11) and ( 12) with the FDTD results, calculated from equations ( 9) and ( 10).The results are presented in Figure 4.

Propagation parameters
We calculated the phase velocity ph V and the attenuation rate  using two different methodsnumerical and analytical.
Analytically the propagation parameters can be calculated using conductivity profile and electric and magnetic altitudes given by equations ( 11) and ( 12), through the following relationship (Kulak and Mlynarczyk, 2011) In case of numerical calculations those parameters were obtained by two approaches: a) similarly to analytical as described above but using electric and magnetic altitudes from FDTD model by equations ( 9), ( 10), and b) directly from the spectra of electric and magnetic field components.Assuming that in our coordinate system the wave propagates in the  direction the following relationship can be written in the frequency domain where the ratio of magnetic field complex amplitudes  H (f) are calculated for two probes "1" and "2" which are located at distances 1

 and 2
 from the source respectively.In this equation the propagation constant    i = + , where  is the attenuation rate (in units Np/m) and  is the phase constant, so the phase velocity is . For further convenient usage we convert the units of attenuation rate to dB/Mm.Similar relationship to (15) but for electric field components r E (f) can be written as well.
Figure 5.The attenuation rate and phase velocity that was calculated analytically and numerically.Using numerical FDTD calculations those parameters were obtained by two approaches: a) using electric and magnetic altitudes from FDTD by equations ( 9), ( 10), and b) directly from the spectra electric and magnetic following the equation ( 15), where we used probes located at 80 0 and 90 0 .
We should note that equation ( 15) for r E and  H spectra calculated in closed cavity gives non monotonic behavior of propagation parameters.This is caused by superposition of wave attenuation along meridian and different amplitudes of the Schumann resonances for different source-observer distance.One of possible solutions for removing the influence of Schumann resonances is to implement a PML along  direction.However, we solved this problem in a different way, transforming the Maxwell equations at angle 0 90   from spherical coordinates to plane equations, changing the behavior of equations from "close" to "open".In this case the waves are unable to propagate around the Earth and therefore the Schumann resonance does not occur.The calculated propagation parameters from all the methods listed above are presented in Figure 5.

Spectral decomposition method
As an additional validation we applied the spectral decomposition technique to Schumann resonance power spectra.Figure 6 presents the resonance frequencies obtained from the spectra measured at different distance from the source, and the decomposed frequencies that do not depend on the source-observer distance.Following the decomposition algorithm (Kulak et al., 2006;Dyrda et al., 2014) the spectrum is approximated by the function where W(f) is the signal power spectrum, s is the white noise component, /  is the color noise term,   is the power parameter of the k-th resonance peak,   is the peak asymmetry parameter,   is the resonance frequency, and   is the resonant mode half-width parameter.Fitting this function to the FDTD spectra allows us to extract the resonance frequencies   of the cavity, which are not equal to Schumann resonance frequencies.The Schumann resonance frequencies obtained from the spectra depend on the distance from the source and represent the superposition of standing and traveling waves (Kulak et al., 2006).To analyze the standing waves separately and reveal the properties of the resonant cavity we are using the spectral decomposition method described in (Kulak et al., 2006;Dyrda et al., 2015).After applying this method, the resonance frequencies become independent of the source-observer distance (see details in Kulak et al. (2006) and Dyrda et al. (2015)).
The decomposition method shows that the solutions for the electric field are symmetric at   (Kulak et al., 2006).

Resonance frequencies of the Earth-ionosphere cavity
The resonance frequencies of the Earth-ionosphere cavity for a given conductivity profile can be calculated by different approaches and the consistency between them can be considered as an additional validation of the model.We have considered the followings ways to calculate the resonant frequencies: 1).Analytically using the conductivity profile.The resonance frequencies in this approach are calculated by solving the following equation (Mushtak and Williams, 2002;Galejs, 1972) where e m h h S / = 2 depending on the characteristic complex altitudes, which are discussed in Section 4.1; 2).Numerically from FDTD model using equation ( 17).We denote these resonant frequencies by n h f ; 3).Using FDTD spectra for 180 r E (see Section 4.3 for details).We denote these frequencies by n r E f ; 4).Using spectral decomposition method (see Section 4.3 for description).We denote these resonant frequencies by n d f .
The obtained resonance frequencies are presented in Table 4.These frequencies we compared with the analytical results and the absolute and relative differences are presented in the corresponding columns.Table 4 Resonance frequencies for the conductivity profile shown in Figure 2 calculated by different approaches described in Section 4, and the relative error obtained by comparison with the analytically obtained resonance frequency.

Discussion
In this study, we used new methods for validation of numerical simulation: Mode n c f [Hz]   n h f [Hz]   n r E f [Hz]   n d f [Hz]    1).We compared the complex electric and magnetic altitudes of the Earth-ionosphere waveguide, referring to twodimensional formalism of electromagnetic wave propagation in the Earth-ionosphere cavity.The two complex altitudes were calculated numerically directly from their definitions ( 9) and ( 10), making use of radial field solutions E(r) and H(r).These altitudes were directly compared with the altitudes obtained with analytical formulas ( 11) and ( 12) for the same conductivity profile of the atmosphere.This allowed us to fully validate the simulation results.
2).We determined the resonance frequencies of the cavity using the decomposition of power spectra described by equation ( 16).The resonance frequencies cannot be directly determined from the spectra obtained using FDTD, because the field generated by the source in each point of the cavity is a superposition of travelling waves propagating directly from the source and the resonance field resulting from the interference of waves propagating around the world.The resulting Schumann resonance frequencies depend on the source-observer distance (Kulak et al., 2006).Close to the source, where the amplitudes of the travelling waves are significant, the resonance frequencies differ by several percent from the Schumann resonance frequencies obtained by FDTD.With the use of the decomposition method we determined the intrinsic resonance frequencies of the cavity, which are the same at each location, independently from the distance to the source.This enabled us to compare the resonance frequencies obtained from the numerical simulation with the frequencies determined directly from analytical formula ( 9).This method should be recommended as a reference for validation of numerical models.
3).We determined the propagation parameters of the Earth-ionosphere waveguide using the FDTD results in two different points of the great circle.Equation ( 15) that we used allows us to determine the phase velocity and the attenuation rate of the Earth-ionosphere cavity.We compared them with the results obtained from analytical formulas ( 13) and ( 14).

Summary and conclusions
In this paper, we analyzed the solutions of Maxwell's equations obtained by the FDTD method for an axisymmetric uniform Earth-ionosphere cavity.We analyzed the propagation of radio waves generated by a short current impulse.
We have constructed an FDTD model in axisymmetric spherical coordinate system with the source implemented at the pole.
We took into account the finite speed of lightning discharge and implemented the time delay between adjacent nodes in the source.
We validated the model thoroughly, comparing the resonance frequencies, propagation parameters and electric and magnetic characteristic altitudes.Since the conductivity profile of the atmosphere has a significant influence on radio wave propagation and resonance frequencies, we validated our model for various conductivity profiles.
We paid a close attention to the verification of accuracy of the FDTD computations and used a new approach for that purpose.It is based on 2D formalism of wave propagation in the Earth-ionosphere cavity and it allowed us to compare the numerical and the analytical solutions.In this approach, the propagation parameters of the Earth-ionosphere waveguide are defined using the electric and magnetic altitudes.These altitudes can be calculated directly, when the vertical conductivity profile of the atmosphere is known.
The obtained analytical solutions were used as reference and compared with the numerical results.First, we compared the solutions for three cases: a perfect cavity, a cavity formed by a perfect ground and a homogenous conducing layer, and a cavity formed by a perfect ground and two-layered upper boundary of the waveguide.As a measure of error between the models, we took the difference between the first three resonance frequencies.The analytical and numerical solutions were in agreement (Tables 1, 2, 3).Next, we compared the results for a continuous conductivity profile.We built a realistic conductivity model, which lower part was based on a recently proposed conductivity profile and the upper part was based on an IRI model.We obtained a good agreement between the resonance frequencies of the cavity and the observed Schumann resonance frequencies (Table 4).
Using our FDTD model, we also calculated the spectral dependence of the phase velocity and the attenuation rate.We showed that the analytical and numerical models are in agreement.
The presented model can be used for studying the propagation of ELF electromagnetic waves generated by lighting discharges of various types with the round-the-world propagation taken into account.
above 110 km.The analytical results for such model can be obtained fromKulak et al. (2013), which deal with multi-layer waveguides.A comparison of analytical and

Figure 1 .
Figure 1.The FDTD spectrum for spherical waveguide with lossless perfect cavity (a), one conducting layer (b) and two layers (c) for electric component at 180

Figure 2 .
Figure 2. A realistic conductivity profile up to 200 km obtained by combining the profile from Kudintseva et al. (2016) (red part) and an IRI profile (blue part).

Figure 3 .
Figure 3.The results of the FDTD analysis are shown, that were obtained using the realistic conductivity profile described in Figure 2. The solutions for electric component Er (upper panels) and magnetic component Hp (bottom panels) are shown in time domain (left panels) and in frequency domain (right panels).Each plot contains the results for two probe position (60 0 and 120 0 ).

Figure 4 .
Figure 4. Real and imaginary parts of characteristic electric and magnetic altitudes obtained analytically (denoted by "A") and from the FDTD model (denoted by "N").
other and only the standing waves remain.Therefore, the resonant peaks measured from the spectrum at  180 =  for the electric field component ( 180 r E ) represent the resonance frequency of the cavity and they are in agreement with the decomposed frequencies.

Figure 6 .
Figure 6.The Schumann resonance frequencies of first three modes obtained at different source-observer distances and the resonance frequencies of the cavity obtained by the decomposition method(Kulak et al., 2006).

Table 1
Comparison of resonance frequencies obtained analytically and numerically for a lossless spherical waveguide.The analytical solutions for this model were obtained byKulak and Mlynarczyk (2013).A comparison of analytical and numerical results for such model is presented in Table2.

Table 3
Comparison of resonance frequencies obtained analytically and numerically for a two-layered spherical waveguide.