# Enhancing Trace Visualizations for
Microservices Performance Analysis

Jessica Leone<EMAIL_ADDRESS>0000-0002-2870-8161 University
of L’Aquila, Italy and  Luca Traini<EMAIL_ADDRESS>0000-0003-3676-0645
University of L’Aquila, Italy

###### Abstract.

Performance analysis of microservices can be a challenging task, as a typical
request to these systems involves multiple Remote Procedure Calls (RPC)
spanning across independent services and machines. Practitioners primarily
rely on distributed tracing tools to closely monitor microservices
performance. These tools enable practitioners to trace, collect, and visualize
RPC workflows and associated events in the context of individual end-to-end
requests. While effective for analyzing individual end-to-end requests,
current distributed tracing visualizations often fall short in providing a
comprehensive understanding of the system’s overall performance. To address
this limitation, we propose a novel visualization approach that enables
aggregate performance analysis of multiple end-to-end requests. Our approach
builds on a previously developed technique for comparing structural
differences of request pairs and extends it for aggregate performance analysis
of sets of requests. This paper presents our proposal and discusses our
preliminary ongoing progress in developing this innovative approach.

Microservices, Distributed Tracing, Performance Analysis

††copyright: none††ccs: Software and its engineering Software
performance††ccs: Software and its engineering Maintaining software††ccs:
Software and its engineering Software evolution††ccs: Human-centered computing
Visual analytics††ccs: Human-centered computing Visualization toolkits

##  1\. Introduction

Microservices have emerged as a paradigm shift in software development,
enabling the efficient testing and deployment of new software features and
improvements (Newman, 2015). Due to their modular nature, microservices are
particularly well-suited for today’s software market, where the ability to
rapidly release software updates and improvements is perceived as a key
competitive advantage (Rubin and Rinard, 2016).

However, despite these advantages, adopting microservices also poses
significant challenges. A main challenge, for example, is ensuring adequate
software performance. Indeed, a proactive approach to performance assurance
(_e.g.,_ pre-production performance testing (Jiang and Hassan, 2015; Laaber
and Leitner, 2018; Traini et al., 2022)) is often impractical in these
contexts, given the inherent complexity of these systems (Sridharan, 2017;
Veeraraghavan et al., 2016), and the significant pressure to deliver fast-to-
market (Rubin and Rinard, 2016; Traini, 2022). Performance behavior of
microservices systems usually _emerges_ in the field (Veeraraghavan et al.,
2016) as a complex interaction of several RPCs that spans multiple independent
services and machines. Furthermore, these systems are constantly undergoing
software changes (_e.g.,_ multiple releases per day) and are exposed to highly
varying workloads (Ardelean et al., 2018), which make them susceptible to
unexpected performance regressions (Veeraraghavan et al., 2016; Traini et al.,
2021).

To aid performance analysis, practitioners typically rely on distributed
tracing (Parker et al., 2020) (_e.g.,_ Jaeger (Technologies, 2023), Dapper
(Sigelman et al., 2010), Canopy(Kaldor et al., 2017)) as a means to monitor
and inspect microservices performance in production. Distributed tracing tools
capture the workflow of causally-related events (_i.e.,_ work done to process
a request) within and among the components of a microservices system
(Sambasivan et al., 2016), and provide visualizations to support performance
analysis of individual end-to-end requests, _e.g.,_ swimlane visualizations
(Davidson et al., 2023; Sigelman et al., 2010; Technologies, 2023).

Although these visualizations are particularly effective for analyzing
individual requests, they are limited when it comes to provide a comprehensive
picture of microservices performance. Understanding request performance in
_aggregate_ (Parker et al., 2020) often requires switching between various
visualization tools, making the process cumbersome and time-consuming
(Davidson et al., 2023). Current distributed tracing tools do not provide out-
of-the-box support for aggregate performance analysis.

In this paper, we propose a novel visualization approach that facilitates
comprehension of the relationship between request characteristics and end-to-
end response time behavior. The proposed approach enables seamless performance
analysis of multiple requests, eliminating the need to switch between multiple
tools and visualizations. Our approach builds upon a previously developed
visualization method for comparing request pairs (Farro, 2018) and extends it
to accommodate various types of aggregate analysis on sets of requests. This
new approach will empower engineers to swiftly identify significant and
recurring request characteristics that are linked to specific end-to-end
response time behavior, thereby enabling a more effective performance analysis
of microservices. Here, we present our proposal and report on our preliminary
ongoing progress in developing this novel visualization approach.

##  2\. Motivation

Figure 1. Swimlane visualization.

Distributed tracing tools provide swimlane visualizations (Kaldor et al.,
2017; Sigelman et al., 2010; Technologies, 2023; Davidson and Mace, 2022) as a
means to investigate end-to-end response time (refer to Fig. 1 for an
illustration). These visualizations are particularly useful for performance
analysis of individual requests. For example, they can be used to investigate
how each RPC contributes to the end-to-end response time of a specific
request, and pinpoint potential bottlenecks.

Despite their usefulness, swimlane visualizations have certain limitations
when it comes to more complex performance analyses. For instance, analyzing
the performance of individual end-to-end requests can be misleading without
considering the appropriate context (Davidson et al., 2023). The response time
of a request can only be considered anomalous when compared to other requests
of the same type (Anand et al., 2020). Additionally, engineers are often more
interested in investigating overall response time trends rather than the
performance of individual requests (Parker et al., 2020). Different end-to-end
response time behaviors can be associated with specific request
characteristics, such as particular RPC execution paths and/or RPC performance
behaviors, and engineers can be interested in identifying these
characteristics to spot potential performance bugs and/or optimization
opportunities (Parker et al., 2020; Krushevskaja and Sandler, 2013;
Cortellessa and Traini, 2020).

Distributed tracing tools lack support for this type of analysis, which
currently requires the use of multiple visualizations and tools in tandem
(Davidson et al., 2023), _e.g.,_ Jaeger (Technologies, 2023) and Kibana
(Elastic, 2023b). A common approach is to first recognize recurrent
performance behaviors that warrant investigation, and then examine individual
requests to characterize relevant performance behaviors. This is usually
achieved by detecting “modes” of the end-to-end response time distribution
(_e.g.,_ using Kibana), which indicate meaningful recurrent performance
behaviors. Then, samples of requests associated with each mode can be
extracted and individually examined (_e.g.,_ using Jaeger’s swimlanes) to
identify unique characteristics that lead to specific performance behaviors,
_i.e.,_ modes.

Unfortunately, this approach can be particularly tedious as it necessitates
the manual examination and comparison of multiple requests across multiple
visualizations and tools. Furthermore, even when successful, it falls short of
providing a satisfactory level of confidence. Identifying specific
characteristics of a particular mode requires verifying that these
characteristics appear _exclusively_ in requests belonging to that mode, which
can be challenging using current tools.

Figure 2. End-to-end response time distribution.

For instance, consider the scenario illustrated in Fig. 2, which depicts the
distribution of the end-to-end response time for a specific type of request,
such as loading the homepage of a website. As depicted in the figure, requests
exhibit three distinct response time behaviors, _i.e.,_ modes. Suppose that
the rightmost mode is characterized by a unique request characteristic, such
as a specific RPC execution path. That is, a specific RPC is invoked in all
requests belonging to the rightmost mode, but not in the others (_e.g.,_ due
to a cache miss). Identifying patterns like this can be particularly
challenging with current distributed tracing tools.

##  3\. Proposal

Figure 3. Dashboard.

We present a novel interactive visualization approach for aggregate
performance analysis of end-to-end requests. Our approach provides a dashboard
comprising two main interactive components: a tree and a histogram (see Fig.
3). The interactive tree offers an overview of the requests’ workflows in
terms of RPC execution paths, while the histogram illustrates the
corresponding end-to-end response time distribution. The key insight behind
this approach is to enable the simultaneous analysis of various aspects of
requests and the relationship among these aspects through a unified dashboard.

The tree visualization proposed in our approach builds upon a previous
visualization technique provided in Jaeger (Technologies, 2023), which enables
the comparison of two requests in terms of RPC execution paths (Farro, 2018).
We extend this technique by proposing a tree visualization that enables the
analysis of an entire set of requests. Each node of the tree represents a RPC
invocation within a specific execution path, where the leftmost node
represents the root RPC, and edges indicate direct RPC invocation. For
example, in Fig. 3, the red node labeled as getcheapestroute indicates the RPC
execution path getbycheapest$\rightarrow$
getleftticket$\rightarrow$getcheapestroute. It is worth noting that when a
particular RPC invokes the same RPC multiple times, this leads to a single
node in the tree. In other words, if the RPC getbycheapest invokes the RPC
getleftticket multiple times, there will be only one child node referring to
getleftticket. The rationale behind this approach is to address the potential
complexity of requests by providing a more concise visualization. An RPC
execution path will appear in the tree only if it is present in at least one
request within the set of requests being analyzed.

Similar to the trace comparison tool of Jaeger (Farro, 2018), we utilize color
encoding to highlight variations within a set of requests. However, in
contrast to this approach, our color encoding scheme denotes the degree of
variance of a particular request characteristic within a group of requests,
rather than the difference between two requests. Furthermore, we employ a
continuous color scale, as opposed to the discrete one used by Jaeger (Farro,
2018). As an example, we can use color encoding to depict the variability in
the number of invocations of a specific RPC execution path within a group of
requests. A white node indicates that the corresponding RPC execution path
appears the same number of times in all the requests under analysis.
Conversely, a red node indicates a substantial change in the number of
invocations of the corresponding RPC execution path from one request to
another. We use standardized measures of dispersion, such as the Coefficient
of Variation (CV) (Everitt, 1998), to encode these differences into colors.
For instance, a CV of 0 results in a white node, while a CV $\geq$ 1 results
in a red node. The shade of color gradually transitions from white to red as
the CV value increases, with 0 $\leq$ CV $\leq$ 1.

This approach facilitates the rapid identification of RPC execution paths that
exhibit divergent request characteristics. An engineer can delve deeper into a
specific RPC execution path by interacting with its corresponding node,
bringing up a new visualization such as the pie chart displayed in Fig. 3.
This allows for a more in-depth examination of the divergent characteristics
observed in that specific RPC execution path. The pie chart in Fig. 3
illustrates, for instance, two distinct behaviors, with 77% of requests not
involving any invocations of the path getbycheapest$\rightarrow$
getleftticket$\rightarrow$getcheapestroute, and 23% of requests involving one
invocation of such path. By interacting with the corresponding slice of the
pie chart, the tool illustrates how this particular request characteristic
affects end-to-end response time by highlighting the corresponding area of the
histogram. For instance, Fig. 3 shows that the rightmost mode of the end-to-
end response time distribution is characterized by the presence of the
particular RPC execution path getbycheapest$\rightarrow$
getleftticket$\rightarrow$getcheapestroute.

Another advantage of this approach is its versatility in analyzing various
request characteristics. This can be accomplished by changing the color
encoding assigned to tree nodes. For instance, the method can be adapted to
examine how different RPC execution times influence the distribution of end-
to-end response time. Or, similarly, it can be utilized to investigate the
impact of different HTTP header values on response time (in these situations,
other statistical metrics such as the Gini index (Gini, 1936) or Entropy
(Shannon, 1948) can be applied). Furthermore, we plan to make the analysis
bidirectional, allowing the engineer to begin their investigation from either
the tree or the histogram to understand response time behavior. If the
engineer chooses to start from the response time distribution, they can select
a specific range using a slider selector. This selection will result in an
updated color scheme on the tree, where the colors will now indicate the
divergence of the selected set of requests compared to all others (using
statistical measures such as Kullback-Leibler divergence (Kullback and
Leibler, 1951)). This approach will enable engineers to identify how certain
modes of response time distribution correlate with specific request
characteristics.

Figure 4. Workflow.

We envision two potential challenges for the development of our visualization
approach. The first challenge concerns the efficiency of our tool. Modern
distributed tracing tools collect huge volumes of traces per day (Kaldor et
al., 2017), which can make the interaction with the proposed visualizations
computationally intensive. In order to enable an efficient generation of
visualizations, we plan to preprocess traces in an optimized format. Figure 4
depicts the workflow of our approach. The traces of the microservices system
are continuously collected by the distributed tracing tool (_i.e.,_ _Trace
Collector_) and persisted in a _Trace Storage_ , such as ElasticSearch
(Elastic, 2023a). The stored traces are then periodically preprocessed and
persisted in dedicated storage in an optimized format. The visual dashboard
application can then directly query the optimized trace storage to enable the
efficient generation of visualizations.

