Multilayer Network Analysis of Nuclear Reactions

The nuclear reaction network is usually studied via precise calculation of differential equation sets, and much research interest has been focused on the characteristics of nuclides, such as half-life and size limit. In this paper, however, we adopt the methods from both multilayer and reaction networks, and obtain a distinctive view by mapping all the nuclear reactions in JINA REACLIB database into a directed network with 4 layers: neutron, proton, 4He and the remainder. The layer names correspond to reaction types decided by the currency particles consumed. This combined approach reveals that, in the remainder layer, the β-stability has high correlation with node degree difference and overlapping coefficient. Moreover, when reaction rates are considered as node strength, we find that, at lower temperatures, nuclide half-life scales reciprocally with its out-strength. The connection between physical properties and topological characteristics may help to explore the boundary of the nuclide chart.

The universe synthesizes lighter particles into heavier ones, and the mechanism which remains a heated area of exploration, includes element nucleosynthesis 1 and nuclear processes [2][3][4][5][6][7][8] . Much attention has been drawn on topics like the limits of nuclear binding and the island of stability [9][10][11] . Those studies help us to understand nuclear reactions taking place both in the distant galaxies and in the reactors on the earth. Detailed studies of nuclide consumption or production through reactions are required with the support of massive nuclear data. The JINA REACLIB database, maintained by the Joint Institute for Nuclear Astrophysics, is a well-known source of thermonuclear reaction rates. It aims to compile a complete set of nuclear reaction rates that is continuously updated and regularly snapshotted 12 . The dataset contains various kinds of nuclear reactions, parameters for calculating reaction rates, Q values, etc., providing rich information for nuclear-related research directions.
Besides acquiring exact nuclide abundance by solving sets of time-dependent differential equations in the database, it is challenging to treat the nuclear reaction system as a complex network to explore its statistical characteristics. The basic idea is introduced from graph theory, which considers the interacting units in a system as nodes and the relationship between two units as an edge, thus the system can be studied as a graph (network). The topic of complex networks has achieved significant advances since the 'small world' 13 and 'scale-free' 14 characteristics were found prevalent in many real world systems such as social connections, the Internet and distributed infrastructures 15 . The structure and dynamics of networks mapped from those systems turn out to be distinct from that of regular or random networks [16][17][18] , and complex networks outperform them in modeling real world systems. For example, the 'scale-free' structure of the Internet, which is hierarchical with many hubs, explains how easy it is for viruses to propagate. These findings help us to understand the systems at more profound levels and have benefited researches in various areas 19 , hopefully including nuclear reactions.
An important branch of the studies on complex networks covered metabolic networks or chemical reaction networks. The topological scaling properties were first analyzed of 43 organisms from different kinds of life in a directed network approach and similarity was found within and outside biological systems 20 . It is believed that 'scale-free' was crucial for the robustness. Moreover, the chemistry of planetary atmospheres and the interstellar medium were compared, and it was demonstrated that only the network of the Earth, coupling with its biosphere, was scale-free 21 . Efforts were made for obtaining more results from different systems like interstellar chemistry [22][23][24] , trying to explain the origins of complexity in large systems or to guide experiments 25 . A majority of the studies map their reaction systems into networks using substrate-product method, which treats reactants and products as nodes and connects them with edges if they are in the same reaction equation. In reference 26 , several mapping methods were evaluated. It was concluded that substrate-product and substance networks could best preserve modularity information. Although effective to a certain degree, as further discussion was made on precision of the methods used 27,28 , the methods for mapping reactions into networks still need improvement to better fit a specific system in order to avoid inaccuracy. Worth mentioning, the bipartite graph is also an efficient way to characterize reactions 29 . It maps the substance and reaction into two kinds of nodes in one network, but the analysis of bipartite graph is complicated. We believe the best way is to extract the most information with a simpler method.
In the present work, we study a directed weighted nuclear reaction network based on substrate-product method, including all nuclides and reactions from the JINA REACLIB database. Additionally, we label each edge with 4 types: neutron, proton, 4 He and the remainder, hence a 4-layered network. The multilayer network provides more precise description of a system, without which one might draw misinforming conclusions 30 . To efficiently analyze the multilayer network, frameworks and models are studied [31][32][33][34][35][36] , and metrics for multilayer structure are proposed such as edge overlap 37,38 , various centralities 39,40 and correlations [41][42][43] , adding a new dimension to complex networks research. Detailed reviews on multilayer network are given in reference 44-46 . According to the dependence between nodes in different layers [47][48][49] , there are generally 3 types of multilayer network: one-to-one, partially dependent and one-to-many. We focus on the one-to-one dependent situation, and the multilayer substrate-product method falls into the multidimensional network category, which is originated from colored edge network. There are by far some attempts on biomedicine in reaction/metabolic systems using multilayer method. Here, the advantage of multilayer network analysis helps us to find some structural characteristics of nuclides relating to their physical property, which is a possible path to address the nuclide stability problem that nuclear physics and astrophysics pay close attention to.

