# On the use of graph theory to interpret the output results from a Monte-Carlo depletion code

Benjamin Dechenaux

► To cite this version:

Benjamin Dechenaux. On the use of graph theory to interpret the output results from a Monte-Carlo depletion code. Journal of Nuclear Engineering and Radiation Science, 2021, 195 (5), pp.538-554. 10.1080/00295639.2020.1847980. hal-03191108

HAL Id: hal-03191108

https://hal.science/hal-03191108v1

Submitted on 6 Apr 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
On the use of graph theory to interpret the output results

from a Monte-Carlo depletion code

Benjamin Dechenaux*,ª

ªInstitut de Radioprotection et de Sûreté Nucléaire
31 avenue de la division Leclerc, 92260 Fontenay-aux-Roses, France

*Email: benjamin.dechenaux@irsn.fr

Number of pages: 31

Number of tables: 1

Number of figures: 17
Abstract

The analysis of the results of a depletion code is often considered a tedious and delicate task for
it requires both the processing of large volumes of information (the time dependent composition
of up to thousands of isotopes) and an extensive knowledge of nuclear reactions and associated
nuclear data. From these observations, dedicated developments have been integrated to the up-
coming version of the Monte-Carlo depletion code VESTA 2.2, in order to implement an innovative
representation of depletion problems. The aim is to provide users with an adaptable and efficient
framework to ease the analysis of the results of the code and facilitate their interpretation. This
effort ultimately culminated in the development of the representation of the isotopic evolution of
a given system as a directed graph.

In this paper, it is shown that the Bateman equation encoded in the VESTA code indeed
possesses a natural interpretation in terms of a directed cyclic graph and it is proposed to explore
some of the insight one can gain from the graph representation of a depletion problem. Starting
from the new capabilities of the code, it is shown how one can build on the wealth of existing
methods of graph theory in order to gain useful information about the nuclear reactions taking
place in a material under irradiation. The graph representation of a depletion problem being
especially simple in activation problems, for then only a limited number of nuclides and reactions
are involved, the graph representation and its associated tools will be used to study the evolution
of the structure materials of a simplified model of the ITER fusion reactor.

**Keywords** — Graph theory, Monte-Carlo, Depletion, VESTA
# I. INTRODUCTION

A graph is undoubtedly the most natural and convenient way to pictorially represent the interactions existing between different objects. As such, it has long been in use in a wide variety of domains, with applications to the modelling of social interactions, road network studies, internet mapping and many more [1]. In the field of nuclear science, a graph can be used to visually represent the web of nuclear reactions occurring between nuclides. It is indeed the natural way nuclear scientists picture nuclear reaction chains, by imagining a set of arrows pointing from one isotope to another, to form an entangled web of nuclear reactions. An example of such system is displayed on Figure 1, which shows the graph associated with the decay chain of $^{222}$Rn, as extracted from the JEFF-3.3 nuclear data library [2]. The graph is built using the NetworkX [3] Python library.

<figure><img src="image_2.png" /><figcaption>Fig. 1. The decay chain of $^{222}$Rn.</figcaption></figure>

Although useful, this type of representation seems to reach its limit when the number of nuclides and reactions involved increases, since the intricacy and number of nuclear reactions and
nuclides to keep track of makes any visual representation too complicated. This situation is often encountered by users of fuel depletion codes, such as VESTA [4], where the number of nuclear reactions and isotopes involved make the processing of the results of the code a complex and delicate matter. This complexity might certainly be the main reason why graph theory was never envisaged as a tool to process the output results from a depletion code, although it has proven to be an efficient computational companion for the numerical integration of the Bateman equations [5], [6]. It is our goal to prove that this effectiveness can also be put to good use by providing users with the possibility to use the tools and representations associated with graph theory for the analysis of the calculation results. Another reason for the limited use of graph theory in this particular context might be the availability of efficient and user friendly tools to manipulate graph objects. The analysis proposed in this paper made extensive use of the NetworkX [3] Python package and the Gephi [7] graph visualisation and analysis software. Relying on robust, state of the art methods, both tools are capable to sustain the complexity inherent to depletion problems. Most of the methods presented in this paper are based on simple algebraic algorithms and can be deployed without having to suffer from long computation times, even for the processing of the largest networks one could encounter in depletion problems.

## II. GRAPH REPRESENTATION OF A DEPLETION PROBLEM

### II.A. Preliminary definitions

In mathematical terms, a graph $G$ is whatever object which is given by the specification of two particular sets :

* a set of *vertices* $V$;

* a set of *edges* $E$, being a subset of $V^2: E \subseteq \{(x, y) | x, y \in V\}$.

It is clear from Figure 1 that edges in the graph will be used to encode for nuclear reactions occurring between nuclides, themselves represented as vertices. But, because nuclear reactions are always occurring in a unique direction (i.e. from a source to a target nuclide), one needs to specify that each pair $(x, y)$ of the set of edges is to be understood as an *ordered* pair, that is to say that the pairs $(x, y)$ and $(y, x)$ must be viewed as two distinct elements. When this additional structure is given, the underlying graph is said to be *directed*.
Note also that in the definition of the set of edges $E$, nothing prevents the appearance of edges of the form $(x, x)$, i.e. "loops" starting and ending at the same vertex. If, on the model of Figure 1, each vertex is to represent a specific isotope of the system, then loops might be used to represent isomeric transitions occurring when an isotope is produced in an excited state. In depletion problems, because isomeric states can have a non-negligible impact on the behaviour of the system, they must however be treated as being completely separate entities than their associated ground states. For this reason, a different vertex must be associated in the graph with each isomeric state of a nuclide. Working at this isomeric state level, there is no notion of vertex possessing self loops. In mathematical jargon, graphs that do not possess self loops are called *simple graphs*.