The second challenge we envision is related to the user experience. When too
much information is displayed on the screen, the user can become overwhelmed
(Davidson et al., 2023). To mitigate this, we intend to keep our dashboard as
minimalistic as possible to avoid placing an excessive burden on the user. For
instance, certain request types may involve a large number of RPCs, which can
result in an extremely large tree structure. To reduce the user’s cognitive
load, we will allow the user to hide the RPCs invoked within a particular
execution path by double-clicking on the corresponding node. This will enable
the user to navigate the tree more easily.

##  4\. Ongoing work

We are currently engaged in the development of a prototype for our
visualization approach. For our pilot study, we are utilizing a dataset of
traces that was generated in a previous work of the second author (Traini and
Cortellessa, 2021). The dataset was created using the widely-used
microservices benchmark system _TrainTicket_ (Zhou et al., 2021; Li et al.,
2021; Zhang et al., 2022; Guo et al., 2020; Zhou et al., 2019), and the load
generator _PPTAM_ (Avritzer et al., 2019). Our technology stack includes the
use of Jaeger (Technologies, 2023) as the distributed tracing collector and
Elasticsearch (Elastic, 2023a) for storing both common and optimized traces.
The dashboard has been developed using D3.js (Bostock, 2023) for
visualizations and Flask (Pallets, 2023) for the backend service.

At present, we have accomplished the implementation of the UI of our working
prototype, along with its key visualization components, namely the interactive
tree and histogram. However, we are still engaged in the ongoing development
of the preprocessing engine and a portion of the server infrastructure. For
this reason, we are unable to report an authentic screenshot of our working
prototype, at this stage. Nonetheless, we do provide a mock-up of the expected
outcome in Fig. 3.

Our next step entails the completion of the first version of our prototype.
This first version will be centered exclusively on the analysis of the number
of invocations of RPC execution paths and their relationship with end-to-end
performance behavior. In subsequent stages, we plan to extend the capability
of the prototype beyond this specific request characteristic by exploring
other metrics such as RPC execution time. Additionally, we intend to enable
bidirectional analyses by supporting investigations initiated from both the
tree and the histogram (currently, our working prototype only supports a
singular direction, _i.e.,_ starting from the tree).

It is noteworthy that our current design only supports the individual
inspection of specific request characteristics and does not enable the
simultaneous analysis of multiple request characteristics in tandem. This
supplementary capability may become crucial for more elaborate performance
analyses. Thus, as a long-term goal, we intend to broaden our approach by
enabling the analysis of combined sets of multiple request characteristics and
their interrelation with end-to-end performance.

Regarding the evaluation of our proposed approach, we plan to continue using
_TrainTicket_ as the reference case study, given its relevance in the software
engineering community. To further enhance the assessment of our proposed
approach, we will consider a variety of workload scenarios, including some
that involve synthetically injected performance bugs (Traini and Cortellessa,
2021). Furthermore, we aim to experiment with real-world distributed traces
from large microservices systems to facilitate the transition of our approach
to practice. For instance, we plan to utilize the dataset of traces recently
shared by Alibaba (Luo et al., 2021).

##  5\. Related work

Previous research on visualization approaches for distributed systems has
traditionally focused on the analysis of individual requests or comparison
between two requests. ShiViz (Beschastnikh et al., 2020) visualizes a logged
distributed execution as a time-space diagram. The study of Sambavisan  _et
al._ (Sambasivan et al., 2013) compares three well-known visualization
approaches in the context of presenting the results of one automated
performance root cause analysis approach (Sambasivan et al., 2011). TraVista
(Anand et al., 2020) builds on the prevalent approach of visualizing a single
request by augmenting this view with information to aid the user in
contextualizing the performance of one request with respect to other requests.
Jaeger (Technologies, 2023) enables the comparison of structural aspects of
two requests (Farro, 2018). Commercial APM tools (Ahmed et al., 2016), such as
Dynatrace (Inc., 2023b), AppDynamics (Inc., 2023a), or Instana (IBM, 2023),
provide features that support aggregated analysis of end-to-end requests. For
instance, Dynatrace’s Service Flow feature (Inc., 2023c) shows aggregate
workflows of end-to-end requests along with their associated characteristics.
As far as we know, current APM tools do not provide dedicated visualizations
to analyze the relationship between requests’ characteristics and end-to-end
performance behavior. Other studies have proposed automated methods to aid the
detection of patterns in requests’ characteristics associated with anomalous
end-to-end performance behavior (Traini and Cortellessa, 2021; Bansal et al.,
2020; Cortellessa and Traini, 2020; Krushevskaja and Sandler, 2013).

Recent research by Davidson and Mace (Davidson and Mace, 2022) emphasizes the
importance of visualization in the context of systems research and encourages
further effort into this topic. Davidson _et al._ (Davidson et al., 2023)
performed a qualitative interview study to identify shortcomings of existing
distributed tracing tools, and expose various open research problems that
involve several domains (including visualization research).

##  6\. Conclusion

In this paper, we presented a novel visualization approach for microservices
performance analysis. Our approach addresses the limitations of current
distributed tracing visualizations by enabling aggregate performance analysis
of multiple end-to-end requests. Further studies are required to evaluate the
effectiveness of our proposal and gather feedback from practitioners.

###### Acknowledgements.

This work is supported by“Territori Aperti” (a project funded by Fondo
Territori Lavoro e Conoscenza CGIL, CSIL and UIL), and by European Union -
NextGenerationEU - National Recovery and Resilience Plan (Piano Nazionale di
Ripresa e Resilienza, PNRR) - Project: “SoBigData.it - Strengthening the
Italian RI for Social Mining and Big Data Analytics” - Prot. IR0000013 -
Avviso n. 3264 del 28/12/2021. We would like to thank Antinisca Di Marco and
Giovanni Stilo for the useful suggestions and comments that were helpful in
improving the paper.

