# Self-Organized Criticality in a Spherically Closed Cellular Automaton:

Modeling Soft Gamma Repeater Bursts Driven by Magnetic Reconnection

###### Abstract

A new cellular automaton (CA) model is presented for the self-organized criticality (SOC) in recurrent bursts of soft gamma repeaters (SGRs), which are interpreted as avalanches of reconnection in the magnetosphere of neutron stars. The nodes of a regular dodecahedron and a truncated icosahedron are adopted as spherically closed grids enclosing a neutron star. It is found that the system enters the SOC state if there are sites where the expectation value of the added perturbation is nonzero. The energy distributions of SOC avalanches in CA simulations are described by a power law with a cutoff, which is consistent with the observations of SGR 180620 and SGR 1900+14. The power-law index is not universal and depends on the amplitude of the perturbation. This result shows that the SOC of SGRs can be illustrated not only by the crust quake model but also by the magnetic reconnection model.

###### pacs:

05.65.+b, 97.60.Jd, 45.90.+t, 94.30.cp^{†}

^{†}preprint:

## I Introduction

Self-organized criticality (SOC) proposed by Bak et al. bak87; bak88 has revealed a wide range of mechanisms of complex behavior occurring in nature. According to the concept of SOC, a nonequilibrium open system evolves spontaneously into a critical state characterized by a power-law distribution of its physical quantity. One of the most well-known examples of an SOC system is earthquakes, where the energy supplied by plate motion is dissipated in the crust. The SOC of earthquakes has been illustrated with a sand-pile model, and an empirical power-law relation between the size and frequency of seismic events has been demonstrated bak89; ofc92; a related model for the propagation of brittle failure had previously shown power laws katz86.

SOC models fit some properties of astrophysical phenomena, for instance, solar flares asch14. The recurrent bursts of soft gamma repeaters (SGRs) have event energy distributions well fitted with a power law gog99; gog00; prka12. It is currently thought that SGRs are associated with ultrastrongly magnetized neutron stars (10 G). They generally undergo the recurrent emission of soft gamma rays with a short duration (0.1 s). These bursts have been suggested to be due to neutron star crust fractures driven by the stress of an evolving magnetic field thdu95, which are called starquakes. Recently, however, Link link13 pointed out that it is difficult to reproduce the typically observed rise time of emission (10 ms) if the energy is deposited deep in the crust or deeper. The trapped seismic energy takes seconds to minutes to reach the stellar surface and cause burst emissions. As a corollary, the energy should be released not inside the star but in the magnetosphere katz82; lyut03; lyut06; kbl07; gillheyl10.

In this paper, we present a new cellular automaton (CA) model that mimics SOC avalanches in SGRs. We first follow the CA model for solar flares proposed by Lu and Hamilton luha91, which is expressed as discretized magnetohydrodynamic (MHD) equations isli98. In their model, solar flares are interpreted as avalanches of many small magnetic reconnections. Similarly, we assume that the energy is released in the magnetosphere of neutron stars through MHD instabilities, such as the tearing mode lyut03; kbl07. In the original model by Lu and Hamilton, a three-dimensional Cartesian coordinate grid of points is used. However, a neutron star may be sufficiently small for the size of avalanches to be comparable to the system size. In fact, since the energy of a magnetic field in a region of volume is , the length scale is cm with G and erg, which is the lowest observed energy of short bursts gog99; gog00. Note that the radius of neutron stars is typically cm. Furthermore, the perturbation rule in the original model by Lu and Hamilton tends to increase the magnetic field in a certain direction. Here, we use the grid with a closed geometry and modify the perturbation rule to satisfy a conservation law for magnetic field.

## Ii Cellular automaton model

### ii.1 grids

We introduce “spherically” closed grids for our CA model to map the surface enclosing a neutron star. To avoid grid anisotropy, we require that (i) the shape is nearly spherical, (ii) all edges are of equal length and (iii) all nodes are equivalent. Among the polyhedra which satisfy these conditions, the nodes of a regular dodecahedron and a truncated icosahedron (the shapes of a soccer ball and fullerene C), which have and 60 nodes, respectively, are adopted as our grids. They have a closed geometry, in marked contrast to CA models for SOC systems studied so far. Incidentally, Schein and Gayed schgay14 recently presented the techniques to construct the icosahedral Goldberg polyhedra, which are nearly spherical grids with more nodes. These grids satisfy the conditions (i) and (ii) but not (iii). Although the investigation with more nodes is interesting, we defer it to future work.

In the CA simulation, we assign the values of for the -th site of the grids, which represent the deviation of the radial component of the magnetic field from the unperturbed background configuration. We set for all sites at the initial time.

### ii.2 step 1: perturbations