Finally, numerical values can be assigned to vertices or edges. We then say that the graph is *weighted*. In the practical application of graph theory to depletion problems, it appears natural to assign to each vertex the composition of the isomeric state to which it is associated and to each edge the reaction rate of the associated nuclear reaction. With this additional structure, two different visions of the problem can be envisaged. The first is to consider the problem as flows being exchanged between objects on the fixed, pre-existing grid of all nuclides and possible nuclear reactions. For most purposes, and for visualisation and analysis purposes in particular, a second point of view must be favoured, which is to envisage the graph itself as a dynamical entity. As the materials' compositions evolve in time, new vertices and edges will appear or disappear in the graph.

The main idea underlying this paper is to try to extract as much information as possible about the behaviour of a material under irradiation by conscientiously analysing snapshots of the simple directed graph representation of the system at different moments in time.

## II.B. Monte-Carlo fuel depletion simulation with VESTA 2.2

VESTA [4] is a code whose object is to provide an efficient, accurate and user friendly interface between Monte-Carlo codes used for neutron transport in matter, such as MCNP 6 [8] or MORET 5 [9], and fuel depletion modules such as ORIGEN [10], FISPACT-II [11] or PHOENIX (VESTA built-in depletion module). This structure makes VESTA a robust and efficient tool to
model the evolution of materials in complex geometries and makes the code applicable in a wide
variety of contexts, from the modelling of commercial BWR and PWR fission reactors to research
facilities, fusion devices and many others.

A typical VESTA simulation consists of an iterative process, in which a simulation of the transport of neutrons is first carried out. This Monte-Carlo simulation is used to determine the reaction rates of the neutron induced reactions that occur in materials being simulated. This information is then used in a second step that solves the Bateman equations encoding the evolution of the materials' isotopic inventories.

The results of a typical simulation are the compositions of the burned materials after each
depletion step. These concretely come in the form of lists of concentrations, representing a large
volume of information that, as stated in introduction, is generally difficult to process. The main
reason being that when analysing such data, one is left unequipped to disentangle the intricacy of
the nuclear reactions occurring in between nuclides and so making sense of where do nuclides come
from or disappear is a complicated matter that generally requires knowledge of nuclear reactions
and a solid experience of depletion problems. Note that even if it is possible, in codes such as
VESTA, to ask for information about any specific nuclear reaction occurring in the system, an
individual treatment has to be performed on each of the thousands of nuclear reactions to fully
embrace the complexity of the problem. From these observations, some efforts were put and are
still ongoing in the development of the latest version 2.2 of the VESTA code and its associated
AURORA 1.4 [12] analysis software, in order to provide the user with powerful and systematic
ways to ease the processing of depletion output results.

First, a new output file, in XML format, is provided by the VESTA code. This output file
contains the information about all the reaction rates of all states within the system for each time
step. Used in conjunction with compositions data, this new output file basically contains an image
of the directed graph of the system. This new output file can directly be read by the AURORA
analysis software which can subsequently either transform and export the graph in the "Graph
Modelling Language" data format compatible with most graph theory tools or can be used to
directly launch the Gephi graph visualisation software from within its graphical interface. In a
second time, the possibility to develop and add custom plugins to the Gephi software has been used
to design specific tools and methods to ease the analysis of a depletion problem. Those developed
plugins will be made publicly available under a GPL licence.

As already stated, being able to visualise data in the form of a graph is not in itself a sufficient condition to gain insight about a simulation result. But now all of the machinery behind graph theory can be deployed to our advantage. In the following, it is proposed to explore some of the tools that can be employed to analyse the results from a typical depletion application example.

### III. AN APPLICATION EXAMPLE IN A FUSION ACTIVATION PROBLEM

To illustrate the usefulness of graph theory as a tool to study the output results from a depletion calculation, it is proposed to apply this technique in a concrete application example. To this aim, a very simplified toy model of the ITER fusion device is presented in this section. This model will be used to simulate the activation of materials that compose a section of the torus of the ITER reactor.

The model comprises concentric, cylindrical and homogeneous layers of materials. The cylinders are one-meter long and reflection conditions are applied to each of their ends. The model is therefore made, for neutron transport purposes, into an infinite cylinder. A picture of one fourth of the model is shown on Figure 2. All of the dimensions and compositions in the model are extracted and adapted from [13].

It is assumed in the following that a fusion plasma is continuously sustained for a whole year in the vacuum chamber of the torus, which corresponds to the innermost part of the system. For irradiation purposes, it is supposed that the plasma can be viewed as a constant source of 14 MeV neutrons, emitted along the cylindrical symmetry axis of the system with homogeneous angular dependence. Considering a nominal plasma power of 500 MW carried away exclusively by 14 MeV neutrons, an estimate of $2.10^{20}$ neutrons emitted by the plasma each second is obtained.