## References

  * (1)
  * Ahmed et al. (2016) Tarek M. Ahmed, Cor-Paul Bezemer, Tse-Hsun Chen, Ahmed E. Hassan, and Weiyi Shang. 2016.  Studying the Effectiveness of Application Performance Management (APM) Tools for Detecting Performance Regressions for Web Applications: An Experience Report. In _Proceedings of the 13th International Conference on Mining Software Repositories_ (Austin, Texas) _(MSR ’16)_. Association for Computing Machinery, New York, NY, USA, 1–12.  https://doi.org/10.1145/2901739.2901774
  * Anand et al. (2020) Vaastav Anand, Matheus Stolet, Thomas Davidson, Ivan Beschastnikh, Tamara Munzner, and Jonathan Mace. 2020.  Aggregate-Driven Trace Visualizations for Performance Debugging.  arXiv: 2010.13681 [cs] Number: arXiv:2010.13681.
  * Ardelean et al. (2018) Dan Ardelean, Amer Diwan, and Chandra Erdman. 2018.  Performance Analysis of Cloud Applications. In _Proceedings of the 15th USENIX Conference on Networked Systems Design and Implementation_ _(NSDI’18)_. USENIX Association, USA, 405–417.  event-place: Renton, WA, USA.
  * Avritzer et al. (2019) Alberto Avritzer, Daniel Menasché, Vilc Rufino, Barbara Russo, Andrea Janes, Vincenzo Ferme, André van Hoorn, and Henning Schulz. 2019.  PPTAM: Production and Performance Testing Based Application Monitoring. In _Companion of the 2019 ACM/SPEC International Conference on Performance Engineering_ (Mumbai, India) _(ICPE ’19)_. Association for Computing Machinery, New York, NY, USA, 39–40.  https://doi.org/10.1145/3302541.3311961
  * Bansal et al. (2020) Chetan Bansal, Sundararajan Renganathan, Ashima Asudani, Olivier Midy, and Mathru Janakiraman. 2020\.  DeCaf: Diagnosing and Triaging Performance Issues in Large-Scale Cloud Services. In _Proceedings of the ACM/IEEE 42nd International Conference on Software Engineering: Software Engineering in Practice_ (Seoul, South Korea) _(ICSE-SEIP ’20)_. Association for Computing Machinery, New York, NY, USA, 201–210.  https://doi.org/10.1145/3377813.3381353
  * Beschastnikh et al. (2020) Ivan Beschastnikh, Perry Liu, Albert Xing, Patty Wang, Yuriy Brun, and Michael D. Ernst. 2020\.  Visualizing Distributed System Executions.  _ACM Trans. Softw. Eng. Methodol._ 29, 2 (March 2020).  https://doi.org/10.1145/3375633 Place: New York, NY, USA Publisher: Association for Computing Machinery.
  * Bostock (2023) Mike Bostock. 2023\.  D3: Data-Driven Documents.  https://d3js.org ”Accessed 2023-01-21 19:45”.
  * Cortellessa and Traini (2020) Vittorio Cortellessa and Luca Traini. 2020.  Detecting Latency Degradation Patterns in Service-Based Systems. In _Proceedings of the ACM/SPEC International Conference on Performance Engineering_ _(ICPE ’20)_. Association for Computing Machinery, New York, NY, USA, 161–172.  https://doi.org/10.1145/3358960.3379126
  * Davidson and Mace (2022) Thomas Davidson and Jonathan Mace. 2022.  See it to believe it? The role of visualisation in systems research. In _Proceedings of the 13th Symposium on Cloud Computing_ _(SoCC ’22)_. Association for Computing Machinery, New York, NY, USA, 419–428.  https://doi.org/10.1145/3542929.3563488
  * Davidson et al. (2023) Thomas Davidson, Emily Wall, and Jonathan Mace. 2023\.  A Qualitative Interview Study of Distributed Tracing Visualisation: A Characterisation of Challenges and Opportunities.  _IEEE Transactions on Visualization and Computer Graphics_ (2023), 1–12.  https://doi.org/10.1109/TVCG.2023.3241596
  * Elastic (2023a) Elastic. 2023a.  Elasticsearch: The Official Distributed Search & Analytics Engine.  https://www.elastic.co/elasticsearch/ ”Accessed 2023-01-15 18:38”.
  * Elastic (2023b) Elastic. 2023b.  Kibana: Explore, Visualize, Discover Data.  https://www.elastic.co/kibana ”Accessed 2023-01-15 17:38”.
  * Everitt (1998) Brian Everitt. 1998\.  _The Cambridge Dictionary of Statistics_.  Cambridge University Press, Cambridge, UK.
  * Farro (2018) Joe Farro. 2018\.  Trace comparisons arrive in Jaeger 1.7.  https://medium.com/jaegertracing/trace-comparisons-arrive-in-jaeger-1-7-a97ad5e2d05d ”Accessed 2023-01-15 18:16”.
  * Gini (1936) Corrado Gini. 1936\.  _On the measure of concentration with special reference to income and statistics_.  Colorado College Publication, General Series.
  * Guo et al. (2020) Xiaofeng Guo, Xin Peng, Hanzhang Wang, Wanxue Li, Huai Jiang, Dan Ding, Tao Xie, and Liangfei Su. 2020\.  Graph-Based Trace Analysis for Microservice Architecture Understanding and Problem Diagnosis. In _Proceedings of the 28th ACM Joint Meeting on European Software Engineering Conference and Symposium on the Foundations of Software Engineering_ (Virtual Event, USA) _(ESEC/FSE 2020)_. Association for Computing Machinery, New York, NY, USA, 1387–1397.  https://doi.org/10.1145/3368089.3417066
  * IBM (2023) IBM. 2023.  Instana: Automated Observability.  https://www.instana.com ”Accessed 2023-01-16 14:55”.
  * Inc. (2023a) Cisco Systems Inc. 2023a.  AppDynamics: Observability Platform.  https://www.appdynamics.com ”Accessed 2023-01-16 14:53”.
  * Inc. (2023b) Dynatrace Inc. 2023b.  Dynatrace: Modern cloud done right.  https://www.dynatrace.com ”Accessed 2023-01-16 14:50”.
  * Inc. (2023c) Dynatrace Inc. 2023c.  Dynatrace. Service Flow.  https://www.dynatrace.com/platform/service-flow/ ”Accessed 2023-02-14 14:10:59”.
  * Jiang and Hassan (2015) Zhen Ming Jiang and Ahmed E. Hassan. 2015.  A Survey on Load Testing of Large-Scale Software Systems.  _IEEE Transactions on Software Engineering_ 41, 11 (2015), 1091–1118.  https://doi.org/10.1109/TSE.2015.2445340
  * Kaldor et al. (2017) Jonathan Kaldor, Jonathan Mace, Micha\l Bejda, Edison Gao, Wiktor Kuropatwa, Joe O’Neill, Kian Win Ong, Bill Schaller, Pingjia Shan, Brendan Viscomi, Vinod Venkataraman, Kaushik Veeraraghavan, and Yee Jiun Song. 2017.  Canopy: An End-to-End Performance Tracing And Analysis System. In _Proceedings of the 26th Symposium on Operating Systems Principles_ _(SOSP ’17)_. Association for Computing Machinery, New York, NY, USA, 34–50.  https://doi.org/10.1145/3132747.3132749
  * Krushevskaja and Sandler (2013) Darja Krushevskaja and Mark Sandler. 2013.  Understanding Latency Variations of Black Box Services. In _Proceedings of the 22nd International Conference on World Wide Web_ (Rio de Janeiro, Brazil) _(WWW ’13)_. Association for Computing Machinery, New York, NY, USA, 703–714.  https://doi.org/10.1145/2488388.2488450
  * Kullback and Leibler (1951) S. Kullback and R. A. Leibler. 1951.  On Information and Sufficiency.  _The Annals of Mathematical Statistics_ 22, 1 (1951), 79 – 86.  https://doi.org/10.1214/aoms/1177729694
  * Laaber and Leitner (2018) Christoph Laaber and Philipp Leitner. 2018.  An Evaluation of Open-Source Software Microbenchmark Suites for Continuous Performance Assessment. In _Proceedings of the 15th International Conference on Mining Software Repositories_ (Gothenburg, Sweden) _(MSR ’18)_. Association for Computing Machinery, New York, NY, USA, 119–130.  https://doi.org/10.1145/3196398.3196407
  * Li et al. (2021) Bowen Li, Xin Peng, Qilin Xiang, Hanzhang Wang, Tao Xie, Jun Sun, and Xuanzhe Liu. 2021.  Enjoy your observability: an industrial survey of microservice tracing and analysis.  _Empirical Software Engineering_ 27, 1 (2021), 25\.  https://doi.org/10.1007/s10664-021-10063-9
  * Luo et al. (2021) Shutian Luo, Huanle Xu, Chengzhi Lu, Kejiang Ye, Guoyao Xu, Liping Zhang, Yu Ding, Jian He, and Chengzhong Xu. 2021.  Characterizing Microservice Dependency and Performance: Alibaba Trace Analysis. In _Proceedings of the ACM Symposium on Cloud Computing_. 412–426.
  * Newman (2015) Sam Newman. 2015\.  _Building Microservices_ (1st ed.).  O’Reilly Media, Inc.
  * Pallets (2023) Pallets. 2023.  Web Development, one drop at time.  https://flask.palletsprojects.com/ ”Accessed 2023-01-21 19:52”.
  * Parker et al. (2020) Austin Parker, Daniel Spoonhower, Jonathan Mace, Ben Sigelman, and Rebecca Isaacs. 2020.  _Distributed Tracing in Practice: Instrumenting, Analyzing, and Debugging Microservices_.  O’Reilly Media, Incorporated.  tex.googlebooks: fgfIyAEACAAJ.
  * Rubin and Rinard (2016) Julia Rubin and Martin Rinard. 2016.  The Challenges of Staying Together While Moving Fast: An Exploratory Study. In _Proceedings of the 38th International Conference on Software Engineering_ _(ICSE ’16)_. Association for Computing Machinery, New York, NY, USA, 982–993.  https://doi.org/10.1145/2884781.2884871
  * Sambasivan et al. (2016) Raja R. Sambasivan, Ilari Shafer, Jonathan Mace, Benjamin H. Sigelman, Rodrigo Fonseca, and Gregory R. Ganger. 2016.  Principled Workflow-centric Tracing of Distributed Systems. In _Proceedings of the Seventh ACM Symposium on Cloud Computing_ _(SoCC ’16)_. ACM, New York, NY, USA, 401–414.  https://doi.org/10.1145/2987550.2987568
  * Sambasivan et al. (2013) R. R. Sambasivan, I. Shafer, M. L. Mazurek, and G. R. Ganger. 2013.  Visualizing Request-Flow Comparison to Aid Performance Diagnosis in Distributed Systems.  _IEEE Transactions on Visualization and Computer Graphics_ 19, 12 (2013), 2466–2475.  https://doi.org/10.1109/TVCG.2013.233
  * Sambasivan et al. (2011) Raja R. Sambasivan, Alice X. Zheng, Michael De Rosa, Elie Krevat, Spencer Whitman, Michael Stroucken, William Wang, Lianghong Xu, and Gregory R. Ganger. 2011\.  Diagnosing Performance Changes by Comparing Request Flows. In _Proceedings of the 8th USENIX Conference on Networked Systems Design and Implementation_ (Boston, MA) _(NSDI’11)_. USENIX Association, USA, 43–56.
  * Shannon (1948) C. E. Shannon. 1948\.  A mathematical theory of communication.  _The Bell System Technical Journal_ 27, 4 (1948), 623–656.  https://doi.org/10.1002/j.1538-7305.1948.tb00917.x
  * Sigelman et al. (2010) Benjamin H. Sigelman, Luiz André Barroso, Mike Burrows, Pat Stephenson, Manoj Plakal, Donald Beaver, Saul Jaspan, and Chandan Shanbhag. 2010.  _Dapper, a Large-Scale Distributed Systems Tracing Infrastructure_.  Technical Report. Google, Inc.
  * Sridharan (2017) Cindy Sridharan. 2017\.  Testing Microservices, the sane way.  https://copyconstruct.medium.com/testing-microservices-the-sane-way-9bb31d158c16 ”Accessed 2023-01-16 11:14”.
  * Technologies (2023) Uber Technologies. 2023\.  Jaeger: Open source, end-to-end distributed tracing.  https://www.jaegertracing.io/ ”Accessed 2023-01-15 14:10:59”.
  * Traini (2022) Luca Traini. 2022\.  Exploring Performance Assurance Practices and Challenges in Agile Software Development: An Ethnographic Study.  _Empirical Software Engineering_ 27, 3 (2022), 74\.  https://doi.org/10.1007/s10664-021-10069-3 ISBN: 1573-7616.
  * Traini and Cortellessa (2021) Luca Traini and Vittorio Cortellessa. 2021.  DeLag: Using Multi-Objective Optimization to Enhance the Detection of Latency Degradation Patterns in Service-based Systems.  https://doi.org/10.48550/ARXIV.2110.11155
  * Traini et al. (2022) Luca Traini, Vittorio Cortellessa, Daniele Di Pompeo, and Michele Tucci. 2022.  Towards effective assessment of steady state performance in Java software: are we there yet?  _Empirical Software Engineering_ 28, 1 (2022), 13\.  https://doi.org/10.1007/s10664-022-10247-x
  * Traini et al. (2021) Luca Traini, Daniele Di Pompeo, Michele Tucci, Bin Lin, Simone Scalabrino, Gabriele Bavota, Michele Lanza, Rocco Oliveto, and Vittorio Cortellessa. 2021.  How Software Refactoring Impacts Execution Time.  _ACM Trans. Softw. Eng. Methodol._ 31, 2, Article 25 (dec 2021), 23 pages.  https://doi.org/10.1145/3485136
  * Veeraraghavan et al. (2016) Kaushik Veeraraghavan, Justin Meza, David Chou, Wonho Kim, Sonia Margulis, Scott Michelson, Rajesh Nishtala, Daniel Obenshain, Dmitri Perelman, and Yee Jiun Song. 2016.  Kraken: Leveraging Live Traffic Tests to Identify and Resolve Resource Utilization Bottlenecks in Large Scale Web Services. In _12th USENIX Symposium on Operating Systems Design and Implementation (OSDI 16)_. USENIX Association, Savannah, GA, 635–651.
  * Zhang et al. (2022) Chenxi Zhang, Xin Peng, Chaofeng Sha, Ke Zhang, Zhenqing Fu, Xiya Wu, Qingwei Lin, and Dongmei Zhang. 2022\.  DeepTraLog: Trace-Log Combined Microservice Anomaly Detection through Graph-based Deep Learning. In _2022 IEEE/ACM 44th International Conference on Software Engineering (ICSE)_. 623–634.  https://doi.org/10.1145/3510003.3510180
  * Zhou et al. (2021) Xiang Zhou, Xin Peng, Tao Xie, Jun Sun, Chao Ji, Wenhai Li, and Dan Ding. 2021\.  Fault Analysis and Debugging of Microservice Systems: Industrial Survey, Benchmark System, and Empirical Study.  _IEEE Transactions on Software Engineering_ 47, 2 (2021), 243–260.  https://doi.org/10.1109/TSE.2018.2887384
  * Zhou et al. (2019) Xiang Zhou, Xin Peng, Tao Xie, Jun Sun, Chao Ji, Dewei Liu, Qilin Xiang, and Chuan He. 2019.  Latent Error Prediction and Fault Localization for Microservice Applications by Learning from System Trace Logs. In _Proceedings of the 2019 27th ACM Joint Meeting on European Software Engineering Conference and Symposium on the Foundations of Software Engineering_ (Tallinn, Estonia) _(ESEC/FSE 2019)_. Association for Computing Machinery, New York, NY, USA, 683–694.  https://doi.org/10.1145/3338906.3338961


††thanks<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
# Similar Scale-invariant Behaviors between Soft Gamma-ray Repeaters and An
Extreme Epoch from FRB 121102

Jun-Jie Wei  Purple Mountain Observatory, Chinese Academy of Sciences, Nanjing
210023, China  School of Astronomy and Space Sciences, University of Science
and Technology of China, Hefei 230026, China  Xue-Feng Wu  Purple Mountain
Observatory, Chinese Academy of Sciences, Nanjing 210023, China  School of
Astronomy and Space Sciences, University of Science and Technology of China,
Hefei 230026, China  Zi-Gao Dai  School of Astronomy and Space Sciences,
University of Science and Technology of China, Hefei 230026, China  School of
Astronomy and Space Science, Nanjing University, Nanjing 210023, China  Fa-Yin
Wang  School of Astronomy and Space Science, Nanjing University, Nanjing
210023, China  Pei Wang  CAS Key Laboratory of FAST, NAOC, Chinese Academy of
Sciences, Beijing 100101, China  Di Li  CAS Key Laboratory of FAST, NAOC,
Chinese Academy of Sciences, Beijing 100101, China  University of Chinese
Academy of Sciences, Beijing 100049, China  NAOC-UKZN Computational
Astrophysics Centre, University of KwaZulu-Natal, Durban 4000, South Africa
Bing Zhang  Department of Physics and Astronomy, University of Nevada, Las
Vegas, Las Vegas, NV 89154, USA

###### Abstract

In this work, we study the statistical properties of soft gamma-/hard X-ray
bursts from the soft gamma-ray repeater (SGR) 1806–20 and SGR J1935+2154 and
of radio bursts from the repeating fast radio burst (FRB) 121102. For SGRs, we
show that the probability density functions for the differences of fluences,
fluxes, and durations at different times have fat tails with a $q$-Gaussian
form. The $q$ values in the $q$-Gaussian distributions are approximately
steady and independent of the temporal interval scale adopted, implying a
scale-invariant structure of SGRs. These features indicate that SGR bursts may
be governed by a self-organizing criticality (SOC) process. The recent
discovery of a Galactic FRB associated with a hard X-ray burst from the
Galactic magnetar SGR J1935+2154 has established the magnetar origin of at
least some FRBs. Very recently, 1652 independent bursts from the repeating FRB
121102 have been detected by the Five-hundred-meter Aperture Spherical radio
Telescope (FAST). Here we also investigate the scale-invariant structure of
FRB 121102 based on the latest observation of FAST, and show that FRB 121102
and SGRs share similar statistical properties. Scale invariance in both FRB
121102 and SGRs can be well explained within the same physical framework of
fractal-diffusive SOC systems.

Magnetars, Soft gamma-ray repeaters, Radio transient sources, Radio bursts

##  I Introduction

Self-organized criticality (SOC) Bak et al. (1987); Aschwanden et al. (2016)
has been found in dynamical behaviors of astrophysics systems. We here examine
two types of transient phenomena, namely, soft gamma-ray repeaters (SGRs) and
fast radio bursts (FRBs), which have been shown to be associated with each
other in at least one case (Zhang, 2020).