As the first step of our CA rule, we add perturbations to . This perturbation corresponds to the bending of magnetic field line outward because is the deviation of the radial component. Since the CA simulation is carried out on the closed surface () in our model, the magnetic field should be perturbed to satisfy a conservation law,

(1) |

where is the outward-pointing unit normal vector of the surface . Therefore, we choose two neighboring sites and add a positive perturbation () to one and a negative perturbation () to the other. The perturbation, , is given as

(2) |

where is a random number obeying a Gaussian distribution with zero mean and unit variance. The amplitude of the perturbation, , is a model parameter as discussed later.

The sites where the perturbation is added are chosen as follows. In our model, we assume the bending of magnetic field line outward similarly to the rise of the magnetic field line on the surface of the Sun. Then positive and negative perturbations appear at the footpoints of the magnetic field line. For simplicity, we deal with two cases for the unperturbed background configuration, namely, poloidal and toroidal magnetic fields, and randomly choose one of the edges considered to be along the background field lines. The perturbations are added to the sites corresponding to both ends and the polarity of the magnetic field line is reflected in the sign of the perturbations. This procedure is illustrated in Figure 1.

The edges along the background field lines are determined as below. At first, for poloidal case, we choose edges connecting nodes with the different latitude. The number of such edges is 20 for a regular dodecahedron with 30 edges and 60 for a truncated icosahedron with 90 edges. Then the polarity is assigned for the end points of each edge reflecting their latitudes. Next, for toroidal case, we choose the edges to construct the circuits. In this process, the edges connecting nodes with the same longitude are not chosen and the number of edges considered is set to be same with the poloidal case. For a regular dodecahedron, circuits with 5 edges, 10 edges and 5 edges are constructed as shown in the bottom right of Figure 1. For a truncated icosahedron, circuits with 5 edges, 15 edges, 20 edges, 15 edges and 5 edges are constructed.

### ii.3 step 2: reconnections

The second step of our CA rule is reconnection. We define the magnetic field stress, which is the difference between the local magnetic field and the average of its three nearest neighbors , as

(3) |

The site is unstable to reconnection when the absolute value of its magnetic stress exceeds some critical value luha91; dimit11,

(4) |

This criterion is reasonable because a large magnetic stresses favors magnetic reconnection dimit11. If a reconnection instability occurs, the magnetic field stress is canceled as follows:

(5) |

resulting in . If the nearby sites become unstable owing to the reconfiguration expressed by (5), additional reconnection events occur, leading to a large avalanche in some cases.

### ii.4 step 3: event energy

When all instabilities in the grid have been relaxed, we return to the first step of adding perturbations. The energy released during the avalanche is defined as

(6) |

where the superscripts (0) and (1) denote the values before and after the avalanche, respectively.

Note the dimensionless quantities in our CA model. The threshold for the reconnection instability is fixed at in this study. Since the cases with same values of are equivalent, we only investigate the dependence on . If a reconnection event occurs at a single site , the released energy is . Therefore, the minimum released energy in our CA model is for .

## Iii Results

### iii.1 poloidal versus toroidal

Now we move on to the results of our CA simulations. In Figure 2, cumulative energy distributions of avalanches, , are plotted for the case of . As can be seen, the models with perturbations under poloidal and toroidal fields have considerably different profiles. In the models with poloidal perturbations, large avalanches with the size of the system occur, which is one of the characteristic features of SOC. In contrast, the models with toroidal perturbations do not reach the SOC state. There are very few bursts with energy while small bursts occur accidentally. In the toroidal models, the time average of the random perturbations is zero at each site (see Figure 1). This result is consistent with that of Lu and Hamilton luha91, where SOC was not found with the random perturbation being symmetric about zero. However, in the poloidal models, the system reaches the critical state while the spatial average of the perturbations is zero. In this case, the time average of each node depends on the latitude. It is positive for one hemisphere but negative for the other hemisphere. This is because, as already stated, the polarity is assigned reflecting the latitude.

Here, we investigate the mixed case of poloidal and toroidal fields. In the step 1 of the CA rule, we choose the poloidal perturbation or toroidal perturbation randomly and the probability of poloidal is set to . Thus the purely poloidal and toroidal cases correspond to and , respectively. The dependence of is shown in Figure 3 for the cases of and the truncated icosahedron model with nodes. As already mentioned, in our CA model, bursts with happen even if the system does not reach the SOC state. Thus the distributions with may include a contamination and we pay attention to the events with . For the case with , where the poloidal and toroidal fields are regarded as comparable, large avalanches are still observed. On the other hand, for the dominantly toroidal case (), it is hard to find the SOC feature.

### iii.2 parameter dependence