In the following, we will focus our attention on the first, innermost material layer of the model (i.e. the blanket layer shown in yellow on figure 2). The precise isotopic inventory of the blanket before irradiation is given on Table I.

A VESTA simulation of the system is performed using the nuclear data library JEFF-3.1 [14]. The simulation is based on MCNP 6.1 for the fixed source neutron transport simulation and the PHOENIX depletion module for simulating the evolution of the composition of materials over time. The composition and (instantaneous) nuclear reaction rates in the blanket are extracted
<figure><img src="image_3.png" /><figcaption>Fig. 2. One fourth view of the geometric model (dimensions not to scale).</figcaption></figure>

<table><caption>TABLE I
Initial composition of the blanket (natural abundances).</caption><thead><tr><th>Isotope</th><th>composition<br>(atom/barn/cm)</th></tr></thead><tbody><tr><td>H</td><td>1.34 10<sup>-2</sup></td></tr><tr><td>B</td><td>4.18 10<sup>-3</sup></td></tr><tr><td>N</td><td>3.03 10<sup>-5</sup></td></tr><tr><td>O</td><td>6.71 10<sup>-3</sup></td></tr><tr><td>Si</td><td>2.11 10<sup>-4</sup></td></tr><tr><td>P</td><td>8.97 10<sup>-6</sup></td></tr><tr><td>Si</td><td>6.35 10<sup>-6</sup></td></tr><tr><td>Cr</td><td>9.33 10<sup>-3</sup></td></tr><tr><td>Mn</td><td>8.90 10<sup>-4</sup></td></tr><tr><td>Fe</td><td>2.86 10<sup>-2</sup></td></tr><tr><td>Co</td><td>2.56 10<sup>-5</sup></td></tr><tr><td>Ni</td><td>7.71 10<sup>-3</sup></td></tr><tr><td>Cu</td><td>3.49 10<sup>-3</sup></td></tr><tr><td>Nb</td><td>3.13 10<sup>-6</sup></td></tr></tbody></table>

at different times and are analysed using the graph representation described in the introduction. Figure 3 shows the graph representing the evolved composition of the blanket after five days of irradiation. The graph is made of 309 vertices, representing all states including isomeric states in
the system and the 1531 nuclear reactions (directed edges) occurring in between them.

<figure><img src="image_4.png" /><figcaption>Fig. 3. Composition of the blanket after 5 days of irradiation.</figcaption></figure>

## IV. PROPERTIES EXTRACTED FROM THE TOPOLOGICAL STRUCTURE OF THE GRAPH

In this section, a focus on the graph of Figure 3 is proposed. The graph represent the state of the system after five days of irradiation. Vertices in red represent isotopes that were initially present in the system while vertices coloured in orange represent isotopes that were created during the neutron irradiation. Edges in blue represent the instantaneous web of nuclear reactions. Vertices are placed on the graph using the Fruchterman-Reingold force directed algorithm. Finally, here and throughout the rest of this paper, quantities suffixed by an underscore (e.g. "Co60_1") are to be understood as isomeric states of an original nuclide.
To analyse the graph of the system, it is proposed to use first the properties that are asso-
ciated with the global "topological" structure of the graph. The goal pursued in this section is to
exploit the basic topological structure of the graph to our advantage, reducing when possible the
dimensionality of the problem.

IV.A. Strongly and weakly connected components

One of the first basic question one could ask when dealing with a graph of any sort as to
do with the number of separate pieces it contains. A graph can indeed be constituted by different
pieces or *connected components* that does not interact with one another. For undirected graphs
(i.e. no notion of direction is defined on the edges), this has a straightforward meaning and the
number of components simply reflects the fact that the object is made up of separate independent
pieces that can be analysed separately. The question is more subtle when working with directed
graphs, for which connectedness can bear two separate meanings, embodied in the notions of *strong*
and *weak* components of a graph. Weakly connected components of a directed graph correspond to
the intuitive notion of separate pieces that the graph is made of. But when the graph is directed,
a situation might arise where regions (i.e. subsets of vertices) can be accessed by following a path
along the edges of the graph but cannot be left because there is no directed way out of the subset
of vertices. In this latter situation, the subset of vertices is said to form a *strongly connected
component* of the graph.

The graph of Figure 3 only has one weakly connected component, i.e. the graph of the
system comes in one piece. But it has 8 strongly connected components, which are highlighted on
Figure 4. In this figure, it is observed that the graph appears to be built mainly around three of
the strongly connected components. All of the other strongly connected components are composed
of "dead ends", i.e. stable nuclides. Stable nuclides indeed satisfy the property that they can be
reached along a path but can never be left. They therefore form strongly connected components
of "unit size".

From figure 4, one can see that the subset comprising of Hydrogen, Helium and Lithium
appears to form a particular substructure and plays an important role in the graph. This family
of isotopes indeed seems to be connected to a lot of other vertices (hence its central position in the
representation) and form a strongly connected component. Those light elements are produced in a
<figure><img src="image_5.png" /><figcaption>Fig. 4. Each nuclide on this figure is colored based on the strongly connected component to which it belongs.</figcaption></figure>

