AIM-E: E-Region Auroral Ionosphere Model

: The auroral oval is the high-latitude region of the ionosphere characterized by strong variability of its chemical composition due to precipitation of energetic particles from the magneto-sphere. The complex nature of magnetospheric processes cause a wide range of dynamic variations in the auroral zone, which are difficult to forecast. Knowledge of electron concentrations in this highly turbulent region is of particular importance because it determines the propagation conditions for the radio waves. In this work we introduce the numerical model of the auroral E-region, which evaluates density variations of the 10 ionospheric species and 39 reactions initiated by both the solar extreme UV radiation and the magnetospheric electron precipitation. The chemical reaction rates differ in more than ten orders of magnitude, resulting in the high stiffness of the ordinary differential equations system considered, which was solved using the high-performance Gear method. The AIM-E model allowed us to calculate the concentration of the neutrals NO, N( 4 S), and N( 2 D), ions N + , N 2+ , NO + , O 2+ , O + ( 4 S), O + ( 2 D), and O + ( 2 P), and electrons Ne, in the whole auroral zone in the 90‒ km altitude range in real time. The model results show good agreement with observational data during both the quiet and disturbed geomagnetic conditions.


Introduction
The Earth's ionosphere E-layer in the auroral zone is one of the primary indicators of solar influence on our planet. Satellite and ground-based ionospheric observations have limited temporal and spatial resolution, therefore the need for mathematical modeling arises. Models are widely used to analyze the ionosphere state when it is necessary to know the ionospheric parameters for any location and time.
Ionosphere numerical modeling provides a quantitative description of the ionospheric parameters and their spatial-temporal variations. These parameters are electron density, ion composition, ion and electron temperatures, ion velocity, and other parameters, depending on the modeling tasks. The ionospheric models can vary greatly in terms of computational accuracy, temporal and spatial resolution, computational costs (processor time, memory, etc.), purpose of use, and accessibility. There are a number of different ionospheric models described in the literature, but it should be noted that there is no universal approach that simultaneously satisfies all the requirements. The decision of which model to use depends on the goals and needs, computing resources, and the input data availability.
The E-region models, in most cases, are part of the global ionosphere or upper atmosphere models. The ionosphere global numerical models are a complex set of programs for solving the nonlinear differential equation systems. It should be noted that not all global models provide reliable results for the E-region at high latitudes. Depending on the given conditions, the global models can evaluate variations of the ionospheric parameters associated with solar activity cycles, seasonal and diurnal variations, and geomagnetic activity. Usually, the global models are computationally expensive and require supercomputers for global simulations. Therefore, such models are adequate for solving theoretical problems or for case studies, but not for the real time ionosphere monitoring.
Empirical ionosphere models are based on aggregated experimental data obtained by various methods. Usually, they are quite simple in calculations and can be used for real time diagnostics. However, the empirical models a priori do not include any dynamics and can demonstrate only a sequence of average states, and can satisfactorily describe only regular patterns of the ionosphere.
Another weakness of empirical models is the nonuniformity of the data coverage used to train the model. The most commonly used empirical model of the ionosphere is the global model International Reference Ionosphere (IRI), recommended by the Committee on Space Research (COSPAR) and the International Union of Radio Science (URSI) [1,2]. IRI-2016 is the latest modification of the model [3]. Like any empirical model, IRI is most accurate in regions with the densest observation coverage. Since most of the ionospheric data were obtained from the European and North American latitudes, these areas are best represented in the IRI. The Northern Hemisphere and continents are better described than the Southern Hemisphere and oceans, again due to obvious differences in data volume [4]. The IRI model can be used primarily for the middle latitudes. At the auroral and polar latitudes, the model based on observations from a small number of ground-based stations and satellites insufficient for an accurate modeling of this highly variable region. At the high latitudes, IRI can be applied only for rough evaluation of background ionospheric parameters, due to both low data coverage and high variability of the region considered [5]. A comparison of the IRI-2007 model with the vertical sounding data [6][7][8] shows that the model does not provide accurate calculation results for ionospheric parameters in disturbed geomagnetic conditions in the auroral zone and cannot be used for the auroral E-region.
There are a dozen models [9][10][11][12] that are used to study local processes in the auroral zone. For example, the Physicochemical Model of the Auroral Ionosphere [10] is a numerical model describing processes of the main excited and ionized atmospheric components interaction during auroral electron precipitation. The time-dependent, one-dimensional auroral model used in [11] is a combination of an electron transport code [13] and an auroral ionosphere kinetic model [14,15] evaluating the key E-region ions and minor neutral species over an altitude range 75-500 km. The thermosphere-ionosphere Fe/Fe + (TIFe) model [16] was developed for the investigation of formation mechanisms of thermospheric Fe layers observed by lidar in Antarctica. Those models are intended mostly for the case studies, and they show good calculation results and accurately describe the local physical and chemical processes of the auroral ionosphere.
We present an E-Region Auroral Ionosphere Model (AIM-E) that calculates the ion composition and electron concentration in the high-latitude E-region, including the auroral zone using as input parameters the auroral electron fluxes: (1) measured directly along the satellite trajectory using directly measured electron precipitation flux; or (2) their distribution and spatio-temporal variation using the empirical precipitation model. The AIM-E ionosphere model's advantage is flexibility of the input form of the photo-ionization source. The spectrum of solar extreme ultraviolet radiation can be specified in two ways: (1) using direct measurements of the photon flux spectrum from the such as TIMED satellite, for example; and (2) with the model EUVAC spectra, parameterized by the value of the daily index F10.7. The high-performance computation method allows the computer time to be reduced; therefore, the one-dimensional model can be used for fast calculations, and it shows the spatial distribution of electron concentration and other ionosphere parameters over the entire auroral zone.

