A SOC based avalanche model to study the magnetosphere-ionosphere energy transfer and AE index fluctuations

ABSTRACT Magnetosphere-ionosphere energy transfer and AE fluctuations are studied using a cellular automata model of terrestrial magnetosphere based on the concept of self-organised criticality (SOC). The model is a SOC-driven dissipative dynamical system with both spatial and temporal degrees of freedom. The input parameter to this model is derived from the real-time values of solar wind ion density and flow speed data. Both the direction and intensity of the real-time values of the BZ component of the interplanetary magnetic field (IMF) are the factors controlling the energy injection into the system. The model produces an output series which can be regarded as a mathematical representation of the AE index. The spectral response of the simulated output follows a 1/fβ power law, demonstrates a breakpoint at f0 = 0.050 mHz (5.5 hours) having slopes βA = 2.2–2.4 for f > f0 and βB = 0.9–1.0 for f < f0, the typical characteristics of the natural AE index. The entire 23rd solar cycle had been studied using the model. The parameter KA plays a significant role in the entire process. KA represents the remaining percentage of the released energy from the previous magnetosphere-ionosphere energy transfer, stored in the ionosphere.


Introduction
Magnetosphere-ionosphere interaction and the subsequent energy transfer is a significant phenomenon in magnetospheric dynamics. During a strong solar wind-magnetosphere coupling, a large amount of solar wind energy enters the geospace. A part of the energy is stored in the magnetotail and another part drives the convection of plasma particles in the magnetosphere. As the solar wind injection continues, the stored energy in the magnetotail reaches an unstable state, triggering magnetic reconnection. A huge amount of energy is released in the ionosphere causing magnetic fluctuations in the auroral zone. The auroral electrojet (AE) index is a global and instantaneous measurement of the auroral zone magnetic activities. The geomagnetic fluctuation in the horizontal component of the Earth's magnetic field H in the auroral region is measured in 10-13 observatories situated around the auroral zone. For each station, the average value of H of the five international quietest days of the month is considered as the base value of the measurement. For normalisation of data, the base value is subtracted from each data of the station. Then, all the normalised data from all the stations are plotted and superimposed on each other. The maximum and minimum deviations of H are termed as Auroral Upper (AU) and Auroral Lower (AL) index, respectively, which form the upper boundary and lower boundary of the envelope. If there are no disturbances from the distant axially symmetric fields or zonal currents, the AU and AL indices are the direct measurement of the maximum eastward and westward electrojet currents at any time. AE index is defined as the maximum total amplitude of the eastward and westward electrojet currents, that is, AE = AU -AL. As the AE index is the difference value of AU and AL indices, it is independent of the ionospheric zonal current or distant axially symmetric fields (Davis and Sugiura 1966). AE index is extensively used in aeronomy, solar-terrestrial physics, geomagnetism, and auroral studies.
The physical meaning of the AE index, the nature of the eastward and westward electrojets, the limitations of AU and AL indices have been discussed (Rostoker 1972(Rostoker , 2002Baumjohann 1982;Kamide and Kokubun 1996;Kamide and Rostoker 2004) along with a detailed study of spatial and temporal distributions of magnetic effects of the electrojets (Allen and Kroehl 1975), statistical analysis (Nakamura et al. 2015) and its relation with the polar cap index (Vennerstrøm et al. 1991;Vassiliadis et al. 1996). It has been also established that the AE index is subject to universal time variation (Davis and Sugiura 1966;Ahn et al. 2000a), seasonal variation (Ahn et al., 2000b;Cliver et al. 2000;Lyatsky et al. 2001;Russell and McPherron 1973;Li 2002, 2006), annual variation (Lyatsky et al. 2001;Pulkkinen et al. 2011), and solar cycle variation (Ahn et al., 2000b). Different dynamical and numerical models have been presented to predict and study auroral electrojets. The analogue model proposed by  and reviewed by (McPherron and Rostoker 1993), the Faraday loop model (Klimas et al. 1992(Klimas et al. , 1994, using linear prediction filter (LPF) technique (Bargatze et al. 1985) which is further modified as local-linear prediction technique by (Vassiliadis et al. 1995), by artificial intelligence (Hernandez et al. 1993;Lundstedt 1997, 2001;Gavrishchaka and Ganguli 2001;Weigel 2003;Chen and Sharma 2006), by stochastic approach (Pulkkinen et al. 2006) and also using solar wind parameters (Li et al. 2007;Luo et al. 2013).
Previously, we studied the characteristic structure of the Dst index, the global measurement of geomagnetic activities in Earth, and realised the series as a positively correlated fractional Brownian motion, displaying long-range correlation (Banerjee et al. 2011). But, we also understood the complexity of the non-linear dynamics of the geomagnetic fluctuations and the difficulties in predicting them. To gain an insight into the magnetospheric dynamics, we focused our study to develop a cellular automata model of terrestrial magnetosphere based on the concept of selforganised criticality and sandpile dynamics. According to the concept, a dissipative, dynamical system with both spatial and temporal degrees of freedom naturally evolved to a self-organised critical state without much specification of the initial conditions. The dynamical behaviour of a pile of sand is the most prominent example of this type of system. Let us start with a grain of sand and continue to add more grains to it, gradually forming a pile. The increment of the slope enhances the characteristic size of the largest avalanches. As a result, the equilibrium state of the sandpile is seriously disturbed. Eventually, when the slope becomes very large, the pile reaches a critical state. Now, adding a single grain of sand further to the pile collapses it. An avalanche of sand is released from the pile, the base area increases, and the system again returns to a state of equilibrium (Bak et al. 1987(Bak et al. , 1988. The concept of self-organised criticality (SOC) and sandpile model can be the basis of an analytical study of magneto-ionospheric dynamics. The solar wind, a stream of energised plasma particles is emitted from the outer atmosphere of the Sun and approaches the Earth. The interplanetary magnetic field (IMF) is trapped in the solar wind. Near the terrestrial space, the solar wind flow speed varies from a minimum of 260 km/sec to a maximum of 750 km/sec while the ion density varies in a much wider range, from 0.1 cm −3 to 100 cm −3 (Russell 2001). The density fluctuation is the primary controller of the variations in the dynamic pressure, which further controls the solar wind-magnetosphere reconnection. When the solar wind interacts with the terrestrial magnetosphere, this supersonic flow creates a standing shock wave or bow shock in the day-side of the magnetosphere and converts it into a subsonic flow. Majority of the energised solar wind plasma particles is heated and then get deflected around the Earth at the bow shock. The day-side magnetosphere is compressed down while the night-side magnetosphere is stretched up to about 100 Earth radii comprising the magnetotail (Nishida 2000;Borovsky and Valdivia 2018). The trapped interplanetary magnetic field (IMF) plays a crucial role in controlling the intensity and duration of solar wind-magnetosphere coupling. For a northward IMF B Z , the cusp shifts away from the equatorial region and further narrows down with the increasing intensity of B Z , decreasing the amount of energy injection into the magnetosphere. But, for a southward IMF B Z , the cusp shifts towards the equatorial region and widens out gradually as the intensity of B Z increases, triggering an injection of a large amount of solar wind energy into the magnetosphere (Lu et al. 2013). The magnetotail becomes a reservoir of this energy. As the solar wind injection continues, the magnetotail grows further by accumulating more and more energy. Eventually, the magnetotail growth reaches a critical point, becomes unstable and magnetic reconnection occurs in the tail. A part of the stored energy is released through flow kinetic energy and plasma heating, producing substorm (Borovsky and Valdivia 2018). Substorm originates deep in the magnetosphere, near the geostationary orbit (Antonova and Ganushkina 2000). It is a short but intense earthward convection of magnetic flux in the magnetotail which injects energised particles in the dipolar region, also substantially increasing the auroral electrojets . The Joule energy extracted from the magnetosphere is dissipated in the ionosphere through auroral electrojets, the fieldaligned currents flowing between the nightside magnetosphere and nightside ionosphere (Strangeway 2013). The subsequent magnetic fluctuations in the auroral region are measured by the AE index.
SOC has long been proposed as a possible explanation of magnetospheric dynamics. Sandpile model was selected as the first example displaying the concept of self-organised criticality, introduced by Bak et al. in their 1987 and1988 papers. Since then, the application of this theory has produced numerous significant analyses of magnetospheric activities by some eminent researchers. Consolini observed SOC-triggered behaviour in the power spectral density and burst size distribution of the AE index (Consolini 1997). Using the sign-singularity analysis, magnetic field fluctuations in the near-tail regions were investigated based on the concept of self-organised criticality and 2ndorder phase transition (Consolini and Lui 1999).
Chapman et al. presented a sandpile cellular automaton model to study magnetospheric dynamics (Chapman et al. 1998) and Klimas presented a firstorder physical model of plasma sheet . Further studies on modelling continued based on the SOC approach with interesting conclusions (Takalo et al. 1999;Uritsky and Semenov 2000). Other notable works discussed in detail the distinctive features of SOC-driven instabilities and their effects in the context of nonlinear dynamics of magnetosphere (Chang 1992(Chang , 1999Uritsky 1996;Uritsky and Pudovkin 1998a;Sitnov et al. 2000). (Dobias and Wanliss 2009) suggested that both the storm and substorm characteristics are consistent with the behaviour of the critical system and follow the fractal point process (FPP).
(Uritsky and Pudovkin 1998b) developed a twodimensional sandpile model of the magnetosphere current sheet to study the AE fluctuations. The model was a rectangular matrix of x and y dimensions. Each element of the matrix was characterised by an amount of energy. Also, a critical threshold of energy was assigned to all of the elements to determine the stability of the element after each energy injection. They investigated the avalanche formations, energy redistribution, and plasma sheet instabilities of the SOC-driven system in reaction to external disturbances. It was suggested that the spatially localised magnetotail instabilities can be regarded as SOC avalanches and the superposition of these avalanches of different sizes finally produces the characteristic low-frequency 1/f -like fluctuations of the AE index (Uritsky and Pudovkin 1998b). Uritsky et al. continued the study of the above two-dimensional sandpile model in their 2001 paper to explain geomagnetic substorms as a self-organised critical dynamic of the perturbed magnetosphere. In their study, the total accumulated energy in the system, as well as the energy dissipated from the system, were revealed to be the two major factors controlling the overall dynamics of the system. Moreover, the spectral characteristics of the model output showed striking similarities with the natural AE fluctuations (Uritsky et al. 2001). Based on the model (Uritsky et al. 2001), we developed a sandpile-like cellular automata model of Earth's magnetosphere in our previous paper (Banerjee et al. 2015). The model is a SOC-driven dissipative dynamical system with both spatial and temporal degrees of freedom. It is a two-dimensional array of finite dimensions and is characterised by energy E. The input energy to this model is derived from the real-time value of solar ion density and flow speed data. Both the direction and intensity of the realtime value of the B Z component of the interplanetary magnetic field (IMF) are the factors controlling the amount of solar wind energy injected into the system at any time. The total accumulated energy of the system is estimated to produce a simulated output representing the Dst index. The spectral characteristics of the simulated output closely follow the natural Dst index, establishing the acceptability of the model (Banerjee et al. 2015). We continued our study with the model investigating the solar wind-magnetosphere interaction, the injection of plasma particles, and the aspects of internal magnetospheric dynamics as a subsequent effect (Banerjee et al. 2019).
In the present paper, we extend our sandpile-like cellular automata model to study the dynamical behaviour of the auroral zone magnetic activities by focusing on the transferred energy from the system. The model, a representation of the Earth's magnetosphere, is a lattice of n × n elements, having spatial and temporal degrees of freedom. Each element is characterised by an energy E, which is analogous to the slope of the sandpile. The input to the model is solar wind energy, derived from the real-time value of solar ion density and flow speed data. The input energy is injected into the system through a set of elements, representing the cusp. Both the direction and intensity of IMF B Z are the factors controlling the width of the cusp, hence the amount of energy injection into the system. The injected energy is altering the potential energy of the elements of the lattice. Each element has a critical value of energy, known as the excitation threshold. As the energy injection continues, the energy gradually piles up in the elements, similar to the storing of energy in the magnetotail. If the energy of any element reaches the point of criticality by exceeding the excitation threshold, spatially localised instabilities are formed, representing the instability formation in the magnetotail. The pile collapses by releasing avalanches of energy in various sizes and shapes. The released energy is distributed among the adjacent elements and the process continues, gradually spreading throughout the lattice. The excess energy reaching the upper or lower margins is transmitted outside the lattice, representing the magnetosphere-ionosphere energy transfer. The amount of transferred energy at any time t is estimated. By some numerical calculations, an output time series is produced from this estimated energy. The simulated output can be regarded as a mathematical representation of the AE index. It is observed that the simulated series exhibits 1/f β -like power spectrum. The powerlaw exponents β A and β B of the power spectral density of the simulated output are evaluated for both the high-and low-frequency regions respectively for all the years of the entire 23 rd solar cycle. Finally, by comparing the values of β A and β B with that of the real-time AE index, it is observed that the simulated output closely follows the natural AE fluctuations depending on the exact value of the parameter K A . K A represents the remaining percentage of the transferred energy of the previous magnetosphereionosphere energy transfer, stored in the ionosphere. In our previous work (Banerjee et al. 2019), the realtime solar wind and IMF B Z data of the 23 rd solar cycle were used as the input to the model to investigate the solar wind-magnetosphere energy transfer. As we are continuing our study on the model, in the present work, we used the same input dataset to analyse the magnetosphere-ionosphere energy transfer.

Method and data
The cellular automata-based sandpile model used here is an extension of the model presented in our previous papers (Banerjee et al. 2015(Banerjee et al. , 2019 which is in turn based on the model presented by (Uritsky et al. 2001). The model, representation of the Earth's magnetosphere, is a finite matrix of n × n elements having spatial and temporal degrees of freedom. Each element is characterised by an energy E, which is analogous to the slope of the sandpile. E has an arbitrary unit. The threshold value of energy for each element is indicated as E TH , the excitation threshold (Uritsky et al. 2001).
The input to the model is solar wind energy.
The input energy dE is estimated using the realtime ion density and flow speed data obeying the equation (Banerjee et al. 2015(Banerjee et al. , 2019 As the solar wind flows towards the Earth, it interacts with the magnetosphere. Depending on the intensity and direction of IMF BZ, a part of the solar wind energy is injected into the magnetosphere through the cusp while the major part of the solar wind energy is deflected at the bow shock and flows across the Earth. A small fraction of this deflected energy penetrates into the magnetosphere. In this model, the factor K represents this small fractional value. Thus, all the elements of the lattice are credited with the energy K× dE at every initial stage. This alters the potential energy of each element as (Banerjee et al. 2015(Banerjee et al. , 2019 As the two-dimensional lattice is analogical to the terrestrial magnetosphere, the cusp width W C in the model is a set of selected elements including and surrounding the centre one, the element at i = n/2, j = n/2. The number of elements in the set, that is, the size of the cusp width W C is controlled by both the direction and intensity of IMF B Z , following the equations (Banerjee et al. 2015): Here B Zmax is the maximum value of northward B Z and K d is the associated proportionality factor. K d is a function of the direction of the IMF B Z . As the energy injection increases for southward B Z and decreases for northward B Z , the numerical value of K d is higher for the southward direction than that of the northward direction. During a solar wind-magnetosphere coupling, the solar wind energy dE is injected into the system through the cusp width W C following the relation (Uritsky et al. 2001;Banerjee et al. 2015Banerjee et al. , 2019, The threshold value of energy for each element is indicated as E TH , the excitation threshold. After the energy injection, if E(i, j) < E TH , the element is stable. If E(i, j) > E TH , the element is unstable and releases four units of energy to return to a stable state. The released energy is distributed among its four adjacent neighbours according to the following equations (Uritsky et al. 2001;Banerjee et al. 2015Banerjee et al. , 2019, A small value of energy E d is dissipated during the distribution process. The energy of the adjacent elements alters as (Uritsky et al. 2001;Banerjee et al. 2015Banerjee et al. , 2019) Figure 1 illustrates the energy transfer process in a two-dimensional lattice. The unstable elements [E (i, j) > E TH ] are marked with black shade, the stable elements [E(i, j) < E TH ] with grey shade, and the elements with zero or negligible energy with white shade. Initially, all the elements of the lattice have zero energy. The input energy is injected into the system through the cusp W C altering its energy above E TH , as shown in Figure 1(a). Thus, the centre element is in black shade. The cusp W C here consists of only one element, the element at the centre of the lattice (i = n/2, j = n/2). In the next Figure 1(b), the unstable element distributes four units of its energy to its four adjacent neighbours. Two of the neighbouring elements, marked in black shade became unstable, receiving the excess energy. In the final Figure 1(c), the two unstable elements further distribute their energy to return to stability. This way, the injected energy is distributed and redistributed throughout the lattice.
With analogy to the Earth, the system is spherical, meaning the elements belonging to the columns j = n and j = 1 are adjacent neighbours. If the energy of any element belonging to the column j = n has crossed the threshold, one unit of its surplus energy is distributed to its neighbour element in column j = 1 and vice versa. The dissipation of energy through the marginal grids is only applicable for the upper and lower marginal rows, not for the marginal left or marginal right columns. Open boundary condition has been considered for the marginal rows of the lattice. After a consecutive distribution and redistribution, when the released energy finally reaches the marginal upper [i = 1, j = (1 to n)] and lower [i = n, j = (1 to n)] grids, it transmits outside the system. Figure 2 displays the state of the elements of the lattice after a consecutive energy injection, distribution, and redistribution. All the elements have a considerable amount of energy while some of them become unstable. In Figure 2(a), in the marginal rows, two of the elements are unstable in both the upper and lower grids. In the next Figure 2(b), the unstable elements in the marginal rows dissipate their excess energy and a part of this energy is released outside the lattice, representing the magnetosphereionosphere energy transfer.
After the consecutive distribution and redistribution, the total internal accumulated energy of the lattice at any time t can be calculated as (Uritsky et al. 2001;Banerjee et al. 2015Banerjee et al. , 2019 where S ij = 1 if E > E TH and S ij = 0 if E < E TH The upper and the lower grids can be considered as the northern and southern polar cusps of the Earth while the energy transfer process is similar to magnetosphere-ionosphere energy coupling. In the actual measurement of the AE index, all the magnetometer stations are located in the northern hemisphere for recording the associated data. Following this instance, we too considered and based our calculations on the total excess energy dissipating only through the upper grid [i = 1, j = (1 to n)] of the lattice. We denote T R as the time delay between the initial energy injection to the model and the beginning of the energy transfer through the upper marginal grid. T R estimates the total time for the energy to reach the boundary after successive distribution and redistribution processes for each input injection. The total amount of released energy outside the lattice at any time t can be estimated as (Banerjee et al. 2019) The released energy generates two currents, I W and I E , the numerical equivalents of the auroral westward and eastward currents in the ionosphere, respectively. The auroral electrojet is estimated as the difference value of these two currents. For the calculation of electrojets, we considered the dissipated energy from any element belonging to the odd columns, that is, j = 1, 3, 5, . . ., (n-1) are contributing to the westward electrojet while the energy from any element belonging to the even columns, that is, j = 2, 4, 6, . . ., n are contributing to the eastward electrojet. Thus, mathematically The unstable element distributes four units of its energy to its four adjacent neighbours. Two of the neighbouring elements became unstable. (c) The two unstable elements further distribute their energy to return to stability. This way, the injected energy is distributed and redistributed throughout the lattice.
and E e t ð Þ¼ X E i; j ð Þ for i ¼ 1; j ¼ 2; 4; 6; . . . ; n (13) As the distribution process continues, more and more excess energy piles up outside the boundary regions of the lattice. The series E w (t) and E e (t) are the estimations of the total accumulated energy released in the ionosphere responsible for the westward and eastward electrojet currents, respectively. The released energy takes time to completely dissipate through the current system in the ionosphere and a part of this energy remains stored in the ionosphere. During the next magnetosphere-ionosphere energy transfer, this stored energy from the previous transfer acts as a base value and adds up with the newly released value of energy. Here, we introduce a parameter, namely K A which determines the remaining part of the released energy of the previous transfer, stored in the ionosphere. K A has a fractional value. The remaining energy in the westward region is whereas the remaining energy in the eastward region is Mathematically, the total westward energy at any time t is and the total eastward energy at any time t is These two energy components, E tw (t) and E te (t) then drive two currents in the opposite direction throughout the auroral region, namely the westward electrojet current and the eastward electrojet current, respectively.
The maximum westward electrojet current can be considered as and the maximum eastward electrojet current as The differential value of these two components is For further refinement, E N is processed by a filter and labelled as E A . Finally, the output time-series E A can be regarded as a numerical representation of the natural auroral electrojet index, AE. For our analysis, we used the hourly averaged AE index, solar wind ion density, flow speed, and B Z component of the interplanetary magnetic field (IMF) data from the year 1997 to the year 2007 of the 23 rd solar cycle. The dimension of the lattice is 50 × 50. Similar to our original study of the model (Banerjee et al. 2015), the numerical value of the various variables of the model are taken as K = 0.0025, local dissipation term E d = 0.05, K d = 0.5 for southward direction, and K d = 0.005 for northward direction. The input to the model is estimated using equation 1.

Data source
Here we used the hourly averaged AE index, solar wind ion density, flow speed, and B Z component of the interplanetary magnetic field (IMF) data from the year 1997 to the year 2007 of the 23 rd solar cycle as extracted from NASA/GSFC's OMNI data set through OMNIWeb. The OMNI data were obtained from the GSFC/SPDF OMNIWeb interface at http:// omniweb.gsfc.nasa.gov (King & Papitashvili 2005).

Result and discussions
The threshold excitation E TH is crucial for the central characteristics of a SOC system as it allows the existence of multiple metastable states across which the avalanches are carried out throughout the system. The potential energy E of any element in the lattice is analogical to the slope of an actual sandpile. An unstable element having energy E(i, j) > E TH releases four units of energy to return to stability. Thus, a minimum value of E TH = 5 is required to keep the potential energy of the element at a positive non-zero value and to avoid the total internal energy of the system reaching zero or negative energy states at any time. To analyse the dynamical behaviour of the system, the model is subjected to three different E TH values, E TH = 5, E TH = 7, and E TH = 9. After the energy injection and distribution, the total accumulated internal energy of the lattice, E Total , the total number of unstable states, S UN, and the total amount of released energy, E R at any time t can be calculated according to Equations 9, 10, and 11 respectively. Figure 3 illustrates the time vs. E Total plot for the three values of E TH for the year 2002. As the threshold value increases, the rate of energy distribution and dissipation decreases while the total internal energy of the system continues to increase gradually. It takes much more time for the system to achieve a meta-stable state, delaying the SOC dynamics and avalanches throughout the lattice to form properly. As the rate of distribution reduces, the elements in the marginal grids start to store and release energy to the outside of the lattice far more lately, thus delaying the energy transfer process representing the magnetosphere-ionosphere energy coupling in the model. We denote T R as the time delay between the initial energy injection to the model and the beginning of energy transfer through the upper marginal grid. It is observed from the analysis, that the value of T R is T R = 1434 hours for E TH = 5, T R = 1934 hours for E TH = 7 and T R = 2520 hours for E TH = 9. Thus, for further analysis of the system, we considered E TH = 5 as the value for the excitation threshold. Figure 5(b), the plot for the simulated timeseries E A , shows the estimation of E A initiating from the value of T R = 1434 hours for E TH = 5 as from this value of T R the marginal grid starts to release energy outside the system. Figure 4 demonstrates the plots for dE, E Total , S UN, and E R for the September-October period of the year 2002. As seen from the figure, dE has large values for the marked timeline. The injected energy piles up in the lattice increasing the number of unstable states S UN , the total energy of the lattice, E Total , reaches a critical point and then gradually returns to As the threshold value increases, the rate of energy distribution and dissipation decreases while the total internal energy of the system continues to increase gradually. It takes much more time for the system to achieve a meta-stable state. a metastable state by releasing a burst of avalanches. The avalanches take a time to reach the upper marginal grid through successive distribution and redistribution processes, thus E R shows high-energy values dissipating from the upper grid after a while.
After the injection of the energy followed by subsequent distribution and redistribution, the model output E A is generated estimating the released energy outside the lattice as discussed in the method section. The output series E A is the numerical representation of the natural AE index. The power spectral density (PSD) of the simulated output E A is calculated and plotted in a log-log graph. The plot demonstrates the characteristic 1/f β behaviour of the natural AE fluctuations along with the spectral break at f 0 . A detailed study of the plot revealed the value of breakpoint f 0 as f 0 = 0.050 mHz (5.5 hours), as shown in Figure 7(b). A power law is fitted separately in both the high (f > f 0 ) and low frequency (f < f 0 ) regions of the plot to determine the power-law coefficients (slope of the spectral response) β A and β B , respectively. Now keeping the other parameter constant, if the value of K A is varied in the range of 0.10 to 1.00, the β values also vary within ranges. As, for a particular value of K A , the β values of the simulated series closely match with that of the natural AE index, they are noted down along with the value of K A , as shown in Table 1. The β values of the real-time AE index are also displayed in Table 1 for a comparative study. Figure 5(a,b) are the time series of the natural AE index and the simulated E A series of the year 2002, respectively. For comparative purposes, a magnified portion of the natural AE index and that of the simulated E A series of the year 2002 are displayed in Figure 6(a, b), respectively. The real-time AE index is estimated as the difference value of AU and AL indices where the AU and AL indices are the direct measurements of the maximum eastward and westward electrojet currents. (Ahn et al., 2000b) studied the variation pattern of the yearly mean AL index and AU index for 20 years and suggested their absolute values are proportional to each other. The maximum is observed for the AU index in summer while for the AL index, it is in equinoctial months. Both the indices exhibit higher values in the descending phase of the solar cycle (Ahn et al., 2000b). In the current model, the simulated AE index is estimated following the same relation as the difference value between the maximum eastward and westward electrojets. In Figure 5, the May-June summer months of the year 2002 is the period for about the time, Time = 2800 hours -4200 hours. As can be seen from Figure 5(a), the values of the real-time AE index are in the higher ranges for these months. The simulated series in Figure 5(b) also demonstrates the same. Again, in Figure 5, the period of about the time, Time = 5800 hours -6600 hours marks the equinoctial month of September. Here, also the values of the simulated series of Figure 5(b) show higher values, similarly to the real-time AE index, as shown in Figure 5(a). Figure 7(a, b) illustrate the log-log plot of the PSD of the real-time AE index and the simulated E A series of the year 2002, respectively. It is observed from the plots that the simulated E A series exhibits the characteristic 1/f β behaviour of the natural AE fluctuations along with the spectral break at f 0 . As can be seen from Figure 7 Table 1. The power-law coefficients (slope of the spectral response) of the power spectral density associated with the real-time AE index and the simulated model output E A for all the years of the 23 rd solar cycle. The spectral break is at f 0 = 0.050 mHz (5.5 hours). β A denotes the value of the slope for f > f 0 and β B denotes the value of the slope for f < f 0 . The parameter K A is associated with the series E A and shows different values for different years.   slopes of the two frequency regions of the 1/f β power spectrum of the natural AE index had long been a study of keen interest. (Tsurutani et al. 1990) estimated the values of the slopes as β A = 2.2 and β B = 0.98 with a breakpoint at f 0 = 0.050 mHz (5.5 hours) for the hourly average AE data of the period of years 1971-1974. (Uritsky and Pudovkin 1998b found out the values as β A = 2.10 and β B = 0.95 with a breakpoint at f 0 = 0.055 mHz (5 hours) for the hourly average AE data of the period of years -1974. (Woodard et al. 2005 compared all the prominent studies (Tsurutani et al. 1990;Consolini et al. 1996;Uritsky and Pudovkin 1998b;Price and Newman 2001;Watkins 2002) investigating the slopes of the two spectral regions and breakpoint of natural AE fluctuations and concluded the values of slopes as β A = 2.4 ± 0.26 and β B = 1.0 ± 0.10. It is observed from Table 1 that for a particular value of the parameter K A , the β values of the simulated series E A are estimated as β A = 2.2-2.4 and β B = 0.9-1.0 with a spectral break at f 0 = 0.050 mHz (5.5 hours), the typical values associated with natural AE index as reported by all these previous works.
The parameter K A plays a significant role in estimating the simulated E A series from the total released energy. As observed from Table 1, for a particular year and a particular value of K A , the β values of the simulated output series E A are in the specified ranges of the same associated with the natural AE index. For the entire 23 rd solar cycle, the value of K A is in the range of K A = 0.68-0.82. It is observed from the result, that the transferred solar wind energy does not dissipate completely through the westward and eastward auroral electrojet currents in an instant, rather a significant part of it remains present in the ionosphere even in the time of the next magnetosphereionosphere energy transfer. The parameter K A has a fractional value. As can be seen in equations (14) and (15), K A is denoting the remaining fraction of the total accumulated released energy of the previous state, reserved in the ionosphere. The excess energy transferred from the magnetosphere to the ionosphere is then being added up with this base value and forms the westward and eastward currents that are the two key factors for measuring the AE index.
The magnetosphere-ionosphere energy transfer and the stored part of this energy in the ionosphere are primarily controlled by the solar wind injection into the magnetosphere. Again, the amount of injected solar wind energy in the magnetosphere varies with the intensity, and duration of solar windmagnetosphere coupling. Also, the solar cycle has a significant effect on energy injection. During solar storms, there is a large deposit of energy into the magnetosphere which changes the normal quiet time dynamics of the magneto-ionosphere system. Figure 4 illustrates such a case of a large amount of solar wind injection into the magnetosphere. As shown in Figure 4(a), the marked period in the figure is distinguished by the continuous injection of a large amount of energy dE into the magnetosphere for hours. Consequently, the total accumulated energy of the lattice, E Total , gradually starts to pile up in the system, altering its state of equilibrium. The value of the number of unstable states, S UN , also increases indicating the instability formation in the lattice. Figure 4(b, c) show the plots for E Total and S UN , respectively. As the large values of energy injection continue, E Total finally reaches a critical point and the pile collapses. The excess energy is released as a burst of avalanches inside the lattice restoring the stability of the system. After the successive distribution and redistribution process, the excess energy of the lattice is finally transferred outside its boundary region representing the magnetosphere-ionosphere energy transfer. As seen from Figure 4(d), the amount of released energy E R increases after a while as a consequence of the continuous injection of large values of dE in the magnetosphere. The released energy in the ionosphere is dissipated through the auroral electrojets, the fieldaligned currents flowing between the nightside magnetosphere and nightside ionosphere. But it takes time to completely dissipate through the current system in the ionosphere and a part of this energy remains stored in the ionosphere. The factor K A estimates the amount of this stored energy.
As observed from the result, for the entire 23 rd solar cycle, the value of K A is slightly varying, K A = 0.68-0.82. In the model, the parameter K A does not estimate the actual value of the stored energy, rather it indicates the percentage value of the transferred energy that remains stored in the ionosphere. As shown in Table  1, the value of K A for the year 2002 is K A = 0.73. This value indicates, at any current state t, 73% of the transferred energy of the previous state (t-1), remains stored in the ionosphere. Thus, the transmitted energy does not dissipate instantly after each magnetosphereionosphere energy transfer, rather a substantial percentage (73%) of this energy remains reserved in the ionosphere during the next magnetosphereionosphere energy transfer. If the solar wind energy injection into the magnetosphere increases at any time t, the amount of total released energy in the ionosphere also increases for this t. Since K A indicates the remaining percentage value of this released energy, the actual amount of the reserved energy in the ionosphere also increases as an effect of the large solar wind injection. Similarly, if the input injection is small, the amount of released energy decreases, decreasing the actual amount of reserved energy in the ionosphere. Thus, the actual amount of the reserved energy in the ionosphere varies with the variation in the intensity of injected solar wind in the magnetosphere and the parameter K A represents its percentage relationship with the total transferred energy. For the 23 rd cycle, the value of K A ranges between K A = 0.68-0.82, indicating that for each year, a substantial percentage of the transferred energy remains stored in the ionosphere.

Conclusions
AE index is a global and instantaneous measurement of the magnetic fluctuations in the Earth's polar region in response to an external perturbation. In this paper, we developed a numerical cellular automata model of Earth's magnetosphere based on the concept of self-organised criticality and sandpile dynamics to study the complex dynamics of the magnetosphere-ionosphere energy transfer process and AE fluctuations. Our model is a dissipative, dynamical n × n two-dimensional system of finite potential energy and open boundary conditions having the real-time values of solar ion density, flow speed, and IMF B Z as the input parameters. It is analogical to the Earth's magnetosphere while the upper and lower margins of the lattice can be considered as the north and south polar cusps of the Earth. The solar wind is a stream of energised plasma particles emitted from the outer atmosphere of the Sun. As the direction of IMF B Z is southward, a strong coupling occurs between the solar wind and terrestrial magnetosphere injecting a significant amount of solar wind energy into the geospace. Gradually, the energy piles up and reaches a self-organised critical condition after which adding up a small amount of energy into the pile can form a spatially localised magnetospheric instability. To maintain the equilibrium, the system redistributes itself and the excess energy is released as an outburst of avalanches of various sizes in the neighbouring regions. As long as there is local instability, the distribution process continues and successively spread over throughout the system, finally transmitting a large amount of energy into the ionosphere through the polar cusps. The transferred solar wind energy causes magnetic fluctuations in the auroral region and the AE index is the global measurement of the intensity of the fluctuations. In our proposed cellular automata model, the marginal grids of the lattice are equivalent to the polar cusps. The excess energy transferred through the upper grid of the lattice is measured and by some mathematical process a simulated time series has been derived which can be considered as a numerical representation of the real-time AE index.
The spectral response of the simulated output series E A follows a 1/f β power law, demonstrates a breakpoint at f 0 = 0.050 mHz (5.5 hours) having slopes β A = 2.2-2.4 for f > f 0 and β B = 0.9-1.0 for f < f 0 , the typical characteristics of natural AE index. It is observed that the parameter K A plays a significant role in the entire process of forming a proper E A time series estimated from the released energy. K A represents the percentage of the released energy from the previous magnetosphere-ionosphere energy transfer, which remains stored in the ionosphere. Its value varies in a small range of K A = 0.68-0.82 for the eleven years of the 23 rd solar cycle, indicating a substantial percentage of the transferred energy remains reserved in the ionosphere for each year.
The excitation threshold E TH is crucial in forming the SOC dynamics of the model. As an unstable element release four units of energy, for a small value of E TH , the total internal energy of the system can achieve zero or negative potential for some of the values of t. In contrast, a large value of E TH causes the piling up of solar wind energy in the system, reducing the rate of energy distribution as well as dissipation. The total internal energy of the system continues to increase gradually and the system takes a much larger time to achieve a metastable state and to form the SOC dynamics properly. Thus, a moderate value of E TH = 5 is optimum for the proposed model.
Overall, it can be concluded that our proposed model is a simple first-order avalanche model of the Earth's magnetosphere where the real-time solar parameters are the inputs. The model generates a simulated output series E A which shows statistical similarity to the real-time AE index. Also, the parameter K A varies over a small range of K A = 0.68-0.82 which suggests a high percentage of the transferred solar wind energy remains reserved in the ionosphere. In our previous work (Banerjee et al. 2015) we presented a SOC-based cellular automata model and focused on the nature of solar wind-magnetosphere energy transfer and its subsequent effect in the magnetosphere. In continuation of this project, the current work is a further refinement of that model to study the magnetosphere-ionosphere energy transfer process. Future work can be focused to develop a composite cellular automata model of the magnetosphere to study all the intricate characteristics of the solar wind-magnetosphere-ionosphere dynamics.
College, Kolkata, India for their encouragement and support for the present work. Finally, we would like to sincerely thank the anonymous reviewers for their most valuable comments and suggestions to improve the quality of this paper.

Disclosure statement
No potential conflict of interest was reported by the author(s).