number of different ways but seem to never contribute to the formation of heavier elements. This
can be better seen on Figure 5, where only light elements are shown, all the other being condensed
into a unique "Other" vertex. It is observed that the system formed of Hydrogen, Helium and
Lithium form an intertwined network of connections but never contribute to the production of
elements heavier than $^{8}$Li.

The production of light elements, Tritium in particular, is crucial for a device like ITER
and a further study of the structure of the subgraph formed by the light elements might reveal
number of interesting physical information. For the sake of the current exercise however, a focus
on the production of heavier elements will be favoured. In this context, the strongly connected
component of light elements does not play any more role in the system and the corresponding
vertices and edges might be deleted from the graph. Doing so reveals a drastic simplification in
<figure><img src="image_6.png" /><figcaption>Fig. 5. The structure of the Hydrogen-Helium-Lithium strongly connected component of the graph.</figcaption></figure>

the graph structure, that could have hardly been anticipated: the resulting graph of Figure 6 is
now made of two weakly connected components.

The explanation for this separation into two components is, here also, perfectly clear from
the figure. The upper component of the graph is composed of the ensemble of the daughter nuclides
originating from a single initial isotope, namely $^{93}$Nb, while the second is mainly driven by the
lighter elements of the iron alloy composing the blanket. In five days of irradiation, there has not
been time for the nuclear reactions to produce a sufficient number of elements for the two subgraphs
to merge. Obviously the two components will not stay separated for the whole irradiation time but
if one is only interested in short irradiation period, the separation of the system into two unrelated
components is a drastic simplification of the system.

IV.B. Identifying bridges in the graph

Related to the notion of connected components of a graph is the notion of *bridges* or *cut-edge*.
A bridge is an edge whose removal increase the number of weakly connected components of the
graph. The graph of the ITER model after 5 days of irradiation does contain bridges but only of
the most trivial sort: they consist of edges whose removal simply disconnect a single nuclide from
the rest of the graph. This simply reflects the fact that there is only one reaction in the system
<figure><img src="image_7.png" /><figcaption>Fig. 6. The initial graph after removing light elements up to ⁸Li.</figcaption></figure>

producing this particular nuclide.

As said before, the graph of Figure 6 will not stay disconnected for the whole irradiation period. New nuclides will appear in the system and with them, new reactions will slowly introduce more and more connections between the two components. The concept of bridgesª is then a convenient mean to characterize the appearance of those connections. The state of the system after 15 days of irradiation is presented on Figure 7. In this situation, the (n, α) nuclear reaction ⁷⁷³⁴Se + n → ⁷⁴³²Ge + α is a bridge between the two previously separated components of the graph. Its dynamical behaviour is therefore crucial to understand the interaction between the two pieces that otherwise can be analysed separately.

<footnote>ªand its generalisation to cut-sets of edges, allowing to search for sets of edges whose removal disconnect the graph.</footnote>
<figure><img src="image_8.png" /><figcaption>Fig. 7. Representation of the graph of the system after 15 days of irradiation. New nuclides that appeared in between the two time steps are shown in green.</figcaption></figure>

IV.C. Modularity classes

Returning now to the graph of the system after 5 days of irradiation, it is supposed one is
only interested in nuclides belonging to the lower, "iron driven", component of the graph. This
component contains 196 nuclides and a bit more than 600 associated nuclear reactions, still much
more than what can be easily analysed. Nevertheless, there is another topological property of the
graph that can be exploited to further break down the system into smaller parts. Indeed, the
elongated nature of the component clearly express a physical aspect of nuclear reactions : most
reactions occur in between close neighbour, i.e. nuclides sharing similar values of atomic number
Z. There is actually many ways in graph theory to formally express this particular idea. In this
section, it is implemented using the notion of *modularity classes* [15]. It consists of a heuristic
grouping of nuclides being densely interconnected. Figure 8 shows the modularity classes obtained for the iron alloy component of the graph.

<figure><img src="image_9.png" /><figcaption>Fig. 8. Modularity classes built for the iron component of the graph.</figcaption></figure>

By using modularity classes, one can manage to perform a further subdivision of the problem into smaller bits, which makes the problem much more tractable. The idea is to perform first a separate analysis in each separate modularity class of the system and then take into account the interactions between adjacent classes. By using modularity grouping of the system, one is then assured that the number of interactions in between different groups is reduced to a minimal amount of reactions.

To conclude this section, let us note two important remarks concerning the use of modularity classes for the analysis of a depletion problem. Firstly, it should be kept in mind that the division of the problem in modularity classes is heuristic and does not bear – a priori – any particular physical meaning other than a rough estimate of the interactions between regions in the graph. In particular, the division is not unique and applying the same algorithm twice might produce different splitting into classes.
The second remark concerns the applicability of modularity algorithms in problems involving fission. In our particular example, the use of modularity classes appears to be an efficient way to subdivide the system because of the particular elongated shape the graph presents, which is a characteristic of activation problems. In problems involving fission, this particular topology will certainly be lost, because fission add edges between nuclides with very different atomic numbers. In this case, modularity will certainly not be the best tool to analyse the problem and one needs to resort to more advanced or custom methods tailored for this particular problem. In the next section, an example of subdivision that could also be employed when fission is occurring in the system will be provided.