Materials and Methods
The main scope of the AIM-E model is real-time diagnostics of the high-latitude ionospheric E-region. AIM-E provides temporal and spatial density variations of the minor neutral species NO, N( 4 S), N( 2 D), ions N + , N2 + , NO + , O2 + , O + ( 4 S), O + ( 2 D), O + ( 2 P) and electrons Ne. Figure 1 shows the general structure of the AIM-E model. Model input parameters are the solar EUV and auroral electron flux for the given date, time, and geographical coordinates, which can be set on the grid or in arbitrary locations. Neutral atmosphere composition and temperature are also input parameters of the model and they were taken from the NRLMSISE-00 model [17]. Figure 1. AIM-E model structure. Input parameters-date, time, and location-are used to calculate the altitude distribution of the ionization rates, utilizing three empirical blocks: NRLMSISE-00 for the neutral composition, EUVAC for the UV solar radiation spectrum and OVATION-Prime for the precipitating electron spectrum. Calculated distribution of ionization rates then go to the aeronomy block-the ODE solver for the system of continuity equation for 10 ionospheric species: neutrals NO, N( 4 S), and N( 2 D) and ions N + , O2 + , N2 + , NO + , O + ( 4 S), O + ( 2 D), and O + ( 2 P).
In quiet geomagnetic conditions during the daytime, the main ionization source, responsible for regular E layer formation, is EUV solar radiation in the wavelength range from 5 to 105 nm. For the night conditions, it is extremely important to take ionization by electrons of magnetospheric origin with energy 10-30 keV into account. During disturbed periods, sporadic ionization by the precipitating electrons plays a dominant role. In the AIM-E model, both the EUV solar radiation and the electron precipitation flux can be specified in two ways: (1) from the direct satellite measurements; and (2) from the empirical, the EUVAC model for solar radiation [18] and OVATION-Prime for electron precipitation [19]. Any other solar radiation EUV and energetic electron precipitation models can be used as AIM-E input source, depending on the modeling tasks. Solar photons and magnetospheric electrons lose energy due to the absorption passing the atmosphere. This fact is taken into account using the extended Chapman function for the solar EUV radiation propagation [20] and the energy dissipation function for the precipitation electrons [21]. Photo-ionization (quv) and energetic electron (qcorp) ionization rates obtained for different ions come into a system of ten ordinary differential equations (ODE) (Appendix A). The initial solution of the ODE system is based on the equilibrium state of the E-layer, which reaches 39 chemical reactions between ionospheric species (Appendix B) in each location under specific solar and geomagnetic activity conditions. The numerical solution of the ODE system for neutrals and ions is implemented using the high-performance Gear method [22] (see Appendix D).

