Cellular automata model of magnetospheric-ionospheric coupling

Abstract. We propose a cellular automata model (CAM) to describe the substorm activity of the magnetospheric-ionospheric system. The state of each cell in the model is described by two numbers that correspond to the energy content in a region of the current sheet in the magnetospheric tail and to the conductivity of the ionospheric domain that is magnetically connected with this region. The driving force of the system is supposed to be provided by the solar wind that is convected along the two boundaries of the system. The energy flux inside is ensured by the penetration of the energy from the solar wind into the array of cells (magnetospheric tail) with a finite velocity. The third boundary (near to the Earth) is closed and the fourth boundary is opened, thereby modeling the flux far away from the tail. The energy dissipation in the system is quite similar to other CAM models, when the energy in a particular cell exceeds some pre-defined threshold, and the part of the energy excess is redistributed between the neighbouring cells. The second number attributed to each cell mimics ionospheric conductivity that can allow for a part of the energy to be shed on field-aligned currents. The feedback between "ionosphere" and "magnetospheric tail" is provided by the change in a part of the energy, which is redistributed in the tail when the threshold is surpassed. The control parameter of the model is the z-component of the interplanetary magnetic field (Bz IMF), "frozen" into the solar wind. To study the internal dynamics of the system at the beginning, this control parameter is taken to be constant. The dynamics of the system undergoes several bifurcations, when the constant varies from - 0.6 to - 6.0. The Bz IMF input results in the periodic transients (activation regions) and the inter-transient period decreases with the decrease of Bz. At the same time the onset of activations in the array shifts towards the "Earth". When the modulus of the Bz IMF exceeds some threshold value, the transition takes place from periodic to chaotic dynamics. In the second part of the work we have chosen as the source the real values of the z-component of the interplanetary magnetic field taken from satellite observations. We have shown that in this case the statistical properties of the transients reproduce the characteristic features observed by Lui et al. (2000). Key words. Magnetospheric physics (magnetosphere-ionosphere interactions) – Space plasma physics (nonlinear phenomena)