# V. SIMPLE PATHS AND THE RECONSTRUCTION OF DECAY CHAINS

Taking advantage of the topological properties of the graph underlying the system allowed us to subdivide the problem into more manageable pieces. The question now remains to analyse the dynamic structure of the nuclear reactions that occur in each of these parts. To do so, it is proposed to resort to fast and robust methods that exist to navigate along *paths* in the system, i.e. ordered sequences of vertices connected by directed edges. This is achieved by the notion of simple path, which is a path in a graph that does not contain the same vertex twice (the path does not pass twice at the same place). The ability to efficiently navigate along paths in a graph, using well known, available methods [16], is one of the most convenient tools that graph theory has to offer. It will make it possible to reconstruct in an efficient and robust manner all of the reaction chains occurring in the system.

## V.A. Identifying the initial nuclides for the production of a given nuclide

In the present case, where one is interested in the behaviour of the evolved composition at early times, one question that naturally arises is to identify, for each nuclide produced during the irradiation, the initial isotopes that contributed the most to its production. In other words, one would like to associate to each nuclide a notion of how likely it was produced following a chain of reactions that originates from a reaction of one of the initial isotopes composing the material. In practice, the association is constructed by following all of the simple paths from the initial isotopes $N_i^{T=0}$ to a nuclide of interest $N_i$.
Using the well known algorithms of graph theory, it is straightforward to obtain all of the simple paths from one vertex to another (path finding algorithms achieve $o(V + E)$ time). In the application sought, the key point is to assign a probability to each individual path from $N_i^{T=0}$ to $N$. To do this, one needs to have a good view of the system and adapt the notion of probability so as to match the question asked about the system.

The most simple solution one could think of would be to associate each edge in the path with the branching ratio of the reaction it encodes. For one path from $N_i^{T=0}$ to $N$, we could then multiply the branching ratio of the edges composing the path, as illustrated on Figure 9. The numbers obtained could then be interpreted as the probability to end up in vertex $N$ by following the reaction chain along a path starting from the isotope $N_i^{T=0}$. Summing over all paths from $N_i^{T=0}$ to $N$, one can then obtain the total probability to produce some nuclide $N$ starting from a position $N_i^{T=0}$ in the graph.

<figure><img src="image_10.png" /><figcaption>Fig. 9. Following the paths from a source isotope A to a target isotope D, we can associate a probability that the production of D is due to the isotope A, by carefully multiplying and summing branching ratios.</figcaption></figure>

The method just built is valid but does not really answer the question raised in the beginning of this section, which was to identify, for a given nuclide $N$, the initial isotope it most likely originate from. To answer this question, one has to follow the paths backward and thus associate to each path a notion of "ingoing" branching ratio. This notion can be defined in a simple manner, by the ratio of the rate of the reaction one is interested in (i.e. the reaction on the path under consideration) to the sum of the reactions rates ending on the nuclide. Multiplying the ingoing
branching ratios on each path from $N_i^{T=0}$ to N, one obtains a measure whose interpretation can be thought of as the likelihood of a nuclide N to originate from the initial isotope $N_i^{T=0}$.

Performing the computation for all pairs ($N_i^{T=0}$, N), we can use the probabilities to associate each nuclide in the system with a unique (or at least reduced number of) initial nuclide whose nuclear reactions is most likely to be responsible for its creation. In our fusion activation example, such association was performed for the iron component of the evolved composition after five days of irradiation and is illustrated on Figure 10. Each nuclide in the system is given a unique color which is associated to the initial isotope that most likely was responsible for its creation.

<figure><img src="image_11.png" /><figcaption>Fig. 10. Each nuclide on this figure is attributed a unique colour depending on the specific initial isotope it most likely originate from. (Initial isotopes are coloured in red).</figcaption></figure>

Figure 10 gives an approximate but enlightening idea of what happened in the system during the first few days of irradiation. One can see for instance, that all of the nuclides in blue in the middle of the graph originate from the same initial nuclides (Chromium 50 and 52) and form an intertwined web of rapid successive nuclear reactions connecting ultimately with the less massive elements of the graph. If one is particularly interested in the production of one of the elements in this region, it can be directly understood how it is influenced by the initial composition of
Chromium. Thus, one gets a rough idea of the decay channels that are involved. It can then be
used to further simplify the system and allow for more specific treatments.

Another comment can be made about the resemblance of Figure 10 with the modularity
classes found on Figure 8. It appears quite natural that modularity classes, in identifying groups
of tightly bonded elements, is able to broadly perform an association between nuclides and an
ensemble of initially present isotopes. Thus the colouring of Figure 8 must prefigure the results
obtained by a much finer association, such as the one presented on Figure 10.

To end this section, note that by associating the nuclides that appeared in the system with
isotopes that were initially present, one is effectively building a particular kind of graph – a *bipartite
graph* – where edges always start in one family of vertices (the initial isotopes) and point to all
of the other nuclides. This type of graph can also be an invaluable tool to visually display and
understand the dynamics of the activation problem. An extract of the bipartite graph, constructed
from the association explained in this section, is presented on Figure 11.