Magnetars are highly magnetized neutron stars that exhibit dramatic
variability over a broad range of timescales Duncan and Thompson (1992);
Kouveliotou et al. (1998) (see also Kaspi and Beloborodov (2017) for a
review). Many magnetars were first observed as SGRs, i.e., sources of repeated
bursts of soft gamma-/hard X-rays, with typical durations of $\sim 0.1-1$ s
and peak luminosities of $\sim 10^{39}-10^{41}$ erg $\rm s^{-1}$. Three
several-minute-long giant flares with peak luminosities up to $\sim
10^{44}-10^{47}$ erg $\rm s^{-1}$ have also been observed from three different
SGRs (Hurley et al., 1999; Palmer et al., 2005). Both bursts and flares are
believed to be powered by the dissipation and decay of the ultrastrong
magnetic fields, either through neutron star crustquakes (Thompson and Duncan,
1995) or magnetic reconnection (Lyutikov, 2003).

Several works have found that the energy distribution of SGR bursts can be
well fitted by a power-law function Cheng et al. (1996); Göǧüş et al. (1999,
2000); Cheng et al. (2020). Power-law size distributions of extreme events can
be explained by SOC systems. The SOC subsystems will self-organize, owing to
some driving force, to a critical state, at which a small local disturbance
can produce an avalanche-like chain reaction of any size within the system
(Bak et al., 1987). The fundamental property of all SOC systems have in common
is the emergence of scale-free power-law size distributions (Aschwanden, 2011,
2012, 2015; Wang and Dai, 2013). SGRs, sharing statistical properties in terms
of the distributions of energy and waiting time with earthquakes, have been
regarded as SOC systems (Cheng et al., 1996). The earthquake-like behavior
also suggests that the energies of SGRs originate from starquakes of magnetars
(Duncan and Thompson, 1992; Thompson and Duncan, 1995).

Besides the power-law size distribution, another hallmark of SOC systems is
the scale invariance of the avalanche size differences (Caruso et al., 2007;
Wang et al., 2015). Caruso et al. Caruso et al. (2007) showed that the
probability density functions (PDFs) for the earthquake energy differences at
different times have fat tails with a $q$-Gaussian form, which could be well
explained by introducing a small-world topology on the Olami-Feder-Christensen
(OFC) model (Olami et al., 1992). In physics, the OFC model is one of the most
popular models displaying SOC. The $q$-Gaussian distribution was then
suggested as an important character for describing the presence of criticality
(Caruso et al., 2007). Moreover, the form of the PDF does not depend on the
time interval adopted for the earthquake energy difference, i.e., the $q$
values in $q$-Gaussian distributions keep nearly constant for different scale
intervals, which indicates that there is a scale-invariant structure in the
energy differences of earthquakes (Caruso et al., 2007; Wang et al., 2015).
Subsequently, Chang et al. Chang et al. (2017) studied 384 X-ray bursts
observed in three active episodes of SGR J1550–5418 and found that this SGR
shows similar scale-invariant behaviors with earthquakes, i.e., the PDFs of
the differences of fluences, peak fluxes, and durations exhibit a common
$q$-Gaussian distribution at different scale intervals, supporting the idea
that SGRs are powered by crustquakes of magnetars.

FRBs are intense millisecond-duration bursts of radio waves occurring in the
universe (Lorimer et al., 2007; Petroff et al., 2019; Cordes and Chatterjee,
2019). So far, more than 600 FRBs have been reported, and over two dozens of
them have been seen to repeat Spitler et al. (2016); CHIME/FRB Collaboration
et al. (2019). Although rapid developments have been made in the FRB research
field, the physical origin of FRBs remains mysterious (Platts et al., 2019;
Zhang, 2020; Xiao et al., 2021). Magnetars have been proposed as the likely
engine to power repeating FRBs (Popov and Postnov, 2010; Katz, 2016; Metzger
et al., 2017; Wang and Yu, 2017; Beloborodov, 2017; Kumar et al., 2017; Yang
and Zhang, 2018; Cheng et al., 2020; Wang et al., 2020). On 2020 April 28, one
FRB-like event was independently detected by the Canadian Hydrogen Intensity
Mapping Experiment (CHIME/FRB Collaboration et al., 2020) and the Survey for
Transient Astronomical Radio Emission 2 (Bochenek et al., 2020) in association
with a hard X-ray burst from the Galactic magnetar SGR J1935+2154 during its
active phase (Li et al., 2021a; Ridnaia et al., 2021; Tavani et al., 2021).
This discovery strongly supports that magnetar engines can produce at least
some extragalactic FRBs.

Recently, by analyzing 93 bursts from the repeating FRB 121102 by the
observation of the Green Bank Telescope at 4–8 GHz (Zhang et al., 2018), Lin &
Sang Lin and Sang (2020) found that FRB 121102 has the property of scale
invariance similar to that of SGR J1550–5418. During an extreme episode of
bursts, FRB 121102 produced 1652 detectable events within a time-span of 62
days (Li et al., 2021b). The peak burst rate reached 117 hr-1. Such a range of
cadence facilitate further investigation into the SOC structures of this
source. Another key feature of this burst set from FRB 121102 is the bimodal
energy distribution of the burst rate.

Among the 16 known SGRs (12 confirmed and 4 candidates) Olausen and Kaspi
(2014), SGR J1550–5418 is the only source to date that has been used to
investigate the scale-invariant property (Chang et al., 2017). So it is
interesting to know whether other SGRs (e.g., SGR 1806–20 and SGR J1935+2154)
share the same property. The magnetar research group at Sabancı University
presented their systematic temporal and broad-band spectral analysis of over
1,500 bursts from SGR J1550–5418, SGR 1900+14, and SGR 1806–20 observed with
the Rossi X-ray Timing Explorer (RXTE) Kırmızıbayrak et al. (2017). In these
public data, SGR 1806–20 has the largest burst sample (924 bursts), which will
be used in our following analysis. Additionally, since SGR J1935+2154 is the
first source that has been reported to power a Galactic FRB (CHIME/FRB
Collaboration et al., 2020; Bochenek et al., 2020), we will also investigate
the physical connection between SGR J1935+2154 and repeating FRBs.

In this paper, we investigate the scale-invariant behaviors of SGR 1806–20,
SGR J1935+2154, and FRB 121102. For SGRs, the PDFs of the differences of
fluences (or total counts), peak fluxes (or peak counts), and durations at
different time scales are shown in Section II. These PDFs can be well fitted
by a $q$-Gaussian function, and the functional form does not depend on the
scale interval adopted. In Section III, we compare the PDFs between SGRs and
repeating FRB 121102. The evidence that may shed light on the bimodal burst
energy distribution of FRB 121102 is discussed in Section IV. Lastly, a brief
summary and discussion are given in Section V.

##  II Scale-invariance in SGRs

###  II.1 Data

The magnetar research group at Sabancı University has constructed an online
database of magnetar bursts detected by the RXTE between 1996 and 2011, which
can be visited at http://magnetars.sabanciuniv.edu. This database contains 924
bursts from SGR 1806–20, representing the largest sample for a single
observation. The trigger time, total counts, peak counts, and duration of each
burst are available in the database.

For SGR J1935+2154, the bursts observed by the Gamma-ray Burst Monitor onboard
the Fermi satellite during the source’s six active episodes from 2014 to 2020
are used (Lin et al., 2020a, b). The number of bursts is 260. The available
information for each burst include trigger time, fluence, and duration, but
not peak flux.

Note that the adopted durations of all bursts from both SGR 1806–20 and SGR
J1935+2154 were determined by Bayesian block algorithm. We estimate the
waiting time between two adjacent bursts through $\Delta t=t_{i+1}-t_{i}$,
where $t_{i+1}$ and $t_{i}$ are the trigger times of the $(i+1)$-th and $i$-th
bursts, respectively. As both SGR 1806–20 and SGR J1935+2154 consist of
several isolated epochs, we discard the waiting time between the last burst of
a certain epoch and the first burst of the next epoch.

Table 1: The mean values of $q$ in the $q$-Gaussian distribution for earthquakes, SGRs, and FRB 121102 Phenomena | Energy range (erg) | Number |  $q$-fluence (energy) |  $q$-flux |  $q$-duration | References
---|---|---|---|---|---|---
earthquakes |  | 400,000 | $1.75\pm 0.15$ | – | – | Caruso et al. (2007)
SGR J1550–5418 |  | 384 | $2.41\pm 0.29$ | $2.40\pm 0.30$ | $2.06\pm 0.23$ | Chang et al. (2017)
SGR 1806–20 |  | 924 | $2.72\pm 0.05$ | $2.65\pm 0.05$ | $2.48\pm 0.05$ | This work
SGR J1935+2154 |  | 260 | $2.78\pm 0.12$ | – | $2.28\pm 0.15$ | This work
FRB 121102 | Overall bursts | 1652 | $2.51\pm 0.03$ | $2.50\pm 0.03$ | $1.42\pm 0.04$ | This work
| with $4\times 10^{36}<E<8\times 10^{39}$ |  |  |  |  |
FRB 121102 | Low-energy bursts | 1253 | $1.63\pm 0.09$ | $1.81\pm 0.05$ | $1.39\pm 0.06$ | This work
| with $4\times 10^{36}<E\leq 2\times 10^{38}$ |  |  |  |  |
| High-energy bursts | 399 | $2.21\pm 0.10$ | $2.34\pm 0.09$ | $1.50\pm 0.16$ | This work
| with $2\times 10^{38}<E<8\times 10^{39}$ |  |  |  |  |

Figure 1: Scale-invariant structure of SGR 1806–20. PDFs of the differences of
total counts (upper-left panel), peak counts (upper-right panel), and
durations (lower-left panel) for $n=1$ (red diamonds), $n=50$ (green circles),
and $n=100$ (blue triangles), and the best-fitting $q$-Gaussian distributions
(solid lines). Lower-right panel: the best-fitting $q$ values in the
$q$-Gaussian distribution as a function of $n$.

Figure 2: Scale-invariant structure of SGR J1935+2154. PDFs of the differences
of fluences (left panel) and durations (middle panel) for $n=1$ (red
diamonds), $n=50$ (green circles), and $n=100$ (blue triangles), and the best-
fitting $q$-Gaussian distributions (solid lines). Right panel: the best-
fitting $q$ values in the $q$-Gaussian distribution as a function of $n$.

###  II.2 PDFs of the Avalanche Size Differences

We are now in position to use the data of SGR 1806–20 and SGR J1935+2154 to
investigate the PDFs of the avalanche size differences. The difference between
two avalanche sizes is given by

$X_{n}=S_{i+n}-S_{i}\;,$ (1)

where $S_{i}$ is the size (fluence, peak flux, duration, or waiting time) of
the $i$-th burst in temporal order and $n$ (an integer) is the temporal
interval scale. In practice, $X_{n}$ is normalized to

$x_{n}=\frac{X_{n}}{\sigma_{X_{n}}}\;,$ (2)

where $\sigma_{X_{n}}$ is the standard deviation of $X_{n}$. We are interest
in the statistical distribution of $x_{n}$.

Figure 1 shows the statistical results of SGR 1806–20. In this plot, we
display the PDFs of the differences of total counts (upper-left panel), peak
counts (upper-right panel), and durations (lower-left panel) for $n=1$
(diamonds), $n=50$ (circles), and $n=100$ (triangles). Following Chang et al.
(2017), the data are binned based on the Freedman-Diaconis rule (Freedman and
Diaconis, 1981). One can see from Figure 1 that these PDFs $P(x_{n})$ have a
sharp peak and fat tails, which are different from Gaussian behaviors. The
sharp peak means that small size differences are most likely to happen, while
the fat tails suggest that there are rare but relatively large size
differences. The large fluctuations presented in the tails are caused by the
incompleteness given from the lack sampling of small magnitude events at the
global scale (Caruso et al., 2007). Additionally, the data points in Figure 1
are almost independent of the temporal interval scale $n$ adopted for the
avalanche size difference, indicating a common form of $P(x_{n})$. Here we use
the Tsallis $q$-Gaussian function (Tsallis, 1988; Tsallis et al., 1998)

$f(x_{n})=\alpha\left[1-\left(1-q\right)x_{n}^{2}/\beta\right]^{\frac{1}{1-q}}$
(3)

to fit $P(x_{n})$, where $\alpha$, $\beta$, and $q$ are free parameters. When
$q\rightarrow 1$, the $q$-Gaussian distribution reduces to the Gaussian
distribution with mean zero and standard deviation $\sigma=\sqrt{\beta/2}$.
Thus, $q\neq 1$ denotes a departure from Gaussian statistics.

The best-fitting parameters ($\alpha$, $\beta$, and $q$) are obtained by
minimizing the $\chi^{2}$ statistics,

$\chi^{2}=\sum_{i}\frac{\left[f(x_{n,i})-P(x_{n,i})\right]^{2}}{\sigma_{P,i}^{2}}\;,$
(4)

