Fast radio bursts trigger aftershocks resembling earthquakes, but not solar flares

The production mechanism of repeating fast radio bursts (FRBs) is still a mystery, and correlations between burst occurrence times and energies may provide important clues to elucidate it. While time correlation studies of FRBs have been mainly performed using wait time distributions, here we report the results of a correlation function analysis of repeating FRBs in the two-dimensional space of time and energy. We analyze nearly 7,000 bursts reported in the literature for the three most active sources of FRB 20121102A, 20201124A, and 20220912A, and find the following characteristics that are universal in the three sources. A clear power-law signal of the correlation function is seen, extending to the typical burst duration ($\sim$ 10 msec) toward shorter time intervals ($\Delta t)$. The correlation function indicates that every single burst has about a 10-60% chance of producing an aftershock at a rate decaying by a power-law as $\propto (\Delta t)^{-p}$ with $p =$ 1.5-2.5, like the Omori-Utsu law of earthquakes. The correlated aftershock rate is stable regardless of source activity changes, and there is no correlation between emitted energy and $\Delta t$. We demonstrate that all these properties are quantitatively common to earthquakes, but different from solar flares in many aspects, by applying the same analysis method for the data on these phenomena. These results suggest that repeater FRBs are a phenomenon in which energy stored in rigid neutron star crusts is released by seismic activity. This may provide a new opportunity for future studies to explore the physical properties of the neutron star crust.


INTRODUCTION
Fast radio bursts (FRBs) are extragalactic transient phenomena that shine in radio wavelengths for short durations lasting only 1-10 milliseconds (see Platts et al. 2019;Zhang 2020;Petroff et al. 2022, for reviews). Some FRB sources are known to produce many bursts repeatedly. Repeaters are thought to be neutron stars, but the causes of bursts and the radiation mechanism are not well understood.
More than several thousand bursts have already been detected from a few repeaters, and a detailed statistical analysis of their occurrence times and energies may provide some information on the burst production mechanism. Previous studies on time correlations of repeating FRBs have been conducted using the distribution of wait times between two successive bursts (Wang & Yu 2017;Oppermann et al. 2018;Wang et al. 2018;Zhang et al. 2018;Li et al. 2019;Gourdji et al. 2019;Wadiasingh & Timokhin 2019;Oostrum et al. 2020;Tabor & Loeb 2020;Aggarwal et al. 2021;Cruces et al. 2021;Li et al. 2021;Zhang et al. 2021;Hewitt et al. 2022;Xu et al. 2022;Zhang et al. 2022;Du et al. 2023;Jahns et al. 2023;Sang & Lin 2023;Wang et al. 2023;Zhang et al. 2023a,b). The wait time distribution is known to be bimodal with a boundary of about 1 second, and the distribution on the long time side can be described by an uncorrelated Poisson process (Cruces et al. 2021; ★ E-mail: totani@astron.s.u-tokyo.ac.jp Hewitt et al. 2022;Jahns et al. 2023). The short-side distribution is thought to reflect a time scale related to the physical activity of the source or radiative processes, but its origin is unknown.
The wait time distribution does not take full advantage of timerelated statistics, because time correlations may exist not only between consecutive bursts but also across other bursts. Therefore in this study, we directly calculate the two-point correlation function in two-dimensional space of occurrence time and emission energy of FRBs, using the method widely used in cosmology to study largescale structures. We collected a wide range of reported observations including a large number of repeating FRB events from the three most active sources of FRB 20121102A, 20201124A, and 20220912A Hewitt et al. 2022;Xu et al. 2022;Zhang et al. 2022;Jahns et al. 2023;Zhang et al. 2023a), and nearly 7,000 bursts in total are analyzed to compute the correlation function.
At least some FRBs are known to occur in magnetars (neutron stars with extremely strong magnetic fields of 10 14 G) (Bochenek et al. 2020; The CHIME/FRB Collaboration 2020), and explosive phenomena in magnetars are believed to be triggered by starquakes in neutron star crusts, which is induced by magnetic energy (Kaspi & Beloborodov 2017). For this reason, various phenomena of neutron stars have been compared to earthquakes and solar flares, and similarities are often discussed (Cheng et al. 1996;GöǧüŞ et al. 1999). Therefore we will use the same method to analyze the time-energy correlation function of earthquake and solar flare data, and examine the similarities with the statistical properties of repeating FRBs.