<figure><img src="image_12.png" /><figcaption>Fig. 11. An extract of the bipartite graph constructed from the association of one isotope initially present to the isotopes that appeared during the irradiation.</figcaption></figure>

Let us state finally that the pairing that was made in this section must be taken with a
bit of caution and do not aim at being something else than an heuristic manner to analyse in a
simplified fashion the nuclear reaction chains. Some more elaborate methods could be put in place,
in particular to take into account the timescales associated to each nuclear reaction taking place
and thus to gain more insight on the dynamics of the system.

It is nevertheless believed that the method devised in this section gives a good idea of what
is happening in the system. It can be used to further subdivide the graph, allowing for the much
finer analysis that is proposed in the next sections.
V.B. Construction and visualisation of the nuclear reaction chains

Following the subdivision of the problem into smaller, more manageable pieces, the next step of the analysis would consist of analysing specific behaviour of substructures of the graph, emphasising on the small groups of nuclides who are playing a predominant roles for the specific applications sought. We will not delve here in a detailed study of the precise understanding of the origin of the structures that we previously found in the graph, because that would take us too far from the main subject of this paper. We rather take advantage of this section to present the last types of graph that we can build from nuclear reaction chains, tree graphs.

Representing nuclear reaction chains as a tree-like structure is perhaps the most practical way to analyze the structure of the production and decay channels of a specific nuclide. Two examples are shown on Figure 12, presenting the production channels associated to the appearance of $^{60}$Co and the disintegration chains of $^{52}$Cr, an isotope for which it was observed in the last section that it occupies an important place in the iron component of the total graph.

<figure><img src="image_13.png" /><figcaption>Fig. 12. Two examples of tree graphs. (a) The production channels of $^{60}$Co. (b) The disintegration channels of $^{52}$Cr.</figcaption></figure>

Although this kind of graph has much to offer for a precise survey of the structure of the nuclear reaction chains, it is in general not a structure that naturally arise in the entangled web of nuclear reactions. This is because, nuclear reaction chains tends to form cycles, that is to say paths which follow closed loops, passing several times through the same position. For this reason,
a specific methodology must be engineered to bring out this type of structure from the total graph of the system. The most practical method is based on the exploitation of simple paths between nuclides, which, because the paths are simple, naturally do not include any cycle. This notion of simple path directly allows us to extract the desired tree-like structure.

# VI. STUDYING DIFFERENCES IN BETWEEN TWO GRAPHS

Because graphs are well defined mathematical objects, simple algebraic operations like addition and subtraction in between two weighted graphs can be given a definite and precise meaning. This is yet another major advantage of graph theory which comes in handy in situations where one wants to compare two systems sharing common properties or irradiation histories. This can happen for instance when one is simulating the same system using two different sets of nuclear data. It can also be an invaluable tool to analyse how the system is evolving in time.

In this last section, it is proposed to highlight some of the features one can extract using simple vertex-wise and edge-wise operations. Those methods will be applied to the ITER activation model, in both of the above mentioned cases, namely the comparison of nuclear data libraries and the comparison of the system at different time steps.

## VI.A. Identifying differences between nuclear data libraries

The results presented thus far were obtained using the JEFF-3.1 nuclear data library. In this section, a comparison of those results with a simulation of the ITER model based on the ENDF/B-VII.1 nuclear data library [17] will be made. For the activated material after five days of neutron irradiation, one ends up with two graphs that can be compared.

The first thing one could do, is to check if the two graphs contain the same number of nodes and edges, that is to say, verify if the irradiation made using two distinct nuclear data libraries involves the same number of nuclides and nuclear reactions. In the present case, the two graphs involve the same number of nuclides but there is one particular reaction that is missing from ENDF/B-VII.1, it is the $\beta^{-} + \alpha$ reaction of $^{17}$N, decaying into $^{13}$C, see Figure 13.

One can see from Figure 13 that the reaction rate associated with this missing reaction is much lower than any of the other production channels of $^{13}$C. One can also show that this channel is also far from being a dominant reaction in the decay of $^{17}$N. So it is concluded that no major
<figure><img src="image_14.png" /><figcaption>Fig. 13. The production tree of ¹³C, shown for the graph constructed from JEFF3.1 (left) and ENDF/B-VII.1 (right). Total reaction rates (in s⁻¹) are displayed for each reaction.</figcaption></figure>

discrepancy will be observed because of the absence of this particular reaction in ENDF/B-VII.1.

For all of the other reactions that the two nuclear data libraries have in common, we can plot their relative differences to survey the overall behaviour of both data sets. The relative differences in reaction rates, obtained for each nuclear reaction within the material are shown on Figure 14.

Note that each point of the graph represents the total disintegration rate for each of nuclides in the system. This rate is proportional to the isotopic concentration of each nuclide "source" of the associated reaction. This remark being said, it should be understood that the highest discrepancies obtained in the graph do not reflect problems in nuclear data but are mainly obtained for short lived nuclides that are produced in very small quantities in the system. Thus, an overall good agreement between the two nuclear data libraries is observed.

## VI.B. Tracking the temporal evolution of the system