Introduction
The Earth's magnetosphere is immersed in the solar wind. It can be considered as an open dynamic system driven by external and internal sources. These sources are provided by the solar wind and ionosphere that can be treated as a part of the magnetospheric system. In general, considered as the dynamic system it is far from equilibrium, and one would assume that the dynamics of such a system can be quite complex. However, when the external forcing is stationary, such a system can possess non-equilibrium stationary states, as it was shown by Prigogin (1955) and Klimontovich (1995). If the system belongs to such a class, its dynamics can be characterized by a small number of parameters; in other words, the trajectory of the system can be placed on the low-dimensional manifold by means of choosing the appropriate set of variables and coordinates. To characterise this property of complex dynamical systems the notion of selforganization is conventionally used. Many different characteristics of the magnetospheric system have been studied in a number of papers using the dynamic system approach. Indeed, it was unambiguously shown that the magnetosphericionospheric system interacting with the solar wind as an external driving force has a low correlation dimension during substorm activity (Vörös, 1991;Takalo et al., 1993;Klimas et al., 1996;Sharma et al., 1997). For a low-dimensional chaotic system the correlation dimension is an estimation of Bz IMF Fig. 1. Scheme of the "magnetospheric" array and boundary conditions.
the dimension of the trajectory attractor. However, in the case of the colored noise produced by a stochastic system the low correlation dimension is also possible (Osborne and Provenzale, 1989). Then, the correlation dimension is only a measure of the fractal dimension of the system trajectory. Another important question related to the dynamics of the magnetosphere is does it behave as a self-organized critical system?
The notion of self-organized criticality (SOC) was proposed by Bak et al. (1988). It was introduced as a universal feature of complex systems in which the interaction has quite a complex character and there are no characteristic scales. All scales interact with each other, and the behaviour of the system on different spatial scales is self-similar. In such a case the correlation length is infinite. In real systems it is supposed to be of the same order or larger than the characteristic length of the system itself. Such systems can be also characterized by a small number of parameters, and many characteristics of the systems of such a class were studied using "sand-pile" type models. An excellent review on this topic was written by Jensen (1998). This approach has found multiple applications, from forest fires to earthquakes. However, for some applications the relevance between such a simple sand-pile model and a description of physical processes is not established (For more details, see the reviews by Vespignani and Zapperi, 1998;and by Jensen, 1998). For magnetospheric-ionospheric system (MIS) there are several indications that it can behave similarly to SOC systems (Consolini, 1997). Mainly, these indications are the power law spectra of fluctuations (Vörös, 1991;Milovanov et al., 1996;Sharma, 1997;Uritsky and Pudovkin, 1998;Takalo et al., 1999;Uritsky et al., 2002). Moreover, Chang (1992Chang ( , 1998 has pointed out that a SOC system may exhibit a low-dimensional behaviour. The simplest avalanche model of magnetospheric dynamics based on SOC paradigm was presented by Chapman et al. (1998).
However, the problem of whether the magnetosphericionospheric system belongs to the SOC-type systems or not, is not so simple. First of all, the driving force for the con-ventional SOC models is supposed to be random and acting in an arbitrarily chosen spatial part of the system. For the MIS the external forcing is ensured by the solar wind. Fluctuations of the solar wind magnetic field and velocity are turbulent, and their statistical properties are far from Gaussian (Burlaga, 1991;Tu and Marsch, 1995). So there is a question of whether the dynamics of the system is determined by its internal rules or whether it is imposed and mainly reflects the statistical properties of the external driver. For instance, according to Klimas et al. (1998), the dynamic behaviour of the D st index is linearly related with certain characteristics of the parameters of the solar wind. Recently, the low-dimensional global dynamics of the magnetosphere and the multi-scale features of the solar wind-magnetosphere coupling during substorms have been combined in the non-equilibrium phase transition model (Sitnov et al., 2000(Sitnov et al., , 2001.
Another important feature of the MIS related to the solar wind is that the driver operates at the systems' boundaries; thus, the internal dynamics can strongly depend on the way of penetration of the energy from the source to the system. Thus, the question of whether the non-equilibrium stationary states can exist in the MIS seem to be more academic than practical. However, the complex dynamic system can be properly identified by making use of its quasi-periodic trajectories while varying its control parameters. Therefore, the study of such classes of trajectories and bifurcations of the dynamic system allows one to better establish the relevance between the model and real system, though it does not allow one to describe properly real observations.
In previous papers Kozelova, 2000, 2002a) we considered the classification of spontaneous and stimulated transients in the SOC-like system modelled by cellular automata controlled by B z IMF. Developing this model, in papers Kozelova, 2002b, 2002c) we included in the model the positive feedback as an analogy of the feedback between the Earth's magnetosphere and ionosphere, which is activated during an explosive phase of substorm and does not work in quiet time. Such a model has allowed us to explain the obtained (Lui et al., 2000) differences of distribution functions of power dissipated by auroral spots for quiet time and during substorms.
Hereby, we make the following step from the qualitative model of a sand pile to the physically more real model. Saving internal local connections from our previous model (see detailed description below), we organize an external influence on the model as a flow of the magnetospheric tail by the solar wind. We consider that the effect from solar wind penetrates inside the system with finite velocity. In addition, now the internal rules in the model are strictly deterministic.
The control parameter of the model is the B z IMF "frozen" into the solar wind. To study the internal dynamics of the system, it is taken to be constant. We have found that the dynamics of the system undergoes several bifurcations when the constant varies from −0.6 to −6.0. At the control by experimental sequence of B z IMF values, the statistical distributions of the characteristics of transients in the model resemble the features of similar distributions observed by Lui et al. (2000) for auroral spots.

Model
The model described below develops the approach used by Kozelova, 2000, 2002a, b, c). The principal element of the model is a sand-pile type cellular automata with special boundary conditions. The current sheet of the magnetospheric tail is represented as a rectangular array that consists of 50 × 100 cells (see Fig. 1). One short boundary is closed, it corresponds to the boundary of the current sheet closest to the Earth; the opposite boundary (tailward) is opened, and the energy and the magnetic field can be transported through it outside the system towards the interplanetary space. Two other boundaries represent the boundaries of the magnetosphere flown around by the solar wind that carries along them the interplanetary magnetic field. This field can penetrate through these last boundaries into the current sheet. The source is described as the convection of the external magnetic field with finite velocity, which is supposed to correspond to the velocity of waves that carry these perturbations inside. This means that the process considered corresponds to the transport of the magnetic field flux into the system by magneto-hydrodynamic (MHD) waves. The "state" of a cell with coordinates (i, j ) at each moment of time t is characterized by two numbers. One is E t (i, j ), which is considered as the energy stored by the cell at a moment t, another number attributed to each cell, and C t (i, j ) is assumed to describe the conductivity in the ionospheric region associated with the same magnetic tube as a cell of the "current sheet" in the "magnetospheric tail".
As is usual in cellular automata models the time is discrete, but it will be presented in conditional units that will be called seconds and minutes. At every step in time in each cell of the array the small portion of energy dE(i, j ) is added. As long as the value of the stored energy E t (i, j ) does not exceed some critical level E max (i, j ), the cell remains steady.
In contrast to the models discussed in Kozelov and Kozelova (2002a, b), now we suppose that the threshold value in each cell is individual and depends on an external control parameter which influences the long boundaries of the rectangular array. The values at the boundaries of the array E max (1, j ) and E max (50, j ) were determined by the value of (−B z ) IMF and were shifted along the boundaries with a velocity of 1 cell/min. The velocity of propagation of the disturbance inside of the array was also assumed to be 1 cell/min, and the value of the disturbance decreased proportionally to the distance from the boundary.
When the threshold level E max (i, j ) is exceeded the cell passes into an active state, and a certain part of the stored energy, E = E t (i, j ) − E min (i, j ), is distributed between four adjacent cells: (1) Here, the choice of the coefficients reflects the nonsymmetry in the direction from the closed to opposite boundary. After that, the energy of the adjacent cells can exceed E max . Then these cells at the following time-step transmit, in turn, the energy to their neighbours. It is supposed that the speed of internal transients in a current sheet is higher than the speed of propagation of an external disturbance. In our calculations 1 time-step for internal processes of reallocation in the system was 10 s.
Observational studies of magnetospheric activity suggest that the magnetospheric-ionospheric coupling plays a critical role in the physical processes leading up to a substorm onset. The local redistribution of energy in the magnetosphere causes a local change in the conductivity of the ionosphere in the same magnetic tube: the particles, diffused by pitchangle, are precipitated in the loss-cone along the magnetic field, and ionise atmospheric gases. Let us take into account that the ionosphere conductivity is proportional to the electron density n e , and the differential equation for electron density in a simple "thin" ionosphere is where Q is ionisation rate, α eff is an effective recombination coefficient. The power index is p = 2 for the E-region ionosphere, which is the most effective for disturbed conditions (bright auroral forms). For the F-region p ≈ 1, and this case corresponds to quiet conditions (soft diffuse precipitation). In our model we want to describe the transition between these situations; therefore, we will use p = 1. Then for conductivity in discrete form we have a rule: Here, (1−a) means a "recombination coefficient"; therefore, conductivity of the ionospheric part of a cell depends on the cell history. Selection of a = 0.2 gives us the characteristic time for recombination of τ ∼ 5 s, because during one time step (10 s), without source term b, the conductivity decreases to a 0.2 part of initial value. Taking into account that Q ∼ E, the second term is determined as In turn, if the ionospheric conductivity exceeds some level C max , then the conditions of energy redistribution in the current sheet are changed (at the expense of the formation of field-aligned currents). In the model we shall take into account this influence as the relation of E min (i, j ) from C t (i, j ). This E min (i, j ) value determines which part of the energy may be reallocated in an active cell with the following time step: Em a x j i Fig. 2. Scheme of the local positive feedback between "magnetospheric" and "ionospheric" parts of a cell: E t (i, j )-stored energy, C t (i, j )-ionospheric conductivity.
Here, k < 1(k = 0.75) and C max = 5 are parameters. The positive feedback arising between magnetospheric and ionospheric parts of a cell is shown schematically in Fig. 2. Each cell (except the boundary cells) has four adjacent cells; therefore, the occurrence of the chain response of energy transmission is possible, and it proceeds as long as there are no active cells in the array. Such a chain response is usually named an "avalanche" by analogy, with the formation of a slope for a sand pile (Bak et al., 1988). In our calculations we consider that for one time-step each cell can change its state only once, so the avalanche may have different durations ("running" avalanches). In addition in the considered model there is no additional internal source of stochasticity (random variable); therefore, in this sense the model is strictly deterministic.
One can see that rules defined by Eqs.
(1)-(4) have several numerical parameters. The values of these internal parameters reflect our view on the magnetosphere-ionosphere processes. However, their change in rather broad limits does not change qualitatively the considered dynamics of the model system.

Model dynamics for the constant external parameter
The dynamics of the model system at the varied external control parameter is rather complex. Therefore, for detailed study of the internal dynamics of the system we shall consider a more simple situation, when the control parameter does not vary during a long time. It is found that the following modes are possible in the model at the constant external parameter.
1. At B z > 0 the energy in the system is not accumulated, the processes of energy redistribution in the "magneto-spheric" part do not occur, the conductivity of the "ionosphere" is equal to zero.
2. At 0 > B z > −0.6 the energy in a system can be accumulated, but processes of energy redistribution in the "magnetospheric" part balance the excess of energy. The dynamical stationary state is established. However, the activity in the "magnetospheric" part is not sufficient for the excess of the threshold of ionospheric activation. Therefore, the positive "magnetospheric-ionospheric" feedback is not working.
3. At −0.6 > B z > −4.0 in some moments the threshold of ionospheric activation is exceeded, after which the redistribution in the "magnetospheric" part becomes more intensive, that leads to the increase of energy which flows out, then the activity is falling, etc. This is a mode of (pseudo-) periodical generation.
The origin and motion of transients may be tracked in keograms that are constructed by the method used in the analysis of television data of aurora registration. For this purpose, from each frame (in our case this is the array of cell state) one row (or column) is selected. A sequence of such rows (columns) for a series of instants gives us a keogram. Several keograms, constructed for the analysis of dynamics of transients in our model at B z = −3 nT, are shown in Fig. 3. The rows and columns in the array used for construction of the keograms are shown in the left panel. The grey scale shows the energy redistributed by cells, for which the threshold of activation of feedback is exceeded. In this sense the keograms are the analogy of the keograms of polar auroras.
In the upper keogram of Fig. 3 one can see that the transient starts at some distance from the closed boundary of the array, and extends and moves in the direction opposite to the closed boundary of the array. It is possible to obtain two obvious characteristics of the periodic mode observed in the system: the mean period and distance of the onset point from the closed boundary (see Fig. 3). Figure 4 shows the dependence of these characteristics on the parameter B z . The smaller periods between transients (activations) correspond to a smaller B z , in addition to the onset of activations approaching the closed boundary (to the "Earth").
The position of the onset point of transient depending on B z in the range of −0.6 > B z > −6.0 may be approximated by the relation where the square brackets mean the integer part of the number, B z is expressed in nT. The discrepancy of approximation does not exceed ±1.
4. At B z < −4.0 the periodicity of generation is obviously broken, which leads to a sharp decrease in the average  period at −5.0 > B z > −4.0. The numerical estimation of the maximum Lyapunov exponent by the method of Rosenstein et al. (1993) has shown that it is positive. This is a chaotic mode (Fig. 5), passing at B z < −10.0 in a "turbulent" mode, at which the feedback will be activated at one time step from energy redistribution from one cell.
Thus, the dynamics of the system at the constant control parameter depends strongly on the value of the control parameter. It is possible to consider the series of transitions from the absence of transients to a periodic and, furthermore, to a chaotic mode of generation as a typical series of bifurcations. It is obvious that the collective, large-scale tran-sients in the model are self-organized; however, they are not scale-free (especially for periodic mode). This is not a selforganized critical state.

Statistics for driving by B z IMF
First, we shall show that the dynamics of the suggested model really have common properties with the dynamics of the magnetospheric-ionospheric system. For this purpose we compare the statistical properties of transients in the model with properties of auroral blobs, which were obtained by Lui et al. (2000). We used the B z values measured by the IMP-8 satellite during 1974 as a control parameter of the model. To obtain the statistical characteristics, the array of the stored energy was checked for the presence of active cells at each time step. Let us note that the cell is considered active, if the threshold value E max is exceeded. The set of active cells was broken into linked clusters. Each cluster was considered as a separate transient. The size (square) and distributed power of each cluster was determined. The distributed power was determined as the sum of energy E distributed in all cells of the given cluster (transient) at the given time step.
The ionosphere activity was used for the relation of each time step to perturbed (substorm) or quiet conditions. If, in the system, there were cells with conductivity above the critical level C max , this moment was considered to be perturbed. If, for all cells of the system, the critical level of conductivity was not exceeded, this moment was considered quiet. As in Lui et al. (2000), we sort the distributions into five bins per decade. The total number of transients is ∼2 × 10 7 for quiet time and ∼10 7 for the perturbed time.
The distributions of the size (square) and power of transients for quiet and perturbed time are shown in Fig. 6. One can see that in all cases in the distributions it is possible to select a region of the power law falling with slope ≈ 1. In addition, the distributions for the perturbed time differ from similar distributions for the quiet time by the presence of a maximum at large values of the size and power. Thus, the distributions of the transient parameters in our model have a number of common distinctive features with the distributions obtained by Lui et al. (2000) for auroral blobs.

Discussion
In spite of the fact that the long periods of constant B z are rather rare, there are some indirect experimental confirmations of relations obtained in the previous section. In particular, Zverev et al. (1979) found the minimum latitude, up to which the auroral oval is spread, as a function of B z IMF. The oval is descended at lower latitudes with the increase of |B z |. This is qualitatively in agreement with Eq. (5).
The transients considered in the previous section are spontaneous (see further discussion). For spontaneous substorms (Dmitrieva and Sergeev, 1983) it was shown that the duration of the preliminary phase of the substorm decreased with the increase of |B z |. This is also qualitatively in agreement with the relation from Fig. 4b.
The periods of constant B z value are usually connected to periods of steady magnetospheric convection, during which substorms are not observed or are extremely rare (Sergeev et al., 1996). The steady magnetospheric convection events arise at rather large |B z | values (B z < −3); the active mesoscale transients are observed in the magnetosphere, but large-scale transients (substorms) do not develop. In our model such a condition has no direct analogy because we have no real convection. However, the chaotic mode in the model has some qualitatively similar features. In reality, it arises at B z values, when instead, large-scale periodic activations appear like more small-scale chaotic ones.
Note that the main goal of this paper is to study the influence of the external parameter (B z IMF) on the model dynamics. However, the characteristics of the pseudoperiodical regime also depend on internal parameters. Now we want to discuss once again the spontaneous and stimulated (triggered) transients. In the presented model at the constant negative B z , all transients obviously are spontaneous, as they are not triggered by any changes in external conditions. As we noted in Kozelova (2000, 2002a) the beginning of such transients is simply the corollary of the finite energetic volume of the system. Varied values of B z can lead to the condition, then the accumulated energy in some cell becomes higher than the threshold value because of the downturn of this threshold. It is natural to name such transients as triggered ones. However, the capability of trigger detection, even within the framework of our simple model, becomes problematic without detailed state information of the system near the position of transient onset. The types of transients and transient triggering in the model will be considered in more detail in a future paper.

Conclusions
An new cellular automata model as an analogy of the dynamical magnetospheric-ionospheric system related to the substorm activity is presented. Saving the internal local connections from our previous model, we suppose that the driving force of the system is provided by the solar wind that is convected along the boundaries of the system. The energy flux inside is ensured by the penetration of the energy from the solar wind into the system with a finite velocity. In addition in the model, we avoid using an additional internal source of stochasticity (random variable); therefore, in this sense the model is strictly deterministic.
The control parameter of the model is the z-component of the interplanetary magnetic field (B z IMF) "frozen" into the solar wind. At a constant control parameter the internal dynamics of the model system is determined by the parameter value. The collective dynamics of the system leads to largescale transients and undergoes several bifurcations when the constant varies from −6.0 to 0. The periodic transients (activation) occur in the system at values of −4.0 < B z < −0.6. Smaller periods correspond to smaller B z ; in addition, the position of activation onset in the array is shifting to the "Earth" with the falling B z value. When the modulus of the B z IMF exceeds some threshold value, the transition takes place from periodic to chaotic dynamics.
It was found that for the model driven by the experimental sequence of B z IMF, the probability distributions of the size and power of transients look like the distributions obtained by Lui et al. (2000) for auroral blobs. All distributions have a power law segment at small values of size and power. During the periods of the active feedback (substorm time for au-roral blobs), the distributions have a characteristic maximum at large values of size and power. The distributions for the moments, when the feedback in all the cells of the system is not active (quiet time), have no maximum at large values of size and power. Thus, the transients in the magnetosphereionosphere system have analogies within the framework of the deterministic model with complicated external driving.