Hereafter, we focus on the purely poloidal models. The cumulative distributions appear to decrease rapidly at high energies owing to finite-size effects. In fact, the cutoff energy of the truncated icosahedron model with nodes is higher than that of the regular dodecahedron model with . A cutoff feature was also found in the observations of SGR 180620 and SGR 1900+14 prka12. While the statistical significance of the observed cutoff is not high, bursts with the highest energy may correspond to avalanches enclosing the central neutron star. On the other hand, the distributions of the truncated icosahedron model and regular dodecahedron model are similar at low energies, where finite-size effects are negligible.

We show the dependence on the amplitude of perturbation in Figure 4, where the cumulative distributions are plotted for the models with the truncated icosahedron () grid and poloidal perturbations. We find that the distributions with are well fitted with power laws. The power-law index, , of the differential distributions is related to the power-law fit of the cumulative distributions as . According to the observations of SGR 180620 and SGR 1900+14, the energy distributions of bursts have been fitted with -1.76 gog99; gog00 and prka12. In our CA simulation, depends on the value of and it is 1.44, 1.57 and 1.74 for , 0.2 and 0.4, respectively. Here, Kolmogorov-Smirnov probabilities of the fitting in the range of are 0.999, 0.999 and 0.964 for , 0.2 and 0.4, respectively. On the other hand, in the range of , Kolmogorov-Smirnov probabilities are 0.000 for all models due to the cutoff.

From Figure 5, we can see that the power-law index increases with the amplitude of perturbation . The most natural interpretation for this is that, for large , the perturbation of one site can grow exclusively and the site tends to become unstable before accumulating perturbations of other sites. As a result, the number of large avalanches, i.e. with a large energy, is reduced. The lack of universality of the power-law index is an interesting feature that was also observed in the model of Olami et al. ofc92.

## Iv Conclusion and discussion

In conclusion, we have shown the SOC behavior of a new CA model with spherically closed grids so as to demonstrate the recurrent bursts of SGRs. We adopted the nodes of a regular dodecahedron and a truncated icosahedron as our CA grids. Perturbations were added to satisfy a conservation law (1) and the sign was determined from the polarity of the unperturbed background magnetic field line. For the configuration of the unperturbed field, both poloidal and toroidal cases were considered. We found that the SOC state is reached only for the poloidal case owing to the existence of sites where the expectation value of the added perturbation is nonzero.

The cumulative burst energy distribution has a cutoff at high energies owing to finite-size effects. Similarly, the cutoff found in the observations of SGRs may be attributed to the compactness of the magnetosphere of neutron stars. In addition to the recurrent bursts, giant flares with enormous energy (10-10 erg) and long duration (100 s) are occasionally observed from SGRs maze79; hurl99; hurl05; palm05; tera05. The energy release in the magnetosphere of neutron stars was proposed early on as an explanation of a giant flare katz82. Taking the large difference in emission energies into account, it is difficult for our CA model to describe the giant flares, which might be produced by global field reconfiguration masa10.

The power-law index, , found in our CA model depends on the amplitude of perturbation, , and ranges from 1.2 to 1.8. This value is consistent with the observations. According to Prieskorn and Kaaret prka12, there is no difference between the power-law indices observed in high- and low-burst-rate regimes. This result may appear to contradict our CA model with dependence, but we found that cases with the same values of are equivalent. If the amplitude of perturbation and the criterion for reconnection depend on the unperturbed background magnetic field likewise (for instance, they are linearly related), and may be almost insensitive to the environment. In any case, further investigation of the relation between the power-law index and burst rate is important.

Comparing with the original model by Lu and Hamilton luha91, the grid and the perturbation rule for magnetic field are modified in this study. On the other hand, the CA rule of reconnection is identical with that in Lu and Hamilton model while the number of nearest neighbors is different. This CA model relates to the induction equation,

(7) |

where and are the plasma velocity and the resistivity, respectively isli98; dimit11. Note that the Laplacian of the magnetic field corresponds to the right hand side of equation (3) while the sign is opposite. Thus, if the reconnection occurs, the variation of is proportional to as in equation (5). On the other hand, the perturbation of this CA model is interpreted as the convective term . As already mentioned, this CA model was originally proposed for solar flares. In contrast, we have applied it to recurrent bursts of SGRs. The modifications for SGRs, such as the relativistic forms of the MHD equations and QED effects, will be interesting issues.

In this paper, we have shown that the SOC of SGRs can be illustrated not only by the crust quake model but also by the magnetic reconnection model. To mimic magnetic reconnections, other CA rules have been examined in the context of solar flares lu93; geo95; vass98; mocha10. Incidentally, we have revised the CA rule for the redistribution from equation (5) to that proposed by Lu et al. lu93, but the qualitative features remained unchanged. It would also be interesting to study the distributions of other SOC parameters (avalanche size, waiting time and so forth) asch14; mcin02, which will be presented elsewhere.