Finally, the temporal evolution of the system during the irradiation period is discussed. From a practical point of view, the temporal evolution is given as a set of graphs, one for each time step of the depletion code. Simple manipulations on the set of graphs can be applied to extract useful
<figure><img src="image_15.png" /><figcaption>Fig. 14. Relative differences (ENDF/B-VII.1 - JEFF 3.1)/ENDFB-VII.1 obtained by comparing the instantaneous reactions rates of the fusion activation problem after five days of irradiation. The abscissa represent the atomic number of the nuclide responsible for the reaction.</figcaption></figure>

information about the temporal behaviour of the system.

In doing so, the effort spent understanding the structure of the system during the first moments of the irradiation history will not go wasted. For now one can use the knowledge accumulated and examine each new structure in light of the basic shapes that are by now familiar. As an illustrative example, Figure 15 show the graph of the system after 195 days of irradiation (like before, the Hydrogen-Helium-Lithium structure was suppressed, to gain visibility). As it was mentioned back in section IV.B, the graph is now made into one, weakly connected component. The isotopes that appeared along the evolution are represented using different colours. On overall, it can be seen that, at least for this simple activation example, the material's composition is still built around the basic outline of the first few days, and is little by little enriched with new isotopes. New isotopes rapidly fill in the gap between the initially disconnected components. Once this is done, the only additional effect of irradiation is the creation of heavy elements, that continues to occur up until the very end of the irradiation period.
<figure><img src="image_16.png" /><figcaption>Fig. 15. The graph of the system (without the Hydrogen-Helium-Lithium component) after 195 days of irradiation. Each isotope is coloured based on the time step at which it first appeared.</figcaption></figure>

To conclude this section, let us stress that by looking only at the graph of Figure 15, but this is also true for any other figure of the paper, one is in reality putting aside some of the most important dynamical aspects of the system. Indeed, the concentrations and reactions rates are constantly evolving in the system and the figure does not seem to be able to capture this aspect of the temporal evolution. To understand the true dynamical behaviour of the system, one must resort to a detailed analysis of the composition of each isotopes, in the same fashion that one would do if one had not had graph theory tools at disposal.

It seems like graph theory is not capable of capturing these aspects of the evolving system. But this is merely due to our emphasising of visual aspects linked with the graph representation of the system. In providing a robust and convenient structure to navigate through, graph represen-
tation also is a most efficient container to use in a detailed analysis of more quantitative aspects
of a depletion problem. Using modern and efficient frameworks, such as the NetworkX Python
package, one is left in a position to grasp much more information of a depletion problem.

VII. GRAPH THEORY FOR DEPLETION INVOLVING FISSILE MATERIALS

The activation problem presented in this paper was carefully selected for its simplicity and
because it allowed to motivate the use of the tools and methods of graph theory in a straightforward
and illustrative manner. One might at this point be tempted to doubt the applicability of those
methods to the more intricate problem of depletion in the presence of fissile isotopes. A detailed
analysis of such type of problem would require a paper on its own. Let us nevertheless conclude
this paper by providing the reader with a glimpse of the graph underlying a realistic depletion
problem involving fission.

The graph presented on Figure 16 shows the composition of a PWR type fuel pin at the end
of irradiation. The graph is made of around 1500 states (including isomeric states) and contains
almost 34,000 nuclear reactions. The strongly connected components of the graph are highlighted
using different colours. As before, the notion of strongly connected components of the graph offer a
broad partition of the total graph into categories. In this particular application, it appears that the
global topology of the graph is layered in successive clustered groups of isotopes. The innermost
region of the graph is built around the fissile isotopes in the system. Around this central part, an
arrangement of successive layers is found, staggered in the successive decays of the fission products.

This staggered structure is best understood by considering the graph in Figure 17, which was
obtained after applying simplifications to the graph associated with the total system. First, all of
the isomeric states have been regrouped with their associated ground state, so that only the latter
remain visible on the figure (the reactions associated with isomeric states have not been neglected
but have been condensed on the associated ground states). Then, all of the reactions accounting
for less than 5% of the branching ratio of each isotopes have been removed from the system.

A clear structure now appears in the graph, which is now divided into two main components.
The first component on the upper right is built of the heaviest, fissile elements of the graph. It also
comprises the first few fission products that further decays into the more stable, lighter elements
which are represented on the lower left component of the graph. It appears quite clearly now that
<figure><img src="image_17.png" /><figcaption>Fig. 16. Graph extracted from a fuel pin at the end of irradiation. Each colour in the graph is associated to a strongly connected component of the graph.</figcaption></figure>

the methods presented throughout the paper could now be applied to the simplified system we
have built. Such methods would allow to further reduce the complexity of the problem and to
answer in a relevant and adapted manner to the questions raised by the use of a fuel depletion
code.
<figure><img src="image_18.png" /><figcaption>Fig. 17. The graph of the system after simplifications.</figcaption></figure>

VIII. CONCLUSION

A broad overview of the most commonly known methods used in graph theory has been achieved and their applicability to the analysis of a depletion problem has been demonstrated. Graph theory indeed is a natural representation to visualise an irradiation process but it is also the right tool to disentangle its associated complex web of nuclear reactions and extract from it sensible and useful information, that would otherwise remain out of reach.
Our attention was focused on some of the most simple and readily available methods that a user could use to understand the global structures underlying an activation problem. The few solutions that have been presented, will figure among the methods that will, in term, be made directly available to users of the AURORA analysis software, for the post-processing of the outputs of the VESTA Monte-Carlo depletion code. To keep the directing line of the paper at the level of a broad overview and to emphasise on the methods that AURORA will natively be provided with, the proposed exploration of graph theory presented in this paper is only partial. Many more tools and methods comprised within graph theory could still be applied and help to further understand the delicate problem of nuclear disintegration chains.