where $\sigma_{P,i}=\sqrt{N_{{\rm bin},i}}/N_{\rm tot}$ is the uncertainty of
the data point,111For clarity, the uncertainties of the data points are not
presented in the figure. with $N_{{\rm bin},i}$ is the event number in the
$i$-th bin and $N_{\rm tot}$ is the total number of $x_{n}$. We use the python
implementation, emcee (Foreman-Mackey et al., 2013), to apply the Markov Chain
Monte Carlo method to derive the best-fitting values and their corresponding
uncertainties for these parameters. In Figure 1, the red, green, and blue
smooth curves correspond to the best-fitting results for $n=1$, $n=50$, and
$n=100$, respectively. We may see that the PDFs of the differences of total
counts, peak counts, and durations are well fitted by the $q$-Gaussian
function, and that the three curves in each panel overlap almost completely.
However, the PDF of the waiting time differences can not be well fitted by
$q$-Gaussian (see also Ref. Chang et al. (2017)), so we do not exhibit it in
the figure. The waiting time shows different behaviors, which may be due to
the discontinuous observations of the telescope.

Furthermore, we compute the PDFs of the differences of total counts, peak
counts, and durations at different scale intervals $1\leq n\leq 100$, and fit
the PDFs with the $q$-Gaussian function. In the lower-right panel of Figure 1,
we plot the best-fitting $q$ values as a function of $n$. We find that the $q$
values are approximately invariant for different scale intervals $n$. The
property of the independence of $q$ values on $n$ is referred to as scale
invariance. The mean values of $q$ for total counts, peak counts, and duration
are $2.72\pm 0.05$, $2.65\pm 0.05$, and $2.48\pm 0.05$, respectively, which
are summarized in Table 1. Here the uncertainty denotes the standard deviation
of $q$. Interestingly, the $q$ values we find here are close to the values
derived from SGR J1550–5418 (Chang et al., 2017), which indicates that there
is a common scale-invariant property in SGRs.

Similar results can be obtained for the data of SGR J1935+2154, see Figure 2.
The $q$ values for this data also keep approximately steady for different
scale intervals $n$. The mean values of $q$ for fluence and duration are
$2.78\pm 0.12$ and $2.28\pm 0.14$, respectively. As shown in Table 1, these
values are well consistent with those of SGR J1550–5418 (Chang et al., 2017)
and SGR 1806–20, supporting a common scale-invariant structure of SGRs again.

##  III Comparison with FRB 121102

Owing to the association between FRB 200428 and an X-ray burst from the
Galactic magnetar SGR J1935+2154 (CHIME/FRB Collaboration et al., 2020;
Bochenek et al., 2020; Li et al., 2021a; Ridnaia et al., 2021; Tavani et al.,
2021), it is natural to consider whether radio bursts of repeating FRBs have a
similar scale-invariant behavior. In this section, we compare the PDFs between
SGRs and repeating FRBs.

Recently, Li et al. Li et al. (2021b) reported the detection of 1652
independent bursts in 59.5 hours spanning 62 days using the Five-hundred-meter
Aperture Spherical radio Telescope (FAST) Nan et al. (2011); Li et al. (2018)
at 1.05-1.45 GHz. Such a uniform sample prevent complex selection effects,
which can be introduced by different instruments and different frequency bands
used in observations. Here we focus on the PDFs for the differences of
energies, fluxes, and durations at different time scales. For a fixed $n$, we
fit the PDF with the $q$-Gaussian function and extract the best-fitting $q$
value. Then we vary $n$ and derive the $q$ values as a function of $n$.

Figure 3: Scale-invariant structure of FRB 121102. PDFs of the differences of
energies (upper-left panel), fluxes (upper-right panel) and durations (lower-
left panel) for $n=1$ (red diamonds), $n=50$ (green circles), and $n=100$
(blue triangles), and the best-fitting $q$-Gaussian distributions (solid
lines). Lower-right panel: the best-fitting $q$ values in the $q$-Gaussian
distribution as a function of $n$.

Figure 3 exhibits some examples of the fits for the sample of FRB 121102. In
this figure, we show the PDFs of the differences of energies (upper-left
panel), fluxes (upper-right panel), and durations (lower-left panel) for $n=1$
(diamonds), $n=50$ (circles), and $n=100$ (triangles). The solid curves stand
for the best-fitting results. One can see that the PDFs can be well fitted by
the $q$-Gaussian function, in good agreement with the results of Lin and Sang
(2020). Similar to the scenario of SGRs, we also find that the $q$-Gaussian
function gives a poor fit to the PDF of the FRB waiting time differences,
which is then not shown in the figure (see also Ref. Lin and Sang (2020)). To
better represent the scale-invariant property in FRB 121102, in the lower-
right panel of Figure 3 we plot the best-fitting $q$ values as a function of
the temporal interval scale $n$. We confirm that the $q$ values are almost
stable and independent of $n$. This implies that there are similar scale-
invariant properties among earthquakes, SGRs, and FRBs. We also list the
average $q$ values of energy ($2.51\pm 0.03$), flux ($2.50\pm 0.03$), and
duration ($1.42\pm 0.04$) for FRB 121102 in Table 1. We emphasize that the
$q$-Gaussian distribution is a generalization of the standard Gaussian
distribution and it reduces to a Gaussian distribution when $q\rightarrow 1$.
Since the average $q$ value of duration is relatively small (closing to $1$),
the peak of the PDF in the lower-left panel of Figure 3 is no longer very
sharp as other PDFs.

##  IV Evidence of a bimodal burst energy distribution of FRB 121102

The burst energy distribution of FRB 121102 can be adequately described as
bimodal (Li et al., 2021b). At the low energy end ($\sim E<2\times 10^{38}$
erg), a log-normal function can describe the distribution reasonably well.
While the distribution at the high energy end ($\sim E>2\times 10^{38}$ erg)
can be well fitted by a generalized Cauchy function. It is interesting to
investigate whether low-energy bursts and high-energy bursts of FRB 121102
have a similar scale-invariant behavior.

We divide the overall 1652 bursts into two components, i.e., 1253 low-energy
bursts with $E\leq 2\times 10^{38}$ erg and 399 high-energy bursts with
$E>2\times 10^{38}$ erg. Then the scale-invariant structure of the two
components are investigated, respectively. As shown in Table 1, the average
$q$ values of energy, peak flux, and duration for the low-energy component are
$q_{E}=1.63\pm 0.09$, $q_{P}=1.81\pm 0.05$, and $q_{T}=1.39\pm 0.06$. The
average $q$ values for the high-energy component are $q_{E}=2.21\pm 0.10$,
$q_{P}=2.34\pm 0.09$, and $q_{T}=1.50\pm 0.16$. Interestingly, the average $q$
value obtained from the PDFs of the energy differences of low-energy bursts is
compatible, within the standard deviations, to that one found from earthquakes
($q=1.75\pm 0.15$) Caruso et al. (2007). Consistent with theoretical models of
FRBs based on magnetar Popov and Postnov (2010); Lyubarsky (2014); Beloborodov
(2017); Metzger et al. (2019), the earthquake-like behavior here may indicate
that the low-energy bursts of FRB 121102 could originate from starquakes of a
magnetar. More interestingly, the average $q$ values of energy, peak flux, and
duration for the high-energy component of FRB 121102 are roughly consistent
with those of SGRs. We next show that these $q$ values can be well understood
within the same statistical framework of a SOC system.

For a fractal-diffusive avalanche model, quantitative values for the size
distributions of SOC parameters (length scales $L$, durations $T$, peak fluxes
$P$, and fluences or energies $E$) are derived from first principles, using
the scale-free probability conjecture, $N(L)dL\propto L^{-S}dL$, for three
Euclidean dimensions $S=1$, 2, and 3. The analytical model predicts the
indices of size distributions for $E$, $T$, and $P$ as (Aschwanden, 2012,
2014)

$\displaystyle\alpha_{E}$ $\displaystyle=$ $\displaystyle 1+(S-1)/(D_{S}+2)$
$\displaystyle\alpha_{T}$ $\displaystyle=$ $\displaystyle(1+S)/2\;,$ (5)
$\displaystyle\alpha_{P}$ $\displaystyle=$ $\displaystyle 2-1/S$

where $D_{S}\approx(1+S)/2$ is the mean fractal dimension. As mentioned above,
Caruso et al. Caruso et al. (2007) presented an analysis method to explain SOC
behavior in the limited number of earthquakes by making use of the return
distributions (i.e., distributions of the avalanche size differences at
different times). They obtained the first strong evidence that the return
distributions appear to have the shape of $q$-Gaussians, standard
distributions arising naturally in nonextensive statistical mechanics
(Tsallis, 1988; Tsallis et al., 1998). Under the assumption that there is no
correlation between the sizes of two events, an exact relation between the
index $\alpha$ of the avalanche size distribution and the $q$ values of the
appropriate $q$-Gaussian has been obtained as (Caruso et al., 2007; Celikoglu
et al., 2010)

$q=\frac{\alpha+2}{\alpha}\;,$ (6)

which is important because it makes the $q$ parameter determined a priori from
one of the well-known indices of the system.

In the framework of the fractal-diffusive SOC model, we obtain the theoretical
$q$ values, i.e., $q_{E}\approx 2.33$, $q_{P}\approx 2.20$, and $q_{T}=2.0$,
according to Equations (IV) and (6) by taking the three-dimensional Euclidean
space $S=3$. The observed $q$ values of SGR J1550–5418, SGR 1806–20, and SGR
J1935+2154 at $1\sigma$ confidence level are in the range of $q_{E}\approx
2.12-2.9$, $q_{P}\approx 2.37-2.70$, and $q_{T}\approx 1.83-2.53$, which agree
well with the predictions of the fractal-diffusive SOC model. We also note
that the derived $q_{E}$ and $q_{P}$ of the high-energy component of FRB
121102 ($q_{E}=2.21\pm 0.10$ and $q_{P}=2.34\pm 0.09$) are in good agreement
with the predictions of the fractal-diffusive SOC model, and the derived
$q_{T}$ ($1.50\pm 0.16$) is also roughly consistent with the model prediction
at $3.1\sigma$ confidence level. Therefore, SGRs and high-energy bursts of FRB
121102 share similar scale-invariant behavior, suggesting that they can be
explained within the same statistical framework of fractal-diffusive SOC
systems with the spatial dimension $S=3$.

The observed durations of the repeating bursts of FRB 121102 are directly used
for the analysis of scale invariance. As is well known, the observed pulse
duration of FRB is generally broadened by the instrumental and astrophysical
sources (Cordes and McLaughlin, 2003). The instrumental pulse broadening
includes the sampling time-scale and the intrachannel dispersion smearing. The
temporal resolution can not be better than the sampling time-scale, which
broadens the pulse. The smearing is caused by intrachannel dispersion. The
astrophysical pulse broadening is the scattering process of the radio emission
in inhomogeneous plasma. Therefore, the observed duration does not directly
represent the intrinsic duration. The intrinsic durations of the repeating
bursts, in principle, should be used for our purpose. In practice, it is
difficult to separate the intrinsic duration from the broadening components.
But note that the distribution of the observed duration differences at
different time intervals is studied in this work. Pulse broadening caused by
smearing and scattering can be approximately deducted from subtracting two
observed durations. That is, the difference between two observed durations can
be approximated as the difference between two intrinsic ones.

##  V Summary and Discussion

Previous studies have shown that the PDFs for the earthquake energy
differences at different time intervals have fat tails with a $q$-Gaussian
shape (Caruso et al., 2007; Wang et al., 2015). Another remarkable feature is
that the $q$ values in the $q$-Gaussian distributions are approximately equal
and independent of the temporal interval scale adopted, indicating a scale-
invariant structure in the energy differences of earthquakes. These
statistical features can be well explained within a dissipative SOC scenario
taking into account long-range interactions (Caruso et al., 2007). That is,
scale invariance is an important character representing the system approaching
to a critical state. Recently, Chang et al. Chang et al. (2017) found that SGR
J1550–5418 has the property of scale invariance similar to that of
earthquakes. However, it is unclear whether other SGRs share the same scale-
invariant structure as SGR J1550–5418 and earthquakes.

In this work, we present statistics of soft gamma-/hard X-ray bursts from SGR
1806–20 and SGR J1935+2154. We find that the two SGRs share earthquake-like
behaviors in terms of the PDFs of the fluence differences, consistent with the
results of Chang et al. (2017). The PDFs of the fluence differences can be
well fitted by a $q$-Gaussian function and the $q$ values keep nearly constant
for different scale intervals. The earthquake-like behaviors support that
there is a common scale-invariant structure in SGRs and that systems
responsible for SGR bursts are in an SOC state. Moreover, we show that the
burst fluxes and durations of the SGRs share similar scale-invariant behaviors
with fluences, i.e., the PDFs of the differences of fluxes and durations at
different time intervals also well follow $q$-Gaussian distributions.