The FRB data sets
We compute correlation functions for the seven data sets (Table 1) of FRBs observed by the two radio telescopes (Arecibo and FAST) from the following three FRB repeaters.
FRB 20121102A is the first discovered (Spitler et al. 2014(Spitler et al. , 2016, highly active, and most well-studied repeater located in a starforming dwarf galaxy at redshift = 0.193 (Bassa et al. 2017;Chatterjee et al. 2017;Tendulkar et al. 2017). The solar system barycentric time and emitted energy of bursts were taken from tables given in the original papers of the L21 , H22 (Hewitt et al. 2022), and J23 (Jahns et al. 2023) data sets. The observation logs (start and end times of each observation) are given in the papers of H22 and J23, and we obtained that for the L21 data directly from the authors. In the J23 data, all 1,027 sub-bursts in the 849 independent events are listed in their table, and we used all sub-bursts in the main analysis.
FRB 20201124A (The CHIME/FRB Collaboration 2021) is another repeater known for its high activity. It is in a Milky Way-like, barred spiral galaxy at = 0.0979 (Xu et al. 2022). Barycentric times of the X22 (Xu et al. 2022) and Z22 (Zhang et al. 2022) data are from the tables in the original papers. Burst energies of X22 were calculated from fluence and signal bandwidth given in their table, using and the corresponding luminosity distance, while for Z22 we used "the energies calculated with the central frequency" given in their table. We obtained the observation logs of the X22 data set directly from the authors. This repeater was extremely active in the Z22 data, with the highest event rate detected from a single FRB source. Particularly large numbers of bursts were detected on days three and four, and therefore, we analyzed the third and fourth days as independent data sets (Z22 D3 and Z22 D4, respectively) to examine the dependence on activity level. FRB 20220912A (McKinven & Chime/FRB Collaboration 2022) is located at a position that is consistent with a likely host galaxy of stellar mass ∼ 10 10 ⊙ at = 0.0771 (Ravi et al. 2022). We took barycentric times and energies of the 1,076 bursts as the Z23 (Zhang et al. 2023a) data from the table in the original paper.
The redshifts of these FRB sources are not large enough for cosmological effects to be significant, and hence no correction for the cosmological time dilation is made.

Correlation function calculations
For the difference of time (Δ ≡ 2 − 1 , 2 > 1 ) and energy (Δ lg ≡ lg 2 − lg 1 , where lg = log 10 ) of a burst pair (burst 1 and 2), the two-point correlation function (Δ , Δ lg ) is defined as the excess of the number of pairs over the uncorrelated case: where¯ is the expected pair number density in the uncorrelated case.
The method of estimating the correlation function is widely used in cosmology to study the spatial correlation of galaxies, and has been described in many references [e.g. §16.7 of Peacock (1999)]. The two-point correlation function is estimated by counting the number of pairs consisting of two objects and comparing them to a random sample made up of uncorrelated objects. To reduce statistical error, random data are usually generated with a much larger number of objects than the observed data, and in this study, the random sample is 100 times larger (10 4 times larger in terms of the number of pairs). The number of all possible pairs involved in a given bin of Δ -Δ lg space is denoted as for the real data sample and for the random sample, where and are appropriately normalized for different sample sizes. Then the correlation function can be estimated as but here we use which is an estimator that has less variance and is most widely used (Landy & Szalay 1993), where is the number of cross-pairs between the data and random samples. (We confirmed that using the former "natural" estimator has little effect on the conclusions of this paper.) The time correlation function (Δ ) is estimated in the same way, but without binning to the energy direction.
The random data to calculate¯ were generated as follows. The observation period for a single FRB data set spans multiple days, but a continuous observation is only a few hours per day. Therefore, we assume that it is a Poisson process with a constant rate and a constant energy distribution within a day. There are occasional interruptions during one day of observation, and these were taken into account in the random data generation according to the observation logs. The energy of the random data was generated assuming a cumulative distribution (< ), which was empirically constructed from the observed energies as (< ) = /(1 + ) ( = 0, . . . , + 1), where (< ) is linearly interpolated between < < +1 . Here, are the energies (in increasing order) of the observed event in one day for = 1, . . . , , and we set 0 = 1 − and +1 = + , where = ( − 1 )/( − 1) is the mean energy separation. Pairs across different days were not considered. Then pair counts within a day were added up for all observed days to compute the correlation function for a single data set.
The simplest way to estimate the correlation function errors is by Poisson error for the pair counts. However, if the correlation is nonzero, the Poisson error is a lower bound to the real error, and furthermore, errors between different bins are correlated. We, therefore, computed the covariance matrix of (Δ ) using the jackknife method (Norberg et al. 2009) by dividing the observation time of each day into 20 segments, where and denote the bins of Δ . The error 1/2 for the -th Δ bin determined by this method is not significantly different from that evaluated by Poisson statistics. Jackknife errors are sometimes smaller than Poisson errors, which is probably due to a sampling bias. Therefore, the larger of Poisson and jackknife errors are conservatively shown at each data point when errors of (Δ ) are presented. The catalog of time and energy obtained from observations is subject to bias due to selection effects. The most obvious effect is that bursts of weak radio fluxes are below the detection limit and do not enter the sample. However, since we are empirically constructing the energy distribution function for the random samples from the observed data, the incompleteness with respect to energy is correctly taken into account here. Another important effect is that when two independent events occur in close proximity in time, they are indistinguishable from sub-bursts within a single event, or the detection sensitivity to fainter bursts is reduced. This effect is, in The best-fit value, 68% confidence level (CL) lower and upper limits of the parameters in the fitting of (Δ ) = (Δ + ) − / − , where in [s] The branching ratio = ∫ (Δ ) (Δ ) When only independent events are counted by grouping sub-bursts together fact, important in the FRB samples analyzed here, and it will be discussed in detail in the next section. Figure 1 shows the event distribution in -space, the wait time distribution, and the pair distribution and correlation function in Δ -Δ lg space for the FRB 20121102A L21 data set (the other six are shown in Figs. A1-A6 in Appendix A).

Results
The time correlation functions (Δ ) for the seven data sets are shown together in Fig. 2 and separately in Fig. 3. No clear signal in the Δ > 1 s region indicates that the uncorrelated Poisson process can describe the data, which confirms previous studies. On the other hand, clear correlation signals are detected for all the data sets at Δ 1 s, corresponding to the shorter side of the bimodal wait time distribution. These signals exhibit power-law behavior in Δ ∼ 0.01-1 s, but are flatter in the shortest Δ region. In the region of Δ 0.01 s, the signal behavior differs among the seven data sets, likely due to different treatments about sub-bursts in a single event.
In general, if the radio flux rises again before falling to noise levels, it is considered a sub-structure within a single event. However, the exact criteria and whether to include sub-bursts in the catalog are not standardized (Hewitt et al. 2022;Jahns et al. 2023). This effect is clearly seen in the FRB 20121102A (J23) data set, where every subburst is explicitly listed in the catalog. When only the first sub-burst of an independent event (849 in total) is used in the analysis, becomes flat below 30 msec, whereas continues to increase smoothly with decreasing Δ until 1 ms when all 1,027 sub-bursts are used (Fig. 3). This suggests that even those classified as sub-bursts might be better treated as independent events.
We fit these signals with the function by minimizing 2 . Ideally, 2 should be calculated using the inverse of the covariance matrix (Norberg et al. 2009), −1 , but it is known that when the sample size is limited, inverting estimated by jackknife can be numerically unstable (Okumura et al. 2008;Pope & Szapudi 2008). In fact, we found that 2 using −1 varies sensitively with changes in the analysis parameters (e.g. binning), and reliable results could not be obtained. Since a precise determination of the model parameter values is not the main purpose of this study, error correlations between different Δ bins were not considered, and the larger of the Poisson and jackknife errors were employed for the error of the -th Δ bin. Perhaps the best way to properly incorporate error correlation is to perform Monte Carlo simulations with realistic theoretical modeling, which is beyond the scope of this study. The upper and lower 1 (68% CL) limits for the three model parameters ( , , and ) were determined by the excess Δ 2 from the minimum, which is expected to follow a 2 distribution with three degrees of freedom.
As seen in the case of FRB 20121102A (J23) (Fig. 3), the artificial effect of sub-burst removal results in a sharp drop of (Δ ) in the smallest Δ region. The drops seen in other data sets are also likely due to this effect, and fitting the functional form of eq. (4) to such a data set may produce biased results. Therefore, we excluded data points at Δ smaller than the peak of (Δ ) from the fit. In the large Δ region ( 1 s), no significant correlation signals are found. However, because of the large number of available pairs in this region, the statistical errors are very small, and consequently the fit becomes sensitive to systematic uncertainties. In particular, our analysis assumes that the burst rate is constant within a day, but if the actual rate is slowly changing on a time scale of about a day, it may appear as a nonzero signal in the correlation function. We are interested in the signal clearly detected at Δ < 1 s, and since it decreases rapidly as Δ increases, the signal is hidden by noise in large Δ regions. Therefore, data at Δ larger than the point at which first becomes negative were excluded from the fit. The best-fit parameter values and their errors of , , and are shown in Table 1, and the best-fit curves are shown in Fig. 3. The power law fits are particularly good for the two data sets (FRB 20121102A L21 and J23) extending to small Δ comparable with typical FRB durations (1-10 msec), with the index ∼ 1.6 or 2.3. The fact that is restricted to a relatively narrow range indicates that these data cannot be fitted by other function forms, such as the exponential function. For the FRB 20201124A (Z22 D3) data, pure power-law without the effect of non-zero is also consistent with the data within 68% CL, and no lower bound on (and corresponding upper bound on ) can be placed.
The data sets other than FRB 20121102A (L21) and (J23) are likely biased toward larger because small Δ pairs have been removed as sub-bursts, thus making the power-law features less clear. However, the power-law with ∼ 2 is within the 1 error range for all FRB data sets except for FRB 20201124A (X22). The value of the X22 data may be significantly larger than those in other sets, but it is due to the only data point with a large error at Δ = 230 msec (Fig. 3), and hence the possibility of statistical fluctuation cannot be ruled out. The lower limits of 2 and 3 are = 3.0 and 2.1, respectively, for the X22 data. It should also be noted that tends to be larger when is large due to degeneracy in the fit (see the next paragraph), and both and may be biased to a greater value because of sub-burst removal.
As for the upper limits on , the data sets other than FRB 20121102A (L21) and (J23) are consistent with extremely large values of > 30 within 1 . This means that (Δ ) can also be fitted with an exponential function, because the fitting function converges into exp(−Δ / ′ ) in the limit of → ∞ and → ∞ with ′ = / , which can be shown as where ≡ /Δ . Our parameter search is limited to ≤ 30, and if = 30 falls within the 1 error region, " > 30" and the corresponding value of are indicated in Table 1 as 1 upper bounds, but this is essentially an exponential fit with exp[−Δ /( / )].
In summary, at least two of the data sets show clear power laws of ∼ 2, and the other data sets are consistent with similar power laws. The power law feature indicates that there is no characteristic time scale for the correlation. It is then natural to assume that the shorter peak of the bimodal wait-time distribution does not reflect the activity duration, but rather correlated aftershocks like earthquakes. In fact, it is known that the aftershock occurrence rate of earthquakes follows a power-law ∝ (Δ + ) − with ∼ 1 (the Omori-Utsu law) (Omori 1894;Utsu 1957Utsu , 1961Utsu et al. 1995;de Arcangelis et al. 2016). This consideration provides a good motivation to investigate the similarity of correlation functions between FRBs and earthquakes.

The earthquake data sets
To quantify the similarity to earthquakes, we used the same analysis method to examine time-energy correlations of earthquakes. The earthquake data used in this study were extracted from the Japan Unified hI-resolution relocated Catalog for Earthquakes (JUICE, Yano et al. 2017). Five data sets of similar sample sizes to FRBs (∼ 1,000 events for each data set) were extracted from various regions, area sizes, and time periods ( Table 2). The three regions of Choshi, Narita, and Izumo were selected to see the regional dependences, and data from the greater Kanto region (including Choshi and Narita) were also analyzed to see the dependence on area size (Fig. 4). To keep the sample size similar to that of FRBs, the time period was adjusted so that the number of events in each data set is approximately 1,000. The rate of earthquakes changed dramatically before and after the huge Tohoku earthquake on March 11, 2011. We analyzed data from Narita before and after this to see the impact of the change in activity. The other three data sets are all before the earthquake. The energy of earthquakes was calculated from the magnitude by the formula = 10 11.8+1.5 erg. Magnitudes in the catalog are discretized with an accuracy of 0.1, and they are shifted randomly with an amplitude of less than ±0.05 to avoid creating pairs with exactly zero energy differences. Unlike FRBs, observations of earthquakes and solar flares were made continuously 24 hours a day. However, to make the analysis as similar as possible to FRBs, we divided a data set into 10 equal partial periods, with each corresponding to a single date of FRB observations. As with FRBs, pairs were looked for only within partial periods, and pair counts were added up over the entire period to compute . Results for the Narita b311 data are shown in Fig. 5, and those of the other four sets are in Figs. A7-A10 in Appendix A.

FRBs versus earthquakes
The time correlation functions of the earthquake data sets are shown in Figs. 6 and 7. They are similar to FRBs in that (Δ ) is a power law of ∝ (Δ + ) − at small Δ but uncorrelated events at a constant rate become dominant in the large Δ region. Denoting as the mean event rate including uncorrelated events (shown in Table 1 and  Table 2), (1 + ) becomes the aftershock rate, i.e., the rate of events occurring Δ after an event, by definition of the correlation function. This is shown in the bottom panel of Figs. 2 (FRBs) and 6 (earthquakes). Since FRB event rates vary on different observation dates within a single data set, we calculated for the entire data set by taking an average weighted by the number of pairs on each date (or partial period in earthquake data).
The expected number of correlated aftershocks to an event (called "the branching ratio" in earthquake studies) is given by ≡ ∫ (Δ ) (Δ ), and if is less than 1, the effect of a single event causing multiple aftershocks is small. The condition for < 1 is roughly < (Δ ) −1 , which is satisfied in both the FRB and earthquake data, as seen in Figs. 2 and 6. The branching ratio was calculated using the best-fit parameter values of (Δ ), and shown in Tables 1 and 2. Since is smaller than 1 in some of the earthquake data sets, the upper bound Δ of the integration for was set to the point where (Δ ) = 1 or 0.1 to avoid divergence. The values of are indeed less than 1, but of order unity ( 0.1) for all FRB and earthquake data sets.
Then these results can be interpreted as each event produces at most one aftershock at a rate consistent with the Omori-Utsu law, ∝ (Δ + ) − , where all events are treated on the same footing and there is no distinction between mainshocks and aftershocks. Such a model is called ETAS (epidemic-type aftershock sequence) in earthquake studies and is known to explain earthquake data well (Ogata 1999;Saichev & Sornette 2006;de Arcangelis et al. 2016).
Further similarities between FRBs and earthquakes can be pointed out. The branching ratio is about 0.1-0.5 for both FRBs and earthquakes (Tables 1 and 2), which is also similar to those found in past earthquake studies (Saichev & Sornette 2006;de Arcangelis et al. 2016). The power law flattens at Δ , and is about 1-10 msec for !"#$%&'%''&%()*%'+ !"#$%&%%&,'%($)-%.+ !"#$%&'%''&%()/%%+ !"#$%&'%''&%()0%.+  FRBs and 0.3-3 min for earthquakes, which are close to the typical duration of the respective phenomena. It should also be noted that the correlated aftershock rate ( ) is stable both in FRBs and earthquakes (bottom panels of Figs. 2 and 6), even though the mean event rate of uncorrelated events varies widely with changes of source activity (and also with the choice of region and its area size for earthquakes). The aftershock rates of the three different FRB sources are not significantly different, but the rate fluctuations are even smaller within each source. This indicates that aftershocks are not caused by the activity of the source as a whole, but by the changes induced by each individual event. Finally, there is little correlation between Δ and Δ lg in either FRBs or earthquakes ( almost constant against Δ lg in Figs. 1, 5, and A1-A10).
In fact, a slight time-magnitude correlation has been reported by detailed statistical tests for earthquakes (Lippiello et al. 2008;de Arcangelis et al. 2016;Zhang et al. 2023b). Asymmetry with respect to Δ lg (more negative Δ lg pairs, i.e., aftershock energy smaller than mainshock) may also be evident in some of our earthquake data sets: Narita (b311, Fig. 5), Choshi (Fig. A8), and Izumo (Fig. A10). Similar asymmetry about Δ lg is not visible in the FRB !"#$%&'()*+()$, data, but it may be due to a narrow dynamic range of caused by detection limits. Careful verification by future studies would be desirable.
The only difference between the two phenomena is the value of the Omori-Utsu index, : FRBs have a larger resulting in the bimodal wait-time distribution, while bimodality is not clearly visible in the earthquake data. It should be noted, however, that of earthquake data varies widely depending on regions, and those of FRBs are marginally within the range for earthquakes ( = 0.6-1.9) reported by past studies (Utsu et al. 1995;Ogata 1999;Wiemer & Katsumata 1999;de Arcangelis et al. 2016). Wang et al. (2018) found that the rate evolution of 14 bursts detected during an observation of about 60 minutes from FRB 20121102A is consistent with the Omori-Utsu law. This result is likely unrelated to the correlations found in this study because the correlation time scale of Wang et al. (2018) is much larger than Δ ∼ 1 s, where we did not find significant correlations. It is difficult to draw strong conclusions from the sample of Wang et al. (2018), not only because of the small statistics but also because of the arbitrariness of where to separate mainshocks and aftershocks. It is possible that the rate variation observed by Wang et al. (2018) does not reflect a systematic aftershock sequence, but rather a change in the source activity during that observation period. Enoto et al. (2017) analyzed the light curves of X-ray outbursts in three magnetars and found that they can be fitted by the Omori-Utsu law, and proposed the possibility that X-ray luminosity is superpositions of many small starquakes. However, the plateau time scale is longer than 10 days, which is many orders of magnitude larger than the time correlations found in this work, suggesting that there is no direct relationship between the two. It should be noted that continuous luminosity evolution can also be explained by continuous evolution of energy release in a single system, rather than by the superposition of many discrete aftershocks. In general, in systems where the energy loss rate is proportional to the power of the system energy ( ∝ ), the evolution of takes the form of the Omori-Utsu formula with = /( − 1) (e.g., = 2 for a rotation-powered pulsar with a constant magnetic field strength).

Comparison with solar flares
The solar flare data were taken from the Hinode catalog (Watanabe et al. 2012). To see the dependence on activity level, data were extracted from two periods (200 days from 2012 Apr. to Oct., and 1,200 days from 2017 Oct. to 2021 Jan.) of high and low solar activity, adjusting the time period so that each period includes about 1,000 events (1,422 and 1,207 events, respectively). The Xray flux was calculated from the flare class in the catalog, and this was used as a proxy for energy . Flare classes are given as two significant digits and, as with the earthquake data, they were randomly shifted by the precision of the significant digits [e.g., X-ray flux of a B4.2 class flare is randomly chosen from the range of (4.15-4.25) × 10 −4 erg cm −2 s −1 ]. These data sets were analyzed in the same way as earthquakes, and the results for the low activity set is shown in Fig. 8 (and the high activity set in Fig. A11 in Appendix A). Time correlation functions for the two sets are shown in Fig. 9. Figure 9 shows that a significant time correlation signal is detected in the region of Δ 100 hr. It is known that flares occurring in different active regions on the solar surface are uncorrelated (Benz Figure 5. The same as Fig. 1, but for the Narita (before 311) earthquake data. 2017), and hence the correlated signals are thought to be between flares occurring in a common active region. We can convert into the correlated aftershock rate , and it is stable regardless of the level of solar activity ( ), which is similar to the behaviors of FRBs and earthquakes. However, reaches at most 5 without showing a clear power law, in sharp contrast to FRBs and earthquakes. The correlated rate is well above (Δ ) −1 , which indicates that the branching ratio is larger than 1 and multiple aftershocks are generated from a single event. Another difference from FRBs or earthquakes is the strong correlation not only in time but also in energy, with peaking at Δ ∼ 10 3 -10 4 s and Δ lg ∼ 0 (lower-right panels of Figs. 8 and A11). This strong time-energy clustering nature is evident in thedistribution shown in these figures (upper-left panels), which are not seen in FRB or earthquake data.
Though similarities between solar flares and magnetars are often discussed, these results indicate a different nature of the time-energy correlation of solar flares compared to FRBs and earthquakes. Interpretation of these results in comparison with the physics and theory of solar flares is outside the scope of this study. Time and/or energy correlations of solar flares have been discussed by analyzing the wait time distributions and cross-correlations between energy and time (either since the last event or to the next event) (Wheatland 2000;Lippiello et al. 2010;Hudson 2020;Carlin et al. 2023). It is not straightforward to compare our results with these, because Δ in our analysis is not wait time but time separation for arbitrary pairs, and we are looking at the energy difference Δ lg of an event pair, rather than the energy itself of one event. Careful comparison of these results to elicit information about the physics of solar flares is an interesting subject for future research. In such studies, it is important to take into account observational biases in the solar flare catalog, such as the obscuration effect (lower detection efficiency for faint flares after a giant flare, Wheatland 2001).

CONCLUSIONS AND DISCUSSION
By examining the correlation functions in time-energy space, we found remarkable similarities between the statistical properties of FRBs and earthquakes, especially the laws on aftershock occurrence. Listing the similarities, (1) the aftershock rate follows a power law of ∝ (Δ + ) − (the Omori-Utsu law), (2) is about the same as the typical event duration, (3) the branching ratio (expected number of aftershocks associated with a single event) is = 0.1-0.6 for both phenomena, (4) the correlated aftershock rates remain stable regardless of the change in the source activity or mean event rate, and (5) there is little correlation between energy and time.
In contrast, the correlation functions for solar flares differ significantly from those of FRBs and earthquakes. The correlation function cannot be fitted by a power law, and the branching ratio is significantly greater than unity, which means that a single event causes multiple aftershocks. There is also a strong correlation in the direction of energy, with flares of similar energy tending to occur in succession. It is well known that the energy distribution of solar flares follows a power law, which is often noted to be similar to that of earthquakes (the Wadati-Gutenberg-Richter law 1 , Wadati 1932;Gutenberg & Richter 1944). However, with respect to correlations in the time-energy space, solar flares do not resemble earthquakes, and by comparison, the similarity between FRBs and earthquakes is remarkable.
Both solar flares and magnetars are caused by magnetic energy, but unlike the Sun, neutron stars are thought to have solid crusts on their surfaces. Therefore, the most natural interpretation of the present results is that FRBs are earthquake-like phenomena that suddenly release the energy stored in neutron star crusts. Other models of FRB repeaters are not immediately dismissed, but any model must be consistent with this time-energy correlation, placing strong constraints on possible models.
In the case of earthquakes, differences in the index are thought to reflect the physical properties of crust and seismic processes (Ogata 1999;Wiemer & Katsumata 1999;de Arcangelis et al. 2016). Then relatively large values of FRBs may provide information about the physical properties of neutron star crusts, and energy production mechanism by seismic processes in them. This suggests a new possibility for future studies to probe the physics of dense nuclear matter by using repeater FRBs.
Two FRB repeaters, including 20121102A, are known to periodically change their activity in cycles of 16 or 160 days (Amiri et al. 2020;Rajwade et al. 2020;Cruces et al. 2021). If these are neutron stars in a binary system, the activity may be increased when the crust is deformed by tidal forces from the companion star at the pericenter. Models of binary origin have been proposed for the periodic FRB activities (Ioka & Zhang 2020;Lyutikov et al. 2020;Barkov & Popov 2022), but most have assumed that the periodicity is due to the absorption of FRB radiation (e.g. by stellar winds from the companion star), rather than intrinsic change of the FRB production activity. A repeater FRB has been found in a globular cluster of the nearby galaxy M81 (Kirsten et al. 2022), and the FRB activity may be stimulated by close encounters with other stars in the cluster.
As far as we know, most previous studies on the time correlation of earthquakes or solar flares have also been based on wait times rather than correlation function (de Arcangelis et al. 2016;Saichev & Sornette 2006;Wheatland 2000). The correlation function method adopted here fully exploits the temporal information of the sample, and the time dependence of the correlated aftershock rate can be seen more directly. This is why the two phenomena (FRBs and earthquakes), that appear to be different in the wait time distribution, turn out to have essentially the same properties when viewed in terms of the correlation function. Applying this method to a wider range of earthquake and solar flare data may provide new insights into these phenomena.
This paper has been typeset from a T E X/L A T E X file prepared by the author.