Also, many more tools inspired by graph theory could be used to gain meaningful insights about depletion problems. For instance, one could use metrics or indices (*centrality measures* for example) to identify in a given system the isotopes or nuclear reactions that are the most important, either because they drive the production of important groups of nuclides or because they induce reactions whose uncertainty could cause a major problem for estimating the composition of other isotopes. Over methods and algorithms could also include the notions of quotient graphs, multi-commodity flows, the study of discrete optimal transport on the graph and the emerging field of graph neural networks, that could reveal itself to be of great interest either to extract sensible information from a graph or construct fast and efficient *ersatz* of a depletion problem, to predict nuclear inventories in situations where rapidity takes precedence over accuracy.
ACKNOWLEDGMENTS

The author would like to thank Mariya Brovchenko for helpful discussions and comments.
REFERENCES

[1] M. NEWMAN, Networks, Oxford University Press (2018).

[2] A. PLOMPEN, O. CABELLOS, C. DE SAINT JEAN, M. FLEMING, A. ALGORA, M. ANGELONE, P. ARCHIER, E. BAUGE, O. BERSILLON, R. BLOKHIN, F. CANTARGI, A. CHEBBOUBI, C. DÍEZ, H. DUARTE, E. DUPONT, J. DYRDA, B. ERASMUS, L. FIORITO, U. FISCHER, and G. ŽEROVNIK, "The joint evaluated fission and fusion nuclear data library, JEFF-3.3," *The European Physical Journal A*, **56** (2020); 10.1140/epja/s10050-020-00141-9.

[3] A. A. HAGBERG, D. A. SCHULT, and P. J. SWART, "Exploring Network Structure, Dynamics, and Function using NetworkX," G. VAROQUAUX, T. VAUGHT, and J. MILLMAN (Editors), *Proceedings of the 7th Python in Science Conference*, 11 – 15, Pasadena, CA USA (2008).

[4] W. HAECK and B. DECHENAUX, "VESTA 2.2 User's manual," PSN-EXP/SNC/2017-251, Institut de Radioprotection et de Sûreté Nucléaire - France (2018).

[5] J. W. EASTWOOD and M. J. GUY, "Pathway and uncertainty prediction in FISPACT-II," *Joint International Conference on Supercomputing in Nuclear Applications + Monte Carlo* (2013).

[6] S. KASSELMANN, O. SCHITTHELM, F. TANTILLO, S. SCHOLTHAUS, C. RÖSSEL, and H.-J. ALLELEIN, "Development of a new nuclide generation and depletion code using a topological solver based on graph theory," *Nuclear Engineering and Design*, **306**, 154 (2016).

[7] M. BASTIAN, S. HEYMANN, and M. JACOMY, "Gephi: An Open Source Software for Exploring and Manipulating Networks," *International AAAI Conference on Weblogs and Social Media* (2009).

[8] C. J. WERNER ET AL., "MCNP6.2 Release Notes," LA-UR-20808, Los Alamos National Laboratory (2018).

[9] B. COCHET ET AL., "Capabilities overview of the MORET 5 Monte Carlo code," *Joint International Conference on Supercomputing in Nuclear Applications and Monte Carlo 2013 (SNA + MC 2013)*, La Cité des Sciences et de l'Industrie, Paris, France (2013).
[10] S. B. LUDWIG and J. P. RENIER, "Standard- and Extended-Burnup PWR and BWR Reactor Models for the ORIGEN2 Computer Code," ORNL/TM-11018, Oak Ridge National Laboratory (1989).

[11] T. J. DONOVAN and Y. DANON, "FISPACT-II: An Advanced Simulation System for Activation, Transmutation and Material Modelling," *Nuclear Data Sheets*, **139**, 77 (2017); https://doi.org/10.13182/NSE03-A2332.

[12] B. DECHENAUX, "AURORA User's manual," PSN-EXP/SNC/2018-180, Institut de Radio-protection et de Sûreté Nucléaire - France (2018).

[13] A. ARAÚJO ET AL., "Flux and dose rate evaluation of ITER system using MCNP5," *Brazilian Journal of Physics*, **40**, 58 (2010).

[14] A. KONING ET AL., *The JEFF-3.1 Nuclear Data Library*, Nuclear Energy Agency, Organisation for Economic Co-operation and Development (2006).

[15] M. E. J. NEWMAN, "Modularity and community structure in networks," *Proceedings of the National Academy of Sciences*, **103**, 23, 8577-8582 (2006); 10.1073/pnas.0601602103., URL http://dx.doi.org/10.1073/pnas.0601602103.

[16] R. SEDGEWICK, *Algorithms in C++ Part 5: Graph Algorithms*, Addison-Wesley Professional (2001).

[17] D. BROWN, "Release of the ENDF/B-VII.1 Evaluated Nuclear Data File," (2012); 10.2172/1043398.