## 1 Introduction  
Very-low-frequency (VLF) electromagnetic waves refer to electromagnetic waves with a frequency range of 3 kHz to 30 kHz. VLF electromagnetic waves have low propagation loss and slow attenuation, and the signal strength and phase are relatively stable. They are often used in submarine communication and geophysical research [1, 2, 3, 4, 5]. In nature, VLF electromagnetic waves can be generated by lightning discharges. The experimental and research results have shown that amplitude-modulated high-frequency electromagnetic waves can be excited to generate ELF/VLF electromagnetic waves by heating the lower ionosphere artificially [5]. The ionosphere is closely related to radio communication, and electromagnetic waves transmitted in it are often affected by changes in the ionosphere [6, 7, 8, 9]. Therefore, it is necessary to study the propagation characteristics of VLF electromagnetic waves in the ionized plasma of the ionosphere for applications such as VLFcommunication, lightning distribution, and ionospheric heating. The electron density in the ionosphere changes rapidly in the vertical direction and changes slowly in the horizontal direction [4,10,11]. Some scholars have conducted research on the related issues. There are three main calculation methods for electromagnetic-wave propagation in the very-low-frequency field. One is the numerical method, the second is the full-wave method (FWM), and the third is the finite-difference time domain (FDTD) method. The authors of ref. [11] numerically studied the amplitude and spectrum of strong electromagnetic pulses penetrating the ionosphere upwards, providing a reference for the design of electromagnetic pulse protection for low-orbit satellites and other aircraft. When simulating the propagation characteristics of electromagnetic waves in the ionosphere, it is usually necessary to divide the ionosphere into many thin layers. For complex anisotropic media and strong instability in the ionosphere, directly using numerical calculation methods results in insufficient accuracy. The authors of ref. [12] studied the amplitude and spectrum of strong electromagnetic pulses penetrating the ionosphere upwards, providing reference for the design of electromagnetic pulse protection for low-orbit satellites and other aircraft. Reference [13] introduces the current frequency-band distribution of satellite communication and compares it, providing a link budget for S-band satellite communication. The authors of reference [14] used the propagation matrix method to study low-frequency seismic electromagnetic radiation signals. The authors of reference [15] simulated the propagation of electromagnetic pulses in the ionosphere using the FDTD method. The authors of reference [16] applied Z-transform to obtain the iterative relationship between the electric displacement vector \(D\) and the electric-field strength \(E\) in FDTD derivation, and analyzed the propagation characteristics of high-frequency radar signals in the ionosphere. When simulating the propagation characteristics of electromagnetic waves in the ionosphere, it is usually necessary to divide the ionosphere into many thin layers. With the continuous increase in ionospheric calculation altitude, the time step of FDTD model is limited by the minimum grid size in space, which leads to high computational complexity and memory-space occupation, resulting in error accumulation and inaccurate calculation results [16, 17, 18]. Therefore, how to quickly and accurately analyze the propagation characteristics of very-low-frequency electromagnetic waves in the ionosphere has become a difficult problem in the fields of computational electromagnetics and radio-wave propagation in recent years [19]. In this paper, a hybrid explicit-implicit time-domain finite difference (HIE-FDTD) method is proposed in the ionosphere.