Abstract
The hypersonic flow field around a reentry vehicle has a significant influence on the ground-vehicle communication as well as on the detection and recognition of the reentry vehicle. Backward scattering characteristics of a reentry vehicle enveloped by a hypersonic flow field are analyzed using a high-order auxiliary differential equation finite difference time-domain (ADE-FDTD) method in this paper. Flow field parameters, including electron density, neutral particle density, and temperature, are obtained by solving the Navier-Stokes (NS) equations numerically. According to the flow field parameters, distributions of the plasma frequency and the collision frequency are then derived. Based on a validity of the physical model and the high-order ADE-FDTD method, backward radar cross sections (RCSs) of a perfect electrical conductor (PEC) sphere enveloped by a hypersonic flow field under different Mach numbers, heights, and incident angles of the electromagnetic (EM) wave are then investigated. Numerical results show that the incident angle of the EM wave exerts noticeable effects on the backward RCS, which is due to an inhomogeneous distribution of the plasma. The flight height and Mach number have significant influences on the distribution of the plasma that they play an important role in the variation of the RCSs. The results presented in this paper provide useful reference data for practical tests in ballistic range or in the high-frequency plasma wind tunnels, where a sphere target is usually used due to its simple shape.
1. Introduction
The development of hypersonic reentry vehicles in the near space, such as shuttles, rockets, and missiles, attracts increasing attention in recent years due to their significant strategic position. When a hypersonic vehicle reenters the atmosphere, severe friction between the aircraft and the ambient air results in a dramatic increase in the local temperature and pressure, which leads to intensified air ionization, and a layer of plasma flow is generated around the reentry vehicle [1]. The plasma flow, which consists of numerous free electrons, ions, and neutral particles, is nonmagnetic, low-temperature, weakly ionized, and quasi-neutral (i.e., electricity neutral on macroscopic-length scales) [2, 3]. It greatly affects the transmission of electromagnetic (EM) waves that it has a very strong interference on the communication between the reentry aircraft and the ground station; it even leads to blackout in certain conditions [4, 5]. Moreover, the reflection, refraction, and absorption of EM waves caused by the plasma sheath or wake change the radar scattering characteristics of the hypersonic vehicle [6–9], which has a significant influence on the detection, tracking, and recognition of the target.
In the past decades, lots of studies have been carried out on the analysis of the formation of plasma flow around a hypersonic vehicle and their interactions with EM waves. Early in the 1960s, a series of hypersonic experiments were performed in the Radio Attenuation Measurements (RAM) program to gain a thorough knowledge of the plasma sheath properties and its effects on the communication systems of the reentry vehicle. Some techniques to alleviate communications blackout were developed [5]. Due to the expensive costs and poor repeatability of the reentry flight experiment, laboratory simulations of hypersonic plasma flows [10], for example, using ballistic range [11] or plasma wind tunnel [12], were developed that the field properties of plasma flow as well as their EM scattering characteristics can be determined [13]. So far, most of the efforts were devoted to the analysis of transmission, reflection, and attenuation of EM waves in the plasma sheath in order to alleviate the interruption of plasma on ground-vehicle communication and to find out an alternative approach to bridge the linkage during the blackout time [14–17]. Limited attention was paid to the analysis of the scattering properties of the vehicle enveloped by a hypersonic plasma flow, which plays an essential role in the detection, tracking, and identification of hypersonic targets.
Due to the lack of experimental data on the plasma flow and the restriction of computing ability, in most of the previous literature, plasma flows with certain electron density profiles were assumed. And the analysis of the scattering by the plasma flow was implemented using approximated approaches, which permit to predict the plasma effects qualitatively. Backscattering cross sections and differential cross sections were calculated for a conducting cylinder coated with an anisotropic and inhomogeneous plasma sheath by Rusch and Yeh [6], where a parabolic electron density profile is chosen such that it approximates the actual profile of the plasma sheath encountered for a reentry vehicle. A spherical model uniformly coated by plasma has been used to calculate the backward cross section of a plasma-clad sphere as a function of the normalized plasma thickness [18]. The scattering characteristics of a blunt cone covered by an inhomogeneous plasma sheath are investigated from S-band to Ku-band frequencies by using the physical optics (PO) method [19], where both the collision frequency and the plasma frequency are assumed to be sine distribution.
With the rapid development of high-performance computers, various efficient algorithms were proposed that they enable us to solve the Navier-Stokes (NS) equations in the computational fluid dynamics (CFD) and the Maxwell equations in the computational electromagnetics (CEM) numerically. EM wave scattering characteristics in time-varying and spatially nonuniform plasma sheaths was investigated by Chen et al. [20] using the finite-difference time-domain (FDTD) method. Nevertheless, the collision frequency of the plasma sheath in [20] was calculated empirically and the variation of RCS is only summarized at the L and S wave bands. RCSs of a 10 cm diameter and 30 cm height metal cone with and without plasma were studied by Chung [21] for the S band and X band using FDTD. However, the simulation parameters were not chosen strictly and a simplified physical model was applied that limited its application in practical cases. Based on the hypersonic flow field data obtained in the repeatable laboratory simulations using ballistic range [11] or plasma wind tunnel [12], more practical physical models of hypersonic flow field can be established using the CFD technique [17]. By adopting a CFD tool to derive the distributions of electron density and atmospheric temperature around the targets, backscattering RCS of a hypersonic perfect electrical conductor (PEC) cone with a plasma sheath were investigated by Qian et al. [22, 23] using the volume-surface integral equation method incorporated with the multilevel fast multipole algorithm (MLFMA). However, only the properties at certain frequencies in the L-band were predicted. To reveal the effects of hypersonic flow field on the scattering characteristics of a reentry vehicle in wide wave bands, backward RCSs of a reentry vehicle enveloped by a hypersonic flow field are analyzed using a high-order auxiliary differential equation finite-difference time-domain (ADE-FDTD) method in this paper, where the hypersonic flow field is established by solving the NS equations numerically under different flight conditions.
The rest of the paper is organized as follows. In Section 2, descriptions of the modeling of a hypersonic flow field and the high-order ADE-FDTD method are briefly revisited. Numerical analysis of EM scattering properties of a PEC sphere enveloped by a hypersonic flow field under different flight conditions is presented in Section 3. Section 4 serves as a conclusion.
2. Method and Theoretical Treatments
2.1. Modeling of the Hypersonic Flow Field
Without any assumption on the distributions of the plasma parameters, the flow field around a hypersonic vehicle is established by numerically solving the NS equations in this paper. The flow field around the hypersonic vehicle is thermochemical nonequilibrium, governed by NS equations with chemical reaction source terms, accounting for the mass, momentum, and energy conservation laws. The Park two-temperature chemical-kinetic model [24] is applied in our simulations, where one temperature is assumed to characterize the translational and rotational energies of neutral, and another temperature for the molecular vibrational, electron translational, and electronic excitation energies. Concerning the chemistry model, the seven-species air model is applied, where the gas in the flow field is chemically regarded as a mixture of partially ionized air with seven species (, , , , , , and ).
The 3D nondimensional Navier-Stokes (NS) equations that include thermochemical nonequilibrium-related terms can be written in the following form:
The vector of conserved quantity is described by where and are the density of species and the density of mixed gas, respectively. , , and denote velocity components in , , and directions. and are the total internal energy per unit mass and vibrational energy per unit mass, respectively.
The inviscid flux vectors , , and in , , and directions are given, respectively, as where and are the pressure and enthalpy.
The viscous flux vectors , , and in , , and directions are expressed, respectively, as where is the Lewis number, is the Prandtl number, and , , and are the heat conductivity coefficient, specific heat capacity, and species mass fraction, respectively. represent the components of the stress tensor, denote the translational-rotational heat in , , and directions, and describe the vibrational heat flux terms in , , and directions.
The source vector has the following form in Park’s two-temperature model: where is the individual species mass formation rate and is the vibrational source term.
The NS equations are solved by the Advection Upstream Splitting Method by Pressure-Based Weight function [25] (AUSMPW+) based on the finite volume method. The diffusion terms are the central difference. The implicit lower-upper symmetric Gauss-Seidel relaxation [26] (LU-SGS) is taken as the time marching algorithm.
In the 1960s, a series of hypersonic experiments were performed to study the communications blackout in the RAM program. Among these flights, a large amount of data was collected during the reentry of the RAM-C II flight. The RAM-C II vehicle was a spherically blunted cone with a 9-degree cone half-angle, 0.1524 m nose radius, and a total length of 1.3 m, whose geometry is shown in Figure 1(a).
(a)
(b)
(c)
To verify the physical models and the CFD calculation schemes used in this paper, temperature distributions along the stagnation streamline and peak electron number density along the body are calculated for a blunt cone using the same parameters of the RAM-C II flight model [27]. The RAM-C II vehicle is assumed to be at 61 km height with a velocity of 23.9 Ma. Distributions of temperature along the stagnation line are displayed in Figure 1(b). Results published in [28], which were obtained by solving the NS equations numerically relying on a multiblock finite volume scheme based on a second-order accurate total variation diminishing, are also plotted in Figure 1(b) for the purpose of comparison. It is worth mentioning that Tchuen and Zeitoun et al. [28] analyzed the flowfield thermally by using a model with four temperatures (i.e., translational–rotational temperature , vibration temperatures of nitrogen and oxygen , , and electron–electronic excitation temperature ). As all the molecules are almost dissociated by the strong heating behind the shock wave, molecules become the only molecular species that dominate the vibrational energy in the shock layer [29]. Thus, the distributions of translational–rotational temperature and nitrogen vibration temperature along the stagnation streamline in [28] are cited in Figure 1(b) for comparison with that of and in Park’s two-temperature model. As shown in Figure 1(b), the temperature distributions along the stagnation streamline predicted using our method are in good agreement with those in [28], which partially indicates the correctness of our method. The peak electron number densities along the body are displayed in Figure 1(c). Experimental data and numerical results given in [28] are also plotted for the purpose of comparison. As shown in Figure 1(c), satisfied agreements are achieved although some differences exist.
Once the flow field parameters, for example, electron concentration, neutral particle density, and temperature, at any point in the calculation domain are determined, the complex relative dielectric constant of nonmagnetic plasma can be expressed using the Drude dispersion model as (the time dependence is described using the convention in this paper) where is the angular frequency of the incident wave. and are the plasma angular frequency and collision frequency of electrons with neutral particles, respectively, which are key parameters in the EM scattering characteristics of the plasma. The plasma angular frequency is calculated by [30] where is the electron mass and is the electron number density.
Concerning , an empirical formula has been used for years to calculate the collision frequency of an aircraft plasma sheath [20, 23]. Nevertheless, this formula was originally applied to the analysis of plasma in the ionosphere at relatively low temperature; a relatively large deviation can occur when this empirical formula was used for air flow in thermochemical nonequilibrium over hypersonic vehicles [28, 31]. Based on the plasma kinetic theory, combined with real flow conditions, the collision frequency of electrons with neutral particles derived from hypersonic air flows in thermal and chemical nonequilibrium gives results [32]. where is the number density of species ( represents , , , , and ), is the Boltzmann constant, and is the collision cross section between electron and species . The collision cross sections obtained from experiments [33] are adopted in our calculations, which are arranged and displayed in Figure 2. The collision cross-section data calculated using reliable numerical fitting results presented in [27, 32] is adopted since experimental data is hardly available publicly.
It is worth mentioning that to enable an efficient numerical calculation, the object covered by a hypersonic flow field is meshed separately in the establishment of flow field and in the computation of EM scattering; that is, the flow field grids employed in CFD and the electromagnetic field grids used in FDTD are inconsistent. The distributions of the EM characteristics (i.e., plasma frequency and collision frequency) of the plasma are obtained using the flow field parameters (i.e., electron density, neutral particle density, and temperature) through converting flow field grids to electromagnetic grids based on an averaging scheme [23]. In order to enable the electromagnetic meshes to reflect the flow field information accurately, not only should the size of the FDTD meshes meet its own convergence requirement but also the grid size of FDTD should be in the same order of magnitude as the smallest grid size of the flow field.
2.2. Auxiliary Differential Equation Finite-Difference Time-Domain Method
To analyze the interactions between EM/light waves and plasma media, various simulation methods were developed, including the Wentzel-Kramers-Brillouin (WKB) method [34], physical optics (PO) method [19], scattering matrix method (SMM) [35], and finite-difference time-domain (FDTD) method [36–39]. Among these methods, the FDTD method is widely applied in analyzing the EM scattering and propagation characteristics of plasma media due to its simplicity and wideband capability. To avoid the complicated transformations between the -domain and time domain in -transform FDTD and complicated convolutions in PLJERC-FDTD, ADE-FDTD was proposed for the scattering characteristics of unmagnetized plasma [38, 39] using the auxiliary differential equation (ADE) scheme. The basic principle of the ADE technique is to transform the frequency-dependent constitutive relation in the time domain and utilize the inverse Fourier transform to obtain an ordinary differential equation. Although different treatments in ADE are required for different kinds of dispersive media, the arithmetic used in ADE is much easier than those used in other FDTD methods toward dispersive media [39]. Furthermore, instead of using the traditional second-order FDTD, the fourth-order spatial central difference approximation technique is applied in this paper in order to reduce the numerical dispersion error and achieve high accuracy [40, 41]. The high-order ADE-FDTD method is briefly introduced as follows.
Maxwell’s equations for isotropic nonmagnetic plasma sheath are described by [42] where and are the electric field and the magnetic field, respectively. is the magnetic permeability in free space. In the frequency domain, the permittivity of a Drude medium is given by where is the electric permittivity in vacuum. is the polarization function, which is given by
Substituting (5) into (4) yields where is the polarization current in the frequency domain
Discretizing (7) by using FDTD (2) and (4), we obtain
Substituting (6) into (8) yields
Applying operator transitions to (10) and implementing an integral on both sides result in
Discretizing (11) gives where
Substituting (12) into the second formula of (9) yields where
To attain the fourth-order spatial discretization of and , the following central finite-difference scheme is employed in Yee’s staggered-grid arrangement where represents the electric components and magnetic field components ; in a three-dimensional FDTD Cartesian lattice, it takes the form as , with , and the spatial and temporal increments, respectively. The temporal discretization is the same as a conventional second-order FDTD. Then (9), (12), (14), and (16) are iterative equations for the and components of the high-order ADE-FDTD.
It should be noted that the complex relative permittivity can be written as where its imaginary part is
It can be known from (10) and (26) that the real and imaginary parts of the relative permittivity of the nonmagnetized plasma are described by
In our high-order ADE-FDTD algorithm, we have , , and for the grids in the plasma region, and , , and for the grids in the PEC sphere region.
A validity of the high-order ADE-FDTD method used in this paper is tested by calculating the backward RCS of a nonmagnetic plasma sphere with radius . The plasma frequency is , and the collision frequency is . The incident wave is assumed to be a Gaussian pulse which takes the form as
The results obtained using the high-order ADE-FDTD method in this paper are displayed in Figure 3, where the results calculated using the Mie theory and using the original ADE-FDTD are also displayed for the purpose of comparison. As shown in Figure 3, a very good agreement is achieved between the results obtained using the high-order ADE-FDTD and the benchmark results obtained using the Mie theory, which partially indicates the correctness of our code. Furthermore, compared to the original ADE-FDTD, results with higher accuracy can be achieved by using the high-order ADE-FDTD.
3. Results and Discussions
In this section, numerical simulations are implemented to analyze the features of a hypersonic flow field under different flight conditions and its EM scattering properties when illuminated by a plane wave with different incident angles. A PEC sphere is taken into the study, with an aim at providing useful reference data for practical tests in ballistic target or in the high-frequency plasma wind tunnels, where a sphere target is used due to its simple shape [13]. The radius of the sphere is . An illustration of a sphere surrounded by a hypersonic flow field is displayed in Figure 4(a); the object is assumed to be illuminated by a plane wave whose incident angle is defined in Figure 4(b).
(a)
(b)
(c)
In the simulation of CFD, the flow field is symmetric along the stagnation line. As a result, the simulation region size is 136.8 × 115.2 × 57.6 mm and a total of 540,000 irregular hexahedron grids have been generated. The freestream air gas is composed of 79% and 21% . According to the hypersonic flow properties, the inlet boundary is set to the freestream condition. The outflow is supersonic, and the zero gradient exit condition is appropriate. The no-slip and no-temperature jump conditions are used. The wall temperature is fixed (i.e., ), and a zero normal pressure gradient is imposed. At the wall which is supposed to be chemically noncatalytic, the flow is generally assumed to be in thermal equilibrium. In the simulation of CEM, the computational grid employed in the FDTD method is uniform with the grid size of 0.4 mm. The simulation region size is , and the total number of cells is 44,783,648. The transverse cross-section plane of the electromagnetic simulation region is shown in Figure 4(c). The boundary conditions are set as follows: connective boundary , output boundary , absorbing boundary , and truncation boundary . The truncation boundary is set to PEC. The uniaxial perfectly matched layer (UPML) [43, 44] with a thickness of 4 mm (i.e., 10 times the grid size) is employed in our simulations. For more details on the handling of boundary in FDTD, we kindly ask the readers to refer to [45, 46].
3.1. Distributions of Plasma Frequency and Collision Frequency
As preconditions for calculating the scattering properties, distributions of plasma frequency and collision frequency of the plasma flow surrounding a PEC sphere under different hypersonic conditions are calculated firstly. The flow field model is established by computational fluid dynamics technique, using the Park two-temperature kinetic model as the temperature model and the seven-species reaction model as the chemistry model. To reveal the influence of flight height and velocity on the formation of the plasma, the flight height is varied in the range of to and the velocity is varied in the range of to . Distributions of plasma frequency and collision frequency under different flight conditions are displayed in Figures 5 and 6, respectively.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
As shown in Figures 5 and 6, generally speaking, both and have the highest values near the head of the sphere, while the values decrease rapidly along the shoulders of the sphere. Compared to the plasma formed in a limited region near the head of the sphere, a quite broad area of plasma can be found in the wake of the sphere.
As the flight height increases from 30 km to 50 km, both the maximum values of and , which occur in the stagnation region, decrease dramatically, especially for which decrease from Hz to Hz. This is reasonable since the components, density, and temperature of the air reduce with an increase in the altitude, resulting in a decrease in the intermolecular collisions and a slow chemical reaction rate. Therefore, the electronic concentration decreases, leading to a decrease in the plasma frequency and the collision frequency.
As the flight velocity increases from 18 Mach to 24 Mach, although the increases of the maximum values of and at the object head are not so prominent, the formation of plasma in the wake of the sphere is greatly enhanced. In general, as the Mach number increases, the compression of gas by shock wave enhances, and more violent chemical reactions occur which strengthen the plasma around the vehicle. Thus, the electron density and the ionization degree increase, which makes more plasma regions to be dense.
Based on the distributions of the plasma frequency and the collision frequency , the complex relative dielectric permittivity of the nonmagnetic plasma in (1) can be rewritten as where we have
As we can see from (18) and (19), if , we have and . For an overdense plasma, that is, , we have , which indicates that EM waves with frequency will be completely reflected by the plasma. While for an underdense plasma, that is, , we have . It means that EM waves with frequency can propagate without attenuation in the plasma. However, we should keep in mind that is always not equal to zero due to the existence of the collision effect in the plasma, which inevitably leads to the energy loss of EM waves.
3.2. Backward Radar Cross Sections
Based on the distributions of the plasma frequency and the collision frequency obtained in the CFD, backward scattering characteristics of a PEC sphere enveloped by a hypersonic flow field are investigated using the high-order ADE-FDTD method in this subsection. The influence of the plasma flow on the backward RCS of the sphere under different flight conditions is analyzed. The incident EM wave applied is assumed to be a Gaussian pulse as defined in (17); the incident angles are defined in Figure 4(b). In our simulations, is fixed and is varied to indicate an arbitrary incident angle of the plane wave. Distributions of backward HH RCS components against frequency for different incident angles, including , , and , are displayed in Figures 7–9, respectively. Frequencies of incident wave in the range of 0.1 GHz~20 GHz are considered. The backward RCS of a static PEC sphere is also plotted for the purpose of comparison.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
As shown in Figure 7, in the low-frequency band, say less than 2 GHz, the plasma flow enhances the backward RCS of the sphere greatly in all the cases under study, while the backward RCS fluctuates with frequency in other wave bands. As displayed in Figure 7(a), for the vehicle flying at 18 Mach, an increase in the flight height from 30 km to 50 km has a small influence on the RCS. While the influence of flight height on the RCS is greatly enhanced as the Mach number increases, as shown in Figure 7(b), large changes occur in the range of 10 GHz~14 GHz with an increase in flight height for a velocity at 24 Mach. At flight height km, the variation of the Mach number changes the RCS greatly, as can be observed in Figure 7(c), but this influence becomes smaller as the flight height increases from 30 km to 50 km, which is mainly due to the decrease in the air density.
As shown in Figure 8, the behaviors of backward RCS for the case of incident wave at are totally different from that at . In most of the wave bands considered here, the plasma enhances the backward RCS of the sphere greatly. This is because compared to the case of a head-on incidence at , the large-scale plasma formed in the wake region contributes directly to the RCS for the case of the side illumination at . It indicates that although the values of and in the wake region are much smaller than those in the stagnation region of the sphere, the plasma wake lays significant influence on the EM scattering properties due to their large-scale distribution.
Compares to the RCS in the case of a head-on incidence at , the RCS increases obviously as the velocity increases, while the influence of the change of the flight height is relatively small. When the flight height of the target increases from 30 km to 50 km at a fixed velocity, the backward RCS of the vehicle enveloped by the flow field decreases with increasing height at most frequency bands. The reason is that the air density and temperature decrease with an increase in altitude, the plasma flow around the sphere becomes weaker, the reflection and absorption effects descend, and the penetration of EM waves in plasma flow enhances with increasing height. When the Mach number increases from 18 to 24 at a fixed flight height, the backward RCS of the target increases steadily in the range of 4 GHz~16 GHz. It can be interpreted that when the Mach number becomes larger, the plasma flow around the sphere becomes thicker because the chemical reaction in the flow field becomes more violent. The electron density increases with the Mach number, resulting in more regions with dense plasma, and the penetrability of the incident EM wave through plasma flow becomes worse.
It is interesting to find out in Figure 8 that as the frequency of the incident wave increases, the RCS of the sphere covered by a plasma flow gradually approaches to that of the original PEC target, which indicates that the plasma does not affect the backward RCS of the vehicle seriously at high-frequency bands. Furthermore, since a thinner layer of plasma is generated in the case of smaller Mach number or larger height, the RCS of the sphere covered by a plasma flow approaches faster to the original PEC target.
The distributions of backward RCS for an EM wave incident at are displayed in Figure 9. Compared to the results of the head-on incidence case displayed in Figure 7, the plasma lays much stronger influence on the RCS in this case, which might be due to the fact that thicker layer and larger scale of the plasma in the wake region contribute more to the RCS for the case of incidence at than that at . Similar to the case of illumination at , significant changes of RCS can be observed in the range of 4 GHz~12 GHz when the velocity is increased from 18 Mach to 24 Mach, as shown in Figures 9(c) and 9(d).
4. Conclusions
In this paper, backward scattering characteristics of a PEC sphere enveloped by a hypersonic flow field is investigated using a high-order ADE-FDTD method. The parameters of the hypersonic flow field around a PEC sphere are derived by solving the NS equations numerically. Then the distributions of plasma frequency and collision frequency are obtained by converting flow field grids to electromagnetic grids based on an average scheme. To reduce the numerical dispersion and improve the computational accuracy, a high-order ADE-FDTD, in which the fourth-order spatial central difference approximation is implemented, is applied to simulate the backward scattering characteristics of the target under different flight conditions. Although the application of the FDTD method to the scattering by arbitrarily shaped objects, such as blunt cone and pointed cone, is straightforward, it is restricted to an analysis of a sphere in this paper, since a sphere target is usually used due to its simple shape for practical tests in ballistic target or in the high-frequency plasma wind tunnels.
The effects of incident angle, flight height, and velocity on the scattering characteristics are discussed in detail. The results show that the hypersonic plasma flow has a significant influence on the backward RCS of the target. For frequencies below 2 GHz, the plasma region enhances the backward RCS in all the cases under study. The backward RCS fluctuates with frequency in other wave bands larger than 2 GHz if a head-on incidence of a plane wave is assumed. Nevertheless, the behaviors of backward RCS for an incident wave at are much different from those at . In most of the wave bands, the plasma enhances the backward RCS of the sphere greatly since a larger scale of the plasma in the wake region contributes directly to the RCS for the case of the side illumination at . The plasma region also lays stronger influence on the RCS for an incident wave at than that from the head of the target at .
Data Availability
The data in .dat form for Figure 2 are obtained from previously reported studies and datasets [33]. The processed data are available from the corresponding author upon request. The .dat data for other figures in this paper used to support the findings of this study are also available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (Grant nos. 61501350 and 61675159), the State Key Development Program of Basic Research of China (Grant no. 2014CB340203), and the National Natural Science Foundation of Shaanxi Province (Grant no. 2018JM6016). This work is also partially supported by the 111 Project (B17035).