Given the association between an FRB-like event and an X-ray burst from the
Galactic magnetar SGR J1935+2154, we also investigate the scale-invariant
structure of the repeating FRB 121102. For FRB 121102, we show that the PDFs
of the differences of energies, fluxes, and durations also exhibit
$q$-Gaussian distributions, with steady $q$ values independent of the time
scale. These properties are very similar to those of earthquakes and SGRs, and
thus all can be attributed to a SOC process. Due to the fact that the burst
energy distribution of FRB 121102 is bimodal, we also investigate the scale-
invariant behaviors of low-energy bursts and high-energy bursts of FRB 121102.
We find that the $q_{E}$, $q_{P}$, and $q_{T}$ values of low-energy bursts are
different from those of high-energy bursts, which indicate that the low-energy
component and the high-energy component may have different physical origins.
Interestingly, the average $q$ value obtained from the PDFs of the energy
differences of low-energy bursts is close to the $q$ value found from
earthquakes. This implies that there may be some similarities between the
energy origins of low-energy bursts of FRB 121102 and earthquakes. More
interestingly, the average $q$ values of energy, peak flux, and duration for
the high-energy component of FRB 121102 are roughly consistent with those of
SGRs. These $q$ values can be well explained within the same statistical
framework of fractal-diffusive SOC systems with the spatial dimension $S=3$.
In the future, much more repeating bursts of FRBs will be detected. The
physical connection between repeating FRBs and SGRs can be further
investigated.

###### Acknowledgements.

This work is partially supported by the National Natural Science Foundation of
China (grant Nos. 11988101, 11725314, U1831122, 12041306, U1831207), the Youth
Innovation Promotion Association (2017366), the Key Research Program of
Frontier Sciences (grant No. ZDBS-LY-7014) of Chinese Academy of Sciences, and
the Major Science and Technology Project of Qinghai Province (2019-ZJ-A10).
ZGD was supported by the National Key Research and Development Program of
China (grant No. 2017YFA0402600), the National SKA Program of China (grant No.
2020SKA0120300), and the National Natural Science Foundation of China (grant
No. 11833003).

## References

  * Bak et al. (1987) P. Bak, C. Tang, and K. Wiesenfeld, Phys. Rev. Lett. 59, 381 (1987).
  * Aschwanden et al. (2016) M. J. Aschwanden, N. B. Crosby, M. Dimitropoulou, M. K. Georgoulis, S. Hergarten, J. McAteer, A. V. Milovanov, S. Mineshige, L. Morales, N. Nishizuka, et al., Space Sci. Rev. 198, 47 (2016), eprint 1403.6528.
  * Zhang (2020) B. Zhang, Nature (London) 587, 45 (2020), eprint 2011.03500.
  * Duncan and Thompson (1992) R. C. Duncan and C. Thompson, Astrophys. J. Lett. 392, L9 (1992).
  * Kouveliotou et al. (1998) C. Kouveliotou, S. Dieters, T. Strohmayer, J. van Paradijs, G. J. Fishman, C. A. Meegan, K. Hurley, J. Kommers, I. Smith, D. Frail, et al., Nature (London) 393, 235 (1998).
  * Kaspi and Beloborodov (2017) V. M. Kaspi and A. M. Beloborodov, Ann. Rev. Astron. Astrophys. 55, 261 (2017), eprint 1703.00068.
  * Hurley et al. (1999) K. Hurley, T. Cline, E. Mazets, S. Barthelmy, P. Butterworth, F. Marshall, D. Palmer, R. Aptekar, S. Golenetskii, V. Il’Inskii, et al., Nature (London) 397, 41 (1999), eprint astro-ph/9811443.
  * Palmer et al. (2005) D. M. Palmer, S. Barthelmy, N. Gehrels, R. M. Kippen, T. Cayton, C. Kouveliotou, D. Eichler, R. A. M. J. Wijers, P. M. Woods, J. Granot, et al., Nature (London) 434, 1107 (2005), eprint astro-ph/0503030.
  * Thompson and Duncan (1995) C. Thompson and R. C. Duncan, Mon. Not. R. Astron. Soc. 275, 255 (1995).
  * Lyutikov (2003) M. Lyutikov, Mon. Not. R. Astron. Soc. 346, 540 (2003), eprint astro-ph/0303384.
  * Cheng et al. (1996) B. Cheng, R. I. Epstein, R. A. Guyer, and A. C. Young, Nature (London) 382, 518 (1996).
  * Göǧüş et al. (1999) E. Göǧüş, P. M. Woods, C. Kouveliotou, J. van Paradijs, M. S. Briggs, R. C. Duncan, and C. Thompson, Astrophys. J. Lett. 526, L93 (1999), eprint astro-ph/9910062.
  * Göǧüş et al. (2000) E. Göǧüş, P. M. Woods, C. Kouveliotou, J. van Paradijs, M. S. Briggs, R. C. Duncan, and C. Thompson, Astrophys. J. Lett. 532, L121 (2000), eprint astro-ph/0002181.
  * Cheng et al. (2020) Y. Cheng, G. Q. Zhang, and F. Y. Wang, Mon. Not. R. Astron. Soc. 491, 1498 (2020), eprint 1910.14201.
  * Aschwanden (2011) M. J. Aschwanden, _Self-Organized Criticality in Astrophysics_ (2011).
  * Aschwanden (2012) M. J. Aschwanden, Astron. & Astrophys. 539, A2 (2012), eprint 1112.4859.
  * Aschwanden (2015) M. J. Aschwanden, Astrophys. J. 814, 19 (2015), eprint 1510.01987.
  * Wang and Dai (2013) F. Y. Wang and Z. G. Dai, Nature Physics 9, 465 (2013), eprint 1308.1253.
  * Caruso et al. (2007) F. Caruso, A. Pluchino, V. Latora, S. Vinciguerra, and A. Rapisarda, Phys. Rev. E 75, 055101 (2007), eprint cond-mat/0606118.
  * Wang et al. (2015) P. Wang, Z. Chang, H. Wang, and H. Lu, European Physical Journal B 88, 206 (2015).
  * Olami et al. (1992) Z. Olami, H. J. S. Feder, and K. Christensen, Phys. Rev. Lett. 68, 1244 (1992).
  * Chang et al. (2017) Z. Chang, H.-N. Lin, Y. Sang, and P. Wang, Chinese Physics C 41, 065104 (2017), eprint 1702.05842.
  * Lorimer et al. (2007) D. R. Lorimer, M. Bailes, M. A. McLaughlin, D. J. Narkevic, and F. Crawford, Science 318, 777 (2007), eprint 0709.4301.
  * Petroff et al. (2019) E. Petroff, J. W. T. Hessels, and D. R. Lorimer, Astron. Astrophys. Rev. 27, 4 (2019), eprint 1904.07947.
  * Cordes and Chatterjee (2019) J. M. Cordes and S. Chatterjee, Ann. Rev. Astron. Astrophys. 57, 417 (2019), eprint 1906.05878.
  * Spitler et al. (2016) L. G. Spitler, P. Scholz, J. W. T. Hessels, S. Bogdanov, A. Brazier, F. Camilo, S. Chatterjee, J. M. Cordes, F. Crawford, J. Deneva, et al., Nature (London) 531, 202 (2016), eprint 1603.00581.
  * CHIME/FRB Collaboration et al. (2019) CHIME/FRB Collaboration, B. C. Andersen, K. Bandura, M. Bhardwaj, P. Boubel, M. M. Boyce, P. J. Boyle, C. Brar, T. Cassanelli, P. Chawla, et al., Astrophys. J. Lett. 885, L24 (2019), eprint 1908.03507.
  * Platts et al. (2019) E. Platts, A. Weltman, A. Walters, S. P. Tendulkar, J. E. B. Gordin, and S. Kandhai, Physics Reports 821, 1 (2019), eprint 1810.05836.
  * Xiao et al. (2021) D. Xiao, F. Wang, and Z. Dai, Science China Physics, Mechanics, and Astronomy 64, 249501 (2021), eprint 2101.04907.
  * Popov and Postnov (2010) S. B. Popov and K. A. Postnov, in _Evolution of Cosmic Objects through their Physical Activity_ , edited by H. A. Harutyunian, A. M. Mickaelian, and Y. Terzian (2010), pp. 129–132, eprint 0710.2006.
  * Katz (2016) J. I. Katz, Astrophys. J. 826, 226 (2016), eprint 1512.04503.
  * Metzger et al. (2017) B. D. Metzger, E. Berger, and B. Margalit, Astrophys. J. 841, 14 (2017), eprint 1701.02370.
  * Wang and Yu (2017) F. Y. Wang and H. Yu, JCAP 2017, 023 (2017), eprint 1604.08676.
  * Beloborodov (2017) A. M. Beloborodov, Astrophys. J. Lett. 843, L26 (2017), eprint 1702.08644.
  * Kumar et al. (2017) P. Kumar, W. Lu, and M. Bhattacharya, Mon. Not. R. Astron. Soc. 468, 2726 (2017), eprint 1703.06139.
  * Yang and Zhang (2018) Y.-P. Yang and B. Zhang, Astrophys. J. 868, 31 (2018), eprint 1712.02702.
  * Wang et al. (2020) F. Y. Wang, Y. Y. Wang, Y.-P. Yang, Y. W. Yu, Z. Y. Zuo, and Z. G. Dai, Astrophys. J. 891, 72 (2020), eprint 2002.03507.
  * CHIME/FRB Collaboration et al. (2020) CHIME/FRB Collaboration, B. C. Andersen, K. M. Bandura, M. Bhardwaj, A. Bij, M. M. Boyce, P. J. Boyle, C. Brar, T. Cassanelli, P. Chawla, et al., Nature (London) 587, 54 (2020), eprint 2005.10324.
  * Bochenek et al. (2020) C. D. Bochenek, V. Ravi, K. V. Belov, G. Hallinan, J. Kocz, S. R. Kulkarni, and D. L. McKenna, Nature (London) 587, 59 (2020), eprint 2005.10828.
  * Li et al. (2021a) C. K. Li, L. Lin, S. L. Xiong, M. Y. Ge, X. B. Li, T. P. Li, F. J. Lu, S. N. Zhang, Y. L. Tuo, Y. Nang, et al., Nature Astronomy 5, 378 (2021a), eprint 2005.11071.
  * Ridnaia et al. (2021) A. Ridnaia, D. Svinkin, D. Frederiks, A. Bykov, S. Popov, R. Aptekar, S. Golenetskii, A. Lysenko, A. Tsvetkova, M. Ulanov, et al., Nature Astronomy 5, 372 (2021), eprint 2005.11178.
  * Tavani et al. (2021) M. Tavani, C. Casentini, A. Ursi, F. Verrecchia, A. Addis, L. A. Antonelli, A. Argan, G. Barbiellini, L. Baroncelli, G. Bernardi, et al., Nature Astronomy 5, 401 (2021), eprint 2005.12164.
  * Zhang et al. (2018) Y. G. Zhang, V. Gajjar, G. Foster, A. Siemion, J. Cordes, C. Law, and Y. Wang, Astrophys. J. 866, 149 (2018), eprint 1809.03043.
  * Lin and Sang (2020) H.-N. Lin and Y. Sang, Mon. Not. R. Astron. Soc. 491, 2156 (2020), eprint 1912.01191.
  * Li et al. (2021b) D. Li, P. Wang, W. W. Zhu, B. Zhang, X. X. Zhang, R. Duan, Y. K. Zhang, Y. Feng, N. Y. Tang, S. Chatterjee, et al., arXiv e-prints arXiv:2107.08205 (2021b), eprint 2107.08205.
  * Olausen and Kaspi (2014) S. A. Olausen and V. M. Kaspi, Astrophys. J. Suppl. Ser. 212, 6 (2014), eprint 1309.4167.
  * Kırmızıbayrak et al. (2017) D. Kırmızıbayrak, S. Şaşmaz Muş, Y. Kaneko, and E. Göğüş, Astrophys. J. Suppl. Ser. 232, 17 (2017), eprint 1708.08067.
  * Lin et al. (2020a) L. Lin, E. Göğüş, O. J. Roberts, C. Kouveliotou, Y. Kaneko, A. J. van der Horst, and G. Younes, Astrophys. J. 893, 156 (2020a), eprint 2003.10582.
  * Lin et al. (2020b) L. Lin, E. Göğüş, O. J. Roberts, M. G. Baring, C. Kouveliotou, Y. Kaneko, A. J. van der Horst, and G. Younes, Astrophys. J. Lett. 902, L43 (2020b), eprint 2010.02767.
  * Freedman and Diaconis (1981) D. Freedman and P. Diaconis, Probability Theory and Related Fields 57, 453 (1981).
  * Tsallis (1988) C. Tsallis, Journal of Statistical Physics 52, 479 (1988).
  * Tsallis et al. (1998) C. Tsallis, R. Mendes, and A. R. Plastino, Physica A Statistical Mechanics and its Applications 261, 534 (1998).
  * Foreman-Mackey et al. (2013) D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, Pub. Astro. Soc. Pacific 125, 306 (2013), eprint 1202.3665.
  * Nan et al. (2011) R. Nan, D. Li, C. Jin, Q. Wang, L. Zhu, W. Zhu, H. Zhang, Y. Yue, and L. Qian, International Journal of Modern Physics D 20, 989 (2011), eprint 1105.3794.
  * Li et al. (2018) D. Li, P. Wang, L. Qian, M. Krco, P. Jiang, Y. Yue, C. Jin, Y. Zhu, Z. Pan, R. Nan, et al., IEEE Microwave Magazine 19, 112 (2018), eprint 1802.03709.
  * Lyubarsky (2014) Y. Lyubarsky, Mon. Not. R. Astron. Soc. 442, L9 (2014), eprint 1401.6674.
  * Metzger et al. (2019) B. D. Metzger, B. Margalit, and L. Sironi, Mon. Not. R. Astron. Soc. 485, 4091 (2019), eprint 1902.01866.
  * Aschwanden (2014) M. J. Aschwanden, Astrophys. J. 782, 54 (2014), eprint 1310.4191.
  * Celikoglu et al. (2010) A. Celikoglu, U. Tirnakli, and S. M. D. Queirós, Phys. Rev. E 82, 021124 (2010), eprint 1003.0361.
  * Cordes and McLaughlin (2003) J. M. Cordes and M. A. McLaughlin, Astrophys. J. 596, 1142 (2003), eprint astro-ph/0304364.