Results
Multilayer reaction network. The dataset used in this paper is the latest version of REACLIBV2.0 containing 8048 nuclides treated as nodes and 82851 reactions mapped into directed edges (pairs from each reactant to each product). Over 3000 nuclides, either exist naturally or are synthesized artificially, whose half-life (t h ) data are complemented from NuDat (http://www.nndc.bnl.gov/nudat2/), while the rest are from theoretical calculation. Take 40 Ca as an example, the recorded reactions involving this nuclide are as follows:  According to the substrate-product method, for the first reaction of 40 Ca, we have 4 edges: one from 40 Ca to n and one from 40 Ca to 39 Ca, and vice versa because of the reverse reaction. However, as neutron, proton and 4 He are the most needed particles in nuclear reactions, more than half of the reactions turn a nuclide into other nuclides by capturing or releasing one or more of these three particles, making their degrees k (number of edges) close to 10 4 . If they are considered as normal nodes like other nuclides, whose degrees are smaller than 10 2 , their role of intermediary material will be lost in the topology and the degree distribution will be highly biased (see Supplementary Figure 1). The nodes with extremely high degrees in reaction systems, often referred to as currency nodes, can make the network difficult to analyze, and they are usually deleted in reaction networks. Here, instead of completely deleting these currency nodes, we label an edge with one of the three particles ('n' for neutron, 'p' for proton or 'h' for 4 He) when a nuclide needs this very particle to produce another nuclide. If none of the three particles are needed as reactant in a reaction, the edges from reactants to products will be labeled as 'the remainder' or 'r' , which is not as simple as the situation for the other 3 layers and some information can be lost. Hence, the first reaction of 40 Ca has only one r-edge from 40 Ca to 39 Ca and one n-edge from 39 Ca to 40 Ca. The corresponding topologies of 40 Ca's 4 layers are shown in Fig. 1.
Note that there are mainly 4 ways to map a reaction into a network: substrate-product network (connecting reactants to products), substrate-substrate network (reactants to reactants and products to products), substance network (the combination of the first and second methods), and reaction-reaction network (taking reactions as nodes and connecting 2 reactions if they share the same substance). Many reactions have multiple substrates while an edge in the network analysis usually connects only one pair of nodes, thus some topological structure could be lost when mapping a system to a network. However, in the case of nuclear reactions, the majority only have one or two reactants (see Supplementary Note 1), so that the reactions are basically about one nuclide to another, making the information lost very limited. Also, if there are multiple products, some products can be decided following the conservation of mass (e.g. the 4 th reaction of 40 Ca, also written as 40 Ca(n,p) 40 K, the information about p is redundant). In this sense, the complete reaction information can be preserved. For other methods, the connection between reactants or products, as in substrate-substrate and substance network, could bring in misleading information, while reaction-reaction network loses too much. So the substrate-product network is our best choice. The topologies of the whole system can be found in Supplementary Figure 2.

Topological characteristics and stable nuclides. The network is usually studied via its adjacency matrix
A. The matrix element a ij = 1 if there is an edge from node i to j, otherwise a ij = 0. In a directed network, in-degree and out-degree are counted separately by summing up A in different directions. As for multilayer network, we have A [α] for layer α, where, in our nuclear reaction network, M = 4 and α ∈ {n, p, h, r}. A is denoted as the aggregated network of those 4 layers, The degree distributions of aggregated network and r-layer on a Z-N plane are shown in Fig. 2. The degrees of neutron-rich nuclides are mostly greater than those of proton-rich ones. While for n-, p-and h-layers, all the degree distributions have a peak at k = 3, the networks are very homogeneous. We also notice that k and k m (averaging degrees of nearest neighbors) are highly correlated, which can be recognized as assortativity and has been found prevalent in many other networks. k [r] and k nn The abundant topological characteristics indicate there is rich information stored in the structure of the aggregated network and r-layer, so we explore them further by subtracting in-degree with out-degree as shown in the lower panels of Fig. 2. We can see that the degree difference of r-layer looks similar to the aggregated network, but there is a clear path-like positioned nodes whose in-degree is larger than out-degree (red dots), which are surrounded by scattered nodes that have smaller in-degree than out-degree (blue dots) along the path. On the other hand, the aggregated network does not possess this kind of pattern. As Fig. 3 shows, which is an enlarged version of panel (f) in Fig. 2, this path fits the beta-stability line smoothly, and we would infer that this pattern describes most the positions of stable nuclides. Also, the remainder of the red part may be related to the r-process of nucleosynthesis.
To determine the characteristics of stable nuclides quantitatively, we compared the degree differences of the aggregated network and r-layer, and found that over 80% of the stable nuclides follow the condition (see the inset of Fig. 3). As in-degree and out-degree correspond to the general ability of production and consumption of a nuclide, the aggregated network describes the total balance, so it is reasonable for k i,in = k i,out to be one of the stability conditions. While in r-layer, where the n, p, 4  [ ] could result that the stable nuclides can be accumulated. As Fig. 4 shows, this condition does not apply to the stable nuclides with number of protons smaller than 15, which is probably because the area has more longer edges according to our measure. Thus we can conclude that for the Z ≥ 15 part, the topology is more regular, while the rest part has more complexity. Furthermore, the stable condition indicates that the degree differences in other 3 layers are not the same for in and out situations, so that the 39      A larger λ indicates higher probability for a reaction to occur, which can be mapped to the notion of edge weight in network theory, describing how easy it is to travel between nodes or how important the edge is. Although k in and k out are highly correlated in binary network (upper and middle panels of Fig. 2), s in and s out have various correlation under different temperature values as the reaction rate is a function of T 9 (see Supplementary Figure 7). This could be the reason why a 4-layered network is able to lay foundation for all  complex systems at atomic level with 3 layers being rather regular. The complexity comes from the diversity of edge weights. Of the temperature range we explored, it is observed that for T 9 ≤ 3, a reciprocal relationship The reason for this phenomenon can be explained that half-life is determined by how easy it is for a nuclide to convert into others, and r-layer contains mainly reactions that are classified as one nuclide to others without any currency particles. Also, it indicates that the life-time of a nuclide does not depend on its input channel by which a nuclide is formed, i.e. the life-time of a certain nucleus has no memory for the formation history of itself.

Discussion
The whole set of nuclear reactions is one of the most fundamental systems in the world. The network constructed is quite different from previously studied ones. The average shortest distance is a bit longer than the network diameter, which indicates that it is not such a small world. The degree distribution is narrow with peaks, far from the power-law often observed. Three layers of the 4-layered network are quite regular. The nodes are spatially constrained, lacking enough long range edges to make the communication or propagation more efficient, yet they can generate much more complicated items like chemical molecules and biological tissues. The temperature-dependent strength may take the credit for the complexity. Despite the special structure, some interesting node characteristics are observed from the complex network point of view.
The multilayer network construction method is effective, even some of the products are deleted, we can still extract important information from the network. As stability is about the balance of in and out corresponding to production and consumption in reactions, the degree difference and node overlapping coefficient are able to identify the majority of stable nuclides, and the relationship between half-life and out-strength in r-layer are discovered. The results are obtained with the help of multidimensional method, linking nuclide characteristics with topological structures.
When we focus on a complicated reaction process (e.g. r-process is about continually capturing neutrons to produce heavier nuclides), it usually occurs only in a certain layer, even if the network structure is simple. Therefore, as we study nucleosynthesis and nuclear processes, based on the primitive results of nuclear reaction network analysis, our future work will focus on a better designed mapping method, the dynamics on the network and more detailed analysis of each layer. Node strength. We calculate the rates of each reaction and assign the values to the edges respectively as edge weight. The rate value ranges from approximately 10 −300 to 10 50 at T 9 = 0.1, and varies as T 9 grows from less than 0.01 up to 10. The weights of the two edges with opposite directions between a pair of nodes are usually not the same as the rates of a reaction and its reverse reaction often differ. For T 9 ≥ 0.01, the calculation of reverse rates should be corrected by including partition functions, as they are calculated via detailed balance without partition functions in the dataset. To understand the characteristics of a node, we derive node strength s i by summing up edge weights from or to all the neighboring nodes as shown in Equation 5.