Photo-Ionization
The photo-ionization rate q (z) of the neutral gas j-th component is the number of photo-ionization acts per unit volume per unit of time at the altitude z, is determined by the expression [23]: where n -the concentration of the gas j, σ -the photoionization cross-section of the gas j by solar radiation at wavelength, λ, σ -the photo-absorption cross-section of the gas j by solar radiation at wavelength, λ, I -the photon flux density at wavelength, λ at the top of the atmosphere, Ch(χ)-the Chapman function.
The photon flux at the top of the atmospheric, I is the number of photons per unit area located perpendicular to the radiation direction per unit of time. The photo-absorption and photoionization rates for the neutral atmospheric components depend on the solar radiation wavelength. The photo-absorption and photo-ionization cross-sections for O, N2, and O2 were taken from [24]. Quantum yields of different states the atomic oxygen ions were taken from [25]. The numerical values of the cross-sections are given in Appendix C.
Calculation of the electron density using the auroral ionosphere model for the daytime can be done using two different methods of solar EUV radiation input: (1) a theoretically calculated EUV spectrum using EUVAC Solar Flux Model [16] parametrized on F10.7 index [26] (see example in Section 3.1); and (2) direct measurements of the solar EUV spectra obtained from spacecraft.

Electron Precipitation Ionization
During the disturbed periods, such as magnetospheric storms and substorms, ionospheric ionization in the auroral zone caused by particle precipitation can be several orders of magnitude higher than photoionization. Thus, it is crucial to take into account the energetic electron precipitation for adequate modeling of ion composition in the high latitude ionosphere.
The ion production rates for ionization of the atmospheric neutrals by auroral electrons are calculated using the approach described in [21].
According to Sergienko and Ivanov method [21] the ionization rate qj(z) (cm −3 s −1 ) of the certain ion j at the altitude z is defined by: where: W(z) is the energy deposition function of the auroral electrons (eV cm −3 s −1 ), Pj(z) is the deposited energy fraction expended in the excitation and ionization of the neutral j, and εj(z) is the "energy cost" for the ionization of the neutral j, that is, the average energy expended in the creation of an electron-ion pair of the neutral j (eV). For mono-energetic electron flux with electron energy of E (eV) the energy deposition function can be expressed in term of the dimensionless dissipation function as: where: z is altitude (km), ρ(z) is the atmospheric mass density (g cm -3 ), s(z) is the distance from the top of ionosphere z0 to the altitude z in mass density units (g cm −2 ), R(E) is the average range of an electron in air (g cm -2 ), Λ is the dimensionless dissipation function, E is the Energy (eV), and Alb(E) is an albedo flux, that is, the dimensionless function describing the part of the total energy contained in the initial electron flux that is back-scattered by the atmosphere. For an arbitrary theoretical or experimental electron differential flux φ(E) at the top of the ionosphere the energy deposition function can be expressed by the integral: φ(E) is the differential auroral electron flux at the top of the ionosphere (700 km) (eV −1 cm −2 s −1 ). All parameters in Equation (3) were taken from [21]. The production rates of the ion species involved in the model are calculated by substituting the Equation (4) into Equation (2) and using the "energy cost" presented in Table 1.

Aeronomy Block
The AIM-E module of ionospheric aeronomy is a solver of a system of ten first-order nonlinear ordinary differential equations (Appendix A). The equations represent the nonstationary one-dimensional continuity equations for ten species included in the model: where N is the concentration of species i, Q , and L are the production and loss rates of species i respectively at the altitude z. The diffusion term is omitted in the continuity equations, because the transport effects are negligible for the E-region ionosphere [27]. The production rates Q and the loss rates L include photo-ionization, ionization, and excitation by the auroral electron precipitation and chemical reactions.
The The system is solved using an iterative procedure, which continues until the state of chemical equilibrium is reached for the given input conditions. The chemical reaction rates differ by more than 12 orders of magnitude. Such a wide range of the coefficient values results in high stiffness of the ODE system. Since the stiff ODE cannot be efficiently solved using explicit numerical schemes, a number of implicit methods were specially designed for such systems. Here we describe the application of the 4-th order iterative implicit scheme from the family of Gear methods [22], which is briefly summarized in Appendix D. Employment of the Gear method with step size adjustment dramatically decreased the computation costs when solving the stiff ODE system.

Results
In this section we present three examples of case studies that present the possible scope of applications of the AIM-E model: (1) monitoring of the regular E layer parameters; (2) calculation of the electron concentration vertical profiles along the satellite trajectory; and (3) evaluation of global electron distribution in the high latitude ionosphere using EUV and particle precipitation models as inputs.