# Spatio-temporal Structure of Excitation and Inhibition Emerges in Spiking
Neural Networks with and without Biologically Plausible Constraints

Balázs Mészáros 1,∗, James Knight 1 and Thomas Nowotny 1

###### Abstract

## 1

We present a Spiking Neural Network (SNN) model that incorporates learnable
synaptic delays using Dilated Convolution with Learnable Spacings (DCLS). We
train this model on the Raw Heidelberg Digits keyword spotting benchmark using
Backpropagation Through Time with surrogate gradients. Analysing the spatio-
temporal structure of synaptic interactions in the network we observe that
after training excitation and inhibition are grouped together both in space
and time. To further enhance the efficiency and biological realism of our
model, we implemented a dynamic pruning strategy that combines DEEP R for
connection removal and RigL for connection reintroduction, ensuring that the
network maintains optimal connectivity throughout training. Additionally, we
incorporated Dale’s Principle, enforcing each neuron to be exclusively
excitatory or inhibitory – aligning our model closer to biological neural
networks. We observed that, after training, the spatio-temporal patterns of
excitation and inhibition appeared in the more biologically plausible model as
well. Our research demonstrates the potential of integrating learnable delays,
dynamic pruning, and biological constraints to develop efficient SNN models
for temporal data processing. Furthermore, our results enhance the
understanding of spatio-temporal dynamics in SNNs – suggesting that the
spatio-temporal features which emerge from training are robust to both pruning
and rewiring processes – providing a solid foundation for future work in
neuromorphic computing applications.

##  2 Keywords:

Spiking Neural Network, Delay Learning, Dynamic Pruning, Receptive Field,
Sparse Connectivity

##  3 Introduction

Spiking Neural Networks (SNNs) are the third generation of artificial neural
networks (Maass, 1997), inspired by the functioning of biological neurons.
Unlike traditional neural networks, which are stateless and process
information through continuous activation values, neurons in SNNs are stateful
and communicate sparsely via binary spikes, mimicking the electrical impulses
observed in biological neurons. This enables SNNs to efficiently process
temporal information, garnering significant interest for tasks involving
sequential data processing. SNNs have demonstrated potential in various
fields, including neuromorphic computing (Furber et al., 2014; Davies et al.,
2018; Merolla et al., 2014), due to their capabilities for low-power and real-
time operation. This makes them a compelling basis for next-generation
intelligent systems.

Recently, learnable synaptic delays (between the emission of a spike and its
arrival at the post-synaptic neuron) have been suggested as a means of
improving the spatio-temporal information processing of SNNs (Hammouamri et
al., 2023). These delays, which can vary across connections, allow neurons to
perform coincidence detections across longer time intervals, enhancing SNN’s
ability to process temporal information (Izhikevich, 2006). In biological
systems, delays encompass axonal, synaptic, and dendritic components and are
modified by processes like myelination to facilitate learning (Bengtsson et
al., 2005) and coincidence detection (Seidl et al., 2010). Furthermore,
neuromorphic hardware such as Intel’s Loihi (Davies et al., 2018) and
SpiNNaker (Furber et al., 2014) incorporate programmable synaptic delays,
enabling SNNs with trained delays to be efficiently deployed for real-time
data processing.

Methods to decrease the number of network parameters were proposed more than
40 years ago (LeCun et al., 1989; Hassibi et al., 1993), and in recent years
overparametrisation in deep neural networks has become a widely acknowledged
problem (Ba and Caruana, 2014). Sparse neural networks significantly reduce
memory usage and energy consumption and are particularly necessary on
neuromorphic systems, which typically have a maximum fan-in per neuron
(Schemmel et al., 2010; Merolla et al., 2014) or limited memory available for
connectivity (Furber et al., 2014; Davies et al., 2018). In order to
efficiently deploy neural networks at the edge, unnecessary connections are
often removed using pruning but, traditional pruning methods often require
training a dense model first (Han et al., 2015), meaning that the size of
models is still limited by the training cost of the original dense model. In
contrast, biological systems dynamically rewire synaptic connections during
learning, suggesting that dynamic pruning and rewiring can enhance neural
network performance and efficiency. Techniques like DEEP R (Bellec et al.,
2018) and RigL (Evci et al., 2021), which integrate pruning and connection
reintroduction during training, aim to emulate this biological efficiency,
potentially offering superior sparsity and accuracy with fewer floating-point
operations (FLOPs).

In this paper, we present our analysis of a spiking neural network trained in
a supervised fashion on the Raw Heidelberg Digits benchmark dataset (Cramer et
al., 2022). We analyze the learnt parameters of the network, then enforce
biological constraints during training (sparse connectivity and Dale’s
principle), and conduct the same analysis on our new, more biologically
plausible network.

##  4 Methods

We consider an SNN with Leaky-Integrate-and-Fire (LIF) neurons,

$\tau\frac{du_{i}^{(l)}}{dt}=-u_{i}^{(l)}(t)+I_{i}^{(l)}(t)$ (1)

where, $u_{i}^{(l)}$ is the membrane potential of the $i$-th neuron in layer
$l$, $\tau=10.05$ is the membrane time constant and $I_{i}^{(l)}(t)$ is the
input current. Neurons emit a spike when their membrane potential reaches the
threshold $\vartheta=1$, yielding the spike train $S_{j}^{(l-1)}$ in layer
$l-1$,

$S_{j}^{(l-1)}(t)=\Theta(u_{j}^{(l-1)}(t)-\vartheta)$ (2)

When a spike is emitted the membrane potential is reset to 0. For our
numerical experiments, we used a linear Euler integration method with timestep
$10ms$, yielding the following discrete update equation

$u_{i}^{(l)}[t]=\left(1-\frac{1}{\tau}\right)u_{i}^{(l)}[t-1]+I_{i}^{(l)}[t-1]$
(3)

Because the Heaviside step function $\Theta$ is non-differentiable, like many
others, we replace it with a $\arctan$ surrogate gradient during the backward
pass of our training (Neftci et al., 2019). To implement the synaptic delay
training, the input current $I_{i}^{(l)}(t)$ is calculated using Dilated
Convolutions with Learnable Spacings (DCLS) (Hammouamri et al., 2023), i.e.
convolving the spike train $S_{j}^{(l-1)}$ from layer $l-1$ with the 1D kernel

$k_{ij}^{(l)}[n]=\frac{w_{ij}^{(l)}}{c}\exp\left(-\frac{1}{2}\left(\frac{n-T_{d}+d_{ij}^{(l)}+1}{\sigma}\right)^{2}\right),\quad
n=0,\ldots,T_{d}$ (4)

where $T_{d}=25$ is the maximum delay, $d_{ij}^{(l)}$ is the synaptic delay
from neuron $j$ to neuron $i$ in layer $l$, $c$ is a normalization term so
that $\sum_{n=0}^{T_{d}}k_{ij}^{(l)}[n]=w_{ij}^{(l)}$ and $\sigma=12.5$ is the
standard deviation of the delay kernel, which is decreased during training. As
the kernel $k_{ij}^{(l)}$ slides through the spike train $S_{j}^{(l-1)}$ at
each timestep $t$ the kernel will have access to spikes in the range of
{$t-T_{d},\ldots,t$}. For more details, we refer the reader to the original
publication (Hammouamri et al., 2023).

For dynamic pruning, we combined two methods: DEEP R (Bellec et al., 2018) for
dropping connections and RigL (Evci et al., 2021) for introducing them. In
DEEP R the synaptic weights are defined as
$w_{ij}^{l}=s_{ij}^{(l)}\text{max}(\theta_{ij}^{(l)},0)$, where $s_{ij}^{(l)}$
is a constant sign value, and $\theta_{ij}^{(l)}$ is the parameter trained
with gradient descent. If $\theta_{ij}^{(l)}$ is not positive, the connection
is considered dormant. We use $L1$ regularization with a strength of $0.1$ to
encourage pruning of unnecessary connections. Although the original DEEP R
method adds noise to the gradients to induce stochasticity, following Bellec
et al. (2020), we omitted this in our experiments. DEEP R maintains a fixed
number of synapses by randomly reactivating synapses throughout training.

RigL introduces synapses based on gradient analysis rather than randomly. In
each iteration, we select the $k$ highest gradients from the inactive
connections:

$\operatornamewithlimits{ArgTopK}_{\theta^{(l)}_{ij}\leq
0}\left(\frac{dL}{d\theta^{(l)}_{ij}},k\right)$ (5)

Since we make the assumption that active connections are positive, the
original RigL method needs to be modified slightly. Instead of picking
synapses to introduce based on the _absolute_ gradient value, we simply pick
based on the highest gradient values, since a negative gradient implies that
gradient descent wants to keep the connection inactive. Because, generally
speaking, larger layers will have more connections, following Evci et al.
(2021), we adjust the target sparsity to be higher (fewer connections) in
larger layers and lower (more connections) in smaller layers, using the Erdős-
Rényi formula (Renyi, 1959):

$p^{(l)}=p\left(1-\frac{n^{(l-1)}+n^{(l)}}{n^{(l-1)}n^{(l)}}\right)$ (6)

where $p$ is the chosen hyperparameter for sparsity (ranging from
{$0.5,\ldots,0.984375$} in our simulations), $n^{(l)}$ is the number of
neurons in layer $l$ and $p^{(l)}$ is the actual enforced sparsity.

Most neural network models do not adhere to Dale’s principle, which states
that neurons are either exclusively excitatory or inhibitory. Computationally,
this means that the signs for each element in a given column of the weight
matrix are the same. Since with DEEP R we have to generate the sign matrix
before training, we can conveniently apply the constraint by generating a
random sign vector $s$ and broadcasting it into a matrix. Unless stated
otherwise, all of our results apply to networks that adhere to Dale’s
principle.

We extended the PyTorch implementation by Hammouamri et al. (2023) for our
experiments. Our architecture included one hidden layer with 256 neurons,
dropout layers with a probability of 0.4, batch normalization, and all layers
had delays in the range of $({0,\ldots,25})$. We used a voltage readout and
cross-entropy loss for training. Our model was trained with the Adam
optimizer, using a OneCycle scheduler for weights and a Cosine Annealing
scheduler for delays.

Figure 1: Toy example of how the “spatio-temporal receptive field” figures
were generated. Note: in this example, we do not adhere to Dale’s law (the
second input neuron has both excitatory and inhibitory connections). Figure 2:
Spatio-temporal receptive fields in the non-constrained network. Red pixels
correspond to excitatory and blue pixels to inhibitory connections. Pixel
transparency corresponds to the synaptic strength. Clear spatio-temporal
patterns of excitation and inhibition emerge in each output neuron after
training the network. Figure 3:  Receptive fields of the sparse Dalean
network. The network was trained with 87.5% sparsity, before the Erdős-Rényi
adjustment. The spatio-temporal grouping of excitation and inhibition seems to
emerge in the more biologically plausible network as well.

##  5 Results

Figure 4: Classification accuracy with various architectures. The baseline
results are from a recurrent spiking neural network (Zenke and Vogels,
2021).(A) Performance comparison based on the number of synapses (x-axis) and
best results on the RawHD test set (y-axis). Results for the dense network
come from single simulations, while sparse Dalean network results are averaged
over three simulations due to initialization sensitivity. The points are the
means and the error bars are the standard deviations. The leftmost example is
not applicable for the sparse network as it represents a fully connected
network. As the number of synapses are decreased, sparse connectivity becomes
beneficial, even with the constraint to adhere to Dale’s principle. (B)
Comparing the effects of structure learning and delay learning. Learning the
structure does not make a huge difference with lower sparsity, but as it
increases the benefit becomes clear. This might not be that surprising, since
these methods were built to train highly sparse models. Delay learning
improves the results with a fixed connectivity matrix, as it was also showed
in (Hammouamri et al., 2023). The benefits of delay learning are not obvious
when the structure is learnt.

In our experiments, we trained our models with the best hyperparameters used
by Hammouamri et al. (2023) for training on the Spiking Heidelberg Digits
(SHD) dataset, but applied them to the Raw Heidelberg Digits dataset (Cramer
et al., 2022). We preprocessed the dataset similarly to Neftci et al. (2019),
creating Mel spectograms of shape 40 (input channels) by 80 (timesteps). The
only additional parameter we needed to tune was the $L1$ regularization
strength. We first trained our model without any sparsity or Dale’s principle-
derived constraints and achieved $94.5$% test accuracy similar to results
reported on SHD by Hammouamri et al.. We then analysed the learnt spatio-
temporal patterns by extracting “spatio-temporal receptive fields” using the
process illustrated in figure 1. We created a panel for each hidden neuron by
placing a color-coded weight value at the learnt delay value on the x-axis,
yielding panels with rows representing input neurons, where each row had a
single either blue or red pixel. We then summed these figures for each neuron
in the output layer, weighted by the learnt hidden-to-output weight, and
adjusted along the x-axis according to the learnt hidden-to-output delay. The
results for the unconstrained network is shown in figure 2. Visually, the
figure clearly illustrates that spatio-temporal patterns of excitation and
inhibition formed after training.

To get a statistical measure of the emerging spatio-temporal patterns we
measured 2D autocorrelation using the method of Moran’s I (Moran, 1950). For
this technique a spatial weighting needs to be defined, we decided to choose
the 8 neighbourhood case (also known as the Queen’s case) to capture a broader
influence. The average value over all output neurons was 0.0901, with the
highest value being 0.1486 for class 14 and the lowest 0.0480 for class 7.
This can be interpreted as spatio-temporal autocorrelations appearing for all
classes, even though they are not very strong (with the highest possible value
being 1).

Next, we trained networks with the biological constraints – utilizing DEEP R
and RigL to enforce both a fixed level of sparsity and Dale’s principle by
making each neuron’s outgoing connections all have the same sign. Visually, it
seems that the spatio-temporal grouping of excitatory and inhibitory neurons
also emerges in this case, as shown in figure 3. We again used Moran’s I with
an 8-neighbourhood matrix, and found an average Moran’s I value of only
0.0184. While this is greater than zero, it is very close to 0 and some output
neurons even have negative values (with a small magnitude). This suggests that
in the biologically constrained case the emergence of coherent spatio-temporal
patterns is less strong, but this is expected, since these networks have very
few connections within the 8 neighbourhood.

Next, we compared the classification performance of the dense network and the
the sparse Dalean network. Starting with a dense network, we progressively
halved the number of connections in subsequent experiments. In the sparse
setting, this corresponds to dropping half of each neuron’s connections on
average whereas, in the dense setting, this corresponds to halving the number
of neurons. The effectiveness of DEEP R and RigL seemed highly dependent on
initialization, so we ran three training sessions for each configuration and
reported the average classification accuracy. As figure 4A illustrates,
initially, the two networks have similar accuracy but, as the number of
synapses decreases, and the individual neurons are forced to compute more in
the dense network, the sparse Dalean network retains better performance.

Finally, we analysed the combined and separate usefulness of dynamic pruning
(i.e. structure learning) and delay learning. We defined four models, one
where the structure is learnt, one where the synaptic delays are learnt, one
where both are learnt and one where neither are (the synaptic weights are
still trained in all cases). For these experiments, we omitted Dale’s
principle, since for the models with fixed connectivity, poor initialisation
could impose significant constraints on the network. We ran three simulations
for each configuration, reporting the average results. Figure 4B demonstrates
that dynamic pruning is highly effective, especially as sparsity increases.
Learning delays is beneficial when the structure is fixed, but its benefits
are less clear when the structure is learnt.

Overall, dynamic pruning seems vital for maintaining high performance in
sparse networks, but we would argue that, when applying such methods, delay
learning might not be necessary.

##  6 Discussion

The computational benefits of Dale’s principle are not obvious but, many
researchers have experimented with methods which enable Dalean networks to
perform similarly to their non-Dalean counterparts (Cornford et al., 2021). In
our setting, Dale’s principle did not seem to impair classification
performance. This may be because, in very sparse networks where neurons only
make a few connections, making all of those connections to have the same sign
is not a very strong constraint.

While delay learning seems useful with sparse connectivity, we found that
learning the structure is more important. When the structure was learnt, delay
learning did not show a clear advantage. This may be because gradient descent
struggles with simultaneously learning both parameters; if a newly connected
synapse has a non-optimized delay, it might be immediately deactivated again.
However, our experiments with using RigL for removing connections (which does
not happen every epoch like Deep R) did not show significant benefits,
challenging this theory. From another perspective, structure learning with
fixed delays _is_ a form of delay learning since, if a synapse has an
ineffective delay, the network can adapt by introducing a more effective
synapse with a new random delay. This method was one of the first proposals
for delay learning in SNNs (Bohte et al., 2002), and is usually referred to as
“delay select”. However, for this method to work in the fully connected
setting, several connections are required between each pre and postsynaptic
neuron, with the number increasing as the delay distribution widens. Ergo,
from the point of view of efficiency, we have an argument for delay learning
methods such as the one based on DCLS, especially in the fully connected
setting. However, we argue that a good baseline of how useful a given delay
learning method is, is to test if in the sparse setting it performs better
than allowing the network to pick a synapse from a pool with a fixed random
delay.

In conclusion, our results show that biologically plausible sparse network
architectures can be efficient for machine learning tasks. Our findings pave
the way for future research into the optimization of SNNs for various
applications, particularly those at the edge that have strong memory and
computational constraints.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any
commercial or financial relationships that could be construed as a potential
conflict of interest.

## Author Contributions

BM – Conceptualization, Investigation, Formal Analysis, Software, Validation,
Visualization, Writing - original draft, Writing - review & editing
JCK – Conceptualization, Supervision, Funding acquisition, Writing - review &
editing
TN – Conceptualization, Supervision, Funding acquisition, Writing - review &
editing

## Funding

This work was supported through a Doctoral Scholarship by the Leverhulme Trust
and the EPSRC (EP/S030964/1, EP/V052241/1).

## Data Availability Statement

The code will be published on Github and generated example networks on the
University of Sussex figshare repository upon acceptance of the manuscript.

## References

  * Ba and Caruana (2014) Ba, J. and Caruana, R. (2014).  Do deep nets really need to be deep?  _Advances in neural information processing systems_ 27
  * Bellec et al. (2018) [Dataset] Bellec, G., Kappel, D., Maass, W., and Legenstein, R. (2018).  Deep rewiring: Training very sparse deep networks
  * Bellec et al. (2020) Bellec, G., Scherr, F., Subramoney, A., Hajek, E., Salaj, D., Legenstein, R., et al. (2020).  A solution to the learning dilemma for recurrent networks of spiking neurons.  _Nature communications_ 11, 3625
  * Bengtsson et al. (2005) Bengtsson, S. L., Nagy, Z., Skare, S., Forsman, L., Forssberg, H., and Ullén, F. (2005).  Extensive piano practicing has regionally specific effects on white matter development.  _Nature neuroscience_ 8, 1148–1150
  * Bohte et al. (2002) Bohte, S. M., Kok, J. N., and La Poutre, H. (2002).  Error-backpropagation in temporally encoded networks of spiking neurons.  _Neurocomputing_ 48, 17–37
  * Cornford et al. (2021) Cornford, J., Kalajdzievski, D., Leite, M., Lamarquette, A., Kullmann, D. M., and Richards, B. A. (2021).  Learning to live with dale’s principle: {ANN}s with separate excitatory and inhibitory units.  In _International Conference on Learning Representations_
  * Cramer et al. (2022) Cramer, B., Stradmann, Y., Schemmel, J., and Zenke, F. (2022).  The heidelberg spiking data sets for the systematic evaluation of spiking neural networks.  _IEEE Transactions on Neural Networks and Learning Systems_ 33, 2744–2757.  10.1109/tnnls.2020.3044364
  * Davies et al. (2018) Davies, M., Srinivasa, N., Lin, T.-H., Chinya, G., Cao, Y., Choday, S. H., et al. (2018).  Loihi: A neuromorphic manycore processor with on-chip learning.  _Ieee Micro_ 38, 82–99
  * Evci et al. (2021) [Dataset] Evci, U., Gale, T., Menick, J., Castro, P. S., and Elsen, E. (2021).  Rigging the lottery: Making all tickets winners
  * Furber et al. (2014) Furber, S. B., Galluppi, F., Temple, S., and Plana, L. A. (2014).  The spinnaker project.  _Proceedings of the IEEE_ 102, 652–665
  * Hammouamri et al. (2023) Hammouamri, I., Khalfaoui-Hassani, I., and Masquelier, T. (2023).  Learning delays in spiking neural networks using dilated convolutions with learnable spacings.  _arXiv preprint arXiv:2306.17670_
  * Han et al. (2015) Han, S., Pool, J., Tran, J., and Dally, W. (2015).  Learning both weights and connections for efficient neural network.  _Advances in neural information processing systems_ 28
  * Hassibi et al. (1993) Hassibi, B., Stork, D. G., and Wolff, G. J. (1993).  Optimal brain surgeon and general network pruning.  In _IEEE international conference on neural networks_ (IEEE), 293–299
  * Izhikevich (2006) Izhikevich, E. M. (2006).  Polychronization: computation with spikes.  _Neural computation_ 18, 245–282
  * LeCun et al. (1989) LeCun, Y., Denker, J., and Solla, S. (1989).  Optimal brain damage.  _Advances in neural information processing systems_ 2
  * Maass (1997) Maass, W. (1997).  Networks of spiking neurons: the third generation of neural network models.  _Neural networks_ 10, 1659–1671
  * Merolla et al. (2014) Merolla, P. A., Arthur, J. V., Alvarez-Icaza, R., Cassidy, A. S., Sawada, J., Akopyan, F., et al. (2014).  A million spiking-neuron integrated circuit with a scalable communication network and interface.  _Science_ 345, 668–673
  * Moran (1950) Moran, P. A. (1950).  Notes on continuous stochastic phenomena.  _Biometrika_ 37, 17–23
  * Neftci et al. (2019) [Dataset] Neftci, E. O., Mostafa, H., and Zenke, F. (2019).  Surrogate gradient learning in spiking neural networks
  * Renyi (1959) Renyi, E. (1959).  On random graph.  _Publicationes Mathematicate_ 6, 290–297
  * Schemmel et al. (2010) Schemmel, J., Brüderle, D., Grübl, A., Hock, M., Meier, K., and Millner, S. (2010).  A wafer-scale neuromorphic hardware system for large-scale neural modeling.  In _2010 ieee international symposium on circuits and systems (iscas)_ (IEEE), 1947–1950
  * Seidl et al. (2010) Seidl, A. H., Rubel, E. W., and Harris, D. M. (2010).  Mechanisms for adjusting interaural time differences to achieve binaural coincidence detection.  _Journal of Neuroscience_ 30, 70–80
  * Zenke and Vogels (2021) Zenke, F. and Vogels, T. P. (2021).  The remarkable robustness of surrogate gradient learning for instilling complex function in spiking neural networks.  _Neural computation_ 33, 899–925