Monitoring of the Regular E-Layer Parameters
During the day, the ionospheric regular E-layer plays an important role in the HF radio wave propagation . Knowing the current conditions and forecast of the regular E-layer is critical to ensuring long-distance logistics on land, at sea, and in the air.
A 3-day period of quiet geomagnetic conditions, 8-10 June 2017, was used to prove the ability of the AIM-E model to predict the main characteristics of regular E layer (mean Kp index for this period is equal 0.7 and mean AE = 67 nT). We compare the measured and calculated critical frequencies for a number of ground-based stations at auroral and subauroral latitudes-Gorkovskaya (GRK), Salekhard (SAH), Lovozero (LOZ), Pevek (PBK), Amderma (AMD), Tiksi (TIK), and Dikson (DIK). Geographic and corrected geomagnetic coordinates of the stations are given in the Table 2, the vertical sounding data provided by the Polar Geophysical Center of the Arctic and Antarctic Research Institute (https://geophys.aari.ru, accessed on 8 June 2021). A negligibility of the sporadic Es layer during this period allowed us to determine 406 values of critical frequency for the regular E layer (foE) from 504 available hourly ionograms. For each time and observation point the AIM-E model with the EUVAC photon flux input was used to calculate the vertical profiles of the electron concentration in the altitude range from 90-130 km with 1 km step size. Using the EUVAC model parametrized on the F10.7 index as the AIM-E model input allowed us to calculate EUV solar radiation fluxes in the wavelength range from 5 nm to 105 nm, including 36 spectral lines.
The critical modeled frequency in MHz was determined from the value of the E-layer maximum electron concentration Ne [m −3 ] as: Figure 2 shows the diurnal variation of the foE measured (red lines) and modeled (blue dashed lines) foE during the considered 3-day period for the selected auroral and subauroral stations. The modeled foE values were in good agreement with measured values at all stations (correlation coefficient R = 0.982). This result shows that the AIM-E model can adequately describe the current state (and even a short-term forecast) of the ionosphere regular E layer, which is especially important in the absence of ionospheric data or during the strong D-region absorption events.

Calculation of the Vertical Electron Concentration Distribution along the Satellite Trajectory
In this section we present the calculation of electron concentration altitude profiles along the REIMEI satellite [28] pass across the auroral oval on 5 December 2007. During this intense electron precipitation event, the REIMEI satellite was in magnetic conjunction with EISCAT radar at 0:36:36 UT (0.610 h), enabling quantitative verification of the AIM-E model. The REIMEI spacecraft has a near-sun-synchronous orbit at 12:50-0:50 LT meridian and its altitude range is 610-670 km. Projection of the spacecraft trajectory along the IGRF-12 magnetic field [29] to the 110 km altitude is shown in Figure 3a. Figure 3b shows the REIMEI satellite differential electron energy flux measured in the energy range from 60 eV to 11   We used the EISCAT Tromso UHF incoherent scattering radar observations (http://portal.eiscat.se/madrigal/, accessed on 6 June 2021) to perform a quantitative verification of the AIM-E results. The radar is located in Tromso, Norway (69°35′ N, 19°13′ E). The EISCAT experiment was carried out simultaneously with the REIMEI satellite passage and the moment of the conjugated measurements was at 0:36:36 UT (0.610 h). The EISCAT UHF radar operated the radar program Arc1 (altitude range: 96-422 km; spatial resolution: 0.9 km; temporal resolution:0.44 s) with the antenna directed to the magnetic zenith. Figure 4 shows a sufficiently good agreement between the altitude profile of the electron concentration calculated with the AIM-E model (red line) and measured by incoherent scattering radar (blue line). The discrepancy between model and radar data can be explained by the fact that the observation of precipitating electrons by REIMEI spacecraft was not performed exactly over the incoherent scatter radar. The spatial distribution of the electron precipitation in the auroral region is very dynamic and has irregular structure, and the difference of even several tens of kilometers may lead to the observed uncertainty in the model results. The results obtained during the quiet geomagnetic conditions (Section 3.1), when solar EUV radiation dominates, indicate that the ODE system solution and the calculation methodology are correct.

Global Distribution of Ionospheric Composition
Using the empirical models of the solar EUV radiation and auroral electron precipitations, the AIM-E model allowed us to obtain spatial distribution of ionospheric parameters regardless of the availability of direct satellite measurements. An example of an electron concentration map, shown in Figure 5, demonstrates an ability of AIM-E to simulate the global distribution of ionospheric composition using modeled input parameters. For these calculations, the empirical EUVAC and OVATION-Prime models were used to provide the input photon and electron fluxes at the top of the ionosphere, respectively. We adapted the OVATION-Prime model, which is based on the DMSP measurements and provides electron precipitation parameters in the high latitudes. Using the OVATION-Prime results as an input for the AIM-E model, we could calculate an instant ion and electron composition everywhere in the high-latitude ionosphere for the wide range of geophysical conditions. The OVATION-Prime model provides only three general precipitation parameters for each grid point-total energy flux, total electron flux, and average electron energy. The latter two quantities are used to reconstruct the energy spectrum of precipitating electrons assuming the Maxwellian distribution.
The calculations were carried out for the northern hemisphere, MLAT = 50-90° N, on a discrete grid MLT × MLAT = 0.25 × 0.25 h in the latitudinal interval in a solar-magnetic coordinate system (SM). The transition between the geographic and SM coordinate systems was performed using GEOPACK [30].
For each grid point, the altitude profile of the electron density was calculated from 90 to 150 km with the 1 km altitude step. Figure 5  During the quiet conditions at 5:00 UT, the electron concentration in the night side part of the auroral oval reaches ~10 11 m −3 due to electron precipitation from the magnetosphere, which is comparable to the ionization due to the EUV radiation on the day side. During the main phase of the storm at 17:00 UT, the situation changes drastically. The auroral oval expands to the low latitudes, and the maximum of electron concentration in the E-region increases by an order of magnitude, which demonstrates the dominant effect of corpuscular ionization in the auroral oval during geomagnetic disturbances. The movie for the diurnal dynamics of the maximum electron concentration during the geomagnetic storm on 17 March 2013 is attached in the supporting materials (Movie S1).
The AIM-E model makes it possible to obtain the spatial distribution of ionospheric parameters in the entire auroral zone under various geomagnetic conditions that may be useful when and where the direct ionospheric observations are absent.
Proper spatial distribution of the ionosphere E-region parameters is extremely important for solving practical and theoretical problems of the radio wave propagation. During the night periods, when the F2 layer degrades, the E-layer plays a major role in the radio communication on the short (up to 2500 km) single-hop radio tracks. This is especially relevant during the polar night. In such cases the global modelling of the E-layer can support the available observation data for more accurate auroral ionosphere monitoring.

Discussion
In this work we present a detailed description of an AIM-E numerical model that calculates the chemical composition of the high-latitude ionosphere in the altitude range from 90 to 150 km. The AIM-E model was developed specifically for the high-latitude ionospheric E-region and takes into account the variability of the magnetospheric electron precipitations that is very important for the active geomagnetic conditions. The model calculates the number densities of the 10 species The model solves the system of stiff ordinary differential equations using the efficient 4th order implicit numerical scheme with variable integration step [22], which extremely reduces computational costs compared to commonly used explicit schemes. The high performance of the AIM-E model allows real-time calculation of the E-region ionosphere composition in the whole auroral zone under different geomagnetic conditions, taking into account high variability of the auroral electron precipitations.
In the photo-ionization module the solar EUV radiation spectrum can be set in two ways, by: (1) relatively rough empirical approximation using the daily F10.7 radio flux, continuously available since 1947 and suitable for 'space climate' investigations; and (2) direct measurements of the photon energy spectrum providing by, e.g., TIMED spacecraft every 97 min since 2002 which can be used for investigation of 'space weather' events, such as solar flares. The precipitating electron ionization module also has two input possibilities: (1) direct spacecraft measurements of the electron energy spectrum available from several low-altitude satellites, e.g., DMSP, NOAA POES, and REIMEI, providing an accurate local solution for vertical profile of chemical composition along the satellite trajectory; and (2) results of empirical model for electron precipitation (here we implement the OVATION-Prime model), which can be used for global scale simulation of the chemical composition in the high-latitude ionosphere.
The model results, calculated with different types of input parameters and for the different geomagnetic conditions, show a good agreement with ground-based ionospheric data (EISCAT UHF incoherent scattering radar and vertical sounding of the ionosphere). The model reproduced the ionospheric response to the geomagnetic substorm with a sufficient accuracy and can be used to describe the large-scale dynamics of the auroral oval during disturbed periods. The high computing performance allowed application of the AIM-E model in regular E-layer monitoring all over the auroral and subauroral zones, which is important for a real-time prediction of the radio wave propagation conditions. The AIM-E model can be applied to various practical tasks using input parameters such as solar EUV radiation and auroral electron flux from in situ measurements by satellites, or from the models, or as their combination depending on the modeling goals.
The AIM-E model was designed to become a useful tool for both operational and scientific applications in the Arctic and Antarctic ionosphere regions. The model can be used in a wide spectrum of scientific problems, including, but not limited to, various aspects of ionospheric chemistry, distribution of conductance, electric fields and currents, calculation of electron density, etc. In the operational mode, the model is suitable for realtime monitoring of the radio waves propagation conditions, which are crucial for the quality of communication in the most unstable auroral zone. The AIM-E model may also be utilized as a part of more complex space weather models as an efficient ionospheric module.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. AIM-E Model ODE Continuity Equations System for Ten Species
The system consists of ODE continuity equations for 10 ionospheric species (NO, N( 4 S), N( 2 D), N + , N2 + , NO + , O2 + , O + ( 4 S), O + ( 2 D), and O + ( 2 P)), including 39 ion-chemical reactions listed in the Appendix B. The system is solved using the high-performance Gear method (Appendix D).

Appendix C. Absorption ( ) and Ionization ( ) Cross-Sections for Neutrals O, N, O2, N2, and NO and Probabilities of the Ionization to the Electronic States of Atomic Oxygen Ions O + ( 4 S), O + ( 2 D), and O + ( 2 P).
The numerical values of the cross-sections for the neutral atmospheric components depending on the solar radiation for 36 spectral lines and channels are presented in the table below. The photo-absorption (σ ) and photo-ionization (σ ) cross-sections for O, N2,and O2 are taken from [24]; for N- [59]; for NO- [60]. Quantum yields of different states the atomic oxygen ions were taken from [25].
The cross-section values are presented in units of megabarns. In order to transform them to m 2 multiply by 10 −22 .

Appendix D. Numerical Scheme for the Continuity Equations System Solving
The ODE system solution is based on the equilibrium state of E-layer that reaches 39 chemical reactions between ionospheric species in each location under specific solar and geomagnetic activity conditions. The chemical reaction rates differ in more than ten orders of magnitude bringing to the high stiffness of the considered ordinary differential equations system, which requires the use of specific numerical methods. In the AIM-E model, the numerical solution of the continuity equations system for neutrals and ions is implemented using the high-performance Gear method [22,61], which implements a high-order implicit backward differentiation scheme. Here we briefly summarize the main points of the numerical method for solving the stiff ODE systems in the predictor-corrector formulation.
Predictor: (l , l , … , l ) for the differential scheme of the q-th order are given in Table A3. The qth order Gear scheme assumes the known solutions in the q previous time steps. In our case of 4th order scheme to start the numerical procedure we have to know Y ⃗ for s = 0,1,2,3. So the first three sets of l coefficient can be used to make the 1st, 2nd, 3rd time step using the lower order scheme. For the 0th step we use the simple forward Euler method. To update the history matrix to the next time step, it should be multiplied on the L×L prediction matrix Â: (ij) = ! !( )! -binomial coefficients. However, representation of the history matrix Z provides the advantage of avoiding the matrix multiplication by using simple replacement with repeated additions: Z , ← Z , + Z , , where n and n + 1 indices denote the current and next time steps, arrow is a replacement operator. This trick reduces the required amount of memory and speeds up the calculations.
Another key advantage of the method in terms of processing time is an adjustment of the step size during the calculations. The accuracy of the solution at every time step is controlled by checking the inequality: -weighted standard deviation of the approximate solution from the exact solution, EWT , = RTOL • Y , + ATOL -error weight vector elements, RTOL and ATOL -local relative and absolute tolerances empirically defined to be 1e-4 and 1e-6, respectively. Summation is going over the number of variables in the system.
The step size adjustment procedure is as follows. If the inequality (A4) is not satisfied, then reduce the integration step size h by 2 times, recalculate the history matrix Z (A3) using the algorithm below, and repeat the iteration procedure (A2). Otherwise, if the inequality (A4) is satisfied, then increase the integration step size h by r: