Xinwei Chai

Tony Ribeiro

Morgan Magnin

Olivier Roux

Katsumi Inoue

Static Analysis and Stochastic Search for Reachability Problem

Keywords: Model checking, Reachability problem, Asynchronous Binary Automata Network, Local Causality Graph, Heuristics, Answer Set Programming

This paper focuses on a major improvement on the analysis of reachability properties in large-scale dynamical biological models. To tackle such models, where classical model checkers fail due to state space explosion led by exhaustive search. Alternative static analysis approaches have been proposed, but they may also fail in certain cases due to non-exhaustive search. In this paper, we introduce a hybrid approach ASPReach, which combines static analysis and stochastic search to break the limits of both approaches. We tackle this issue on a modeling framework we recently introduced, Asynchronous Binary Automata Network (ABAN). We show that ASPReach is able to analyze efficiently some reachability properties which could not be solved by existing methods. We studied also various cases from biological literature, emphasizing the merits of our approach in terms of conclusiveness and performance.

Introduction

With increasing quantities of available data provided by new technologies, e.g. DNA microarray [START_REF] Marx | Biology: The big challenges of big data[END_REF], there is a growing need for expressive modelings and their related high-performance analytic tools. Among them, works on concurrent systems have been of interest in systems biology for over a decade [START_REF] Bockmayr | Using hybrid concurrent constraint programming to model dynamic biological systems[END_REF][START_REF] Bortolussi | Modeling biological systems in stochastic concurrent constraint programming[END_REF][START_REF] Wiley | Computational modeling of the EGF-receptor system: a paradigm for systems biology[END_REF]. If model validation is a major concern, one of the main challenges nowadays is predicting the behavior of these systems.

Reachability problem on formal models is a critical challenge where both validation problems (whether the model satisfies the a priori knowledge) and prediction problems (properties to be discovered) meet. From a formal point of view, numerous biological properties in computational models can be transformed to reachability properties. For example, the reachability of state 0/1 of a could represent the activation/inhibition of certain gene or synthesis of a protein, while initial state could represent initial observation in an experiment. If the reachability of a certain state contradicts with the a priori knowledge, one can modify the model and/or design a new experiment to verify whether there is an error in the a priori knowledge or former observation. Also, reachability analysis is of help to medicine design: for example if one wants to prevent the carcinogenesis of a cell (target state), one possible solution is to find the critical pathways towards the target state and design a medicine to cut them in order to keep the cell healthy.

In the domain of model checking, reachability has been of great interest for over 30 years [START_REF] Clarke | The birth of model checking[END_REF][START_REF] Edmund | 2 5 years of model checking[END_REF]. Various modeling frameworks and semantics in bioinformatics have been studied: Boolean network [START_REF] Akutsu | Control of boolean networks: Hardness results and algorithms for tree structured networks[END_REF], Petri nets [START_REF] Mayr | An algorithm for the general Petri net reachability problem[END_REF][START_REF] Esparza | Reachability in live and safe free-choice Petri nets is NPcomplete[END_REF], timedautomata [START_REF] Daws | Model checking of real-time reachability properties using abstractions[END_REF][START_REF] Woźna | Checking reachability properties for timed automata via SAT[END_REF]. These approaches rely on global search and thus face state explosion problem as the state space grows exponentially with the number of variables. In [START_REF] James L Peterson | Petri nets[END_REF], it has been shown that the reachability problem of Petri net is exponential time-hard and exponential space-hard, and this conclusion does not change even under some specific conditions [START_REF] Esparza | Reachability in live and safe free-choice Petri nets is NPcomplete[END_REF]. For 1-safe Petri nets, the complexity of reachability analysis is generally PSPACE-complete [START_REF] Cheng | Complexity results for 1-safe nets[END_REF]. Li et al. [START_REF] Li | On reachability and controllability of switched boolean control networks[END_REF][START_REF] Li | Stability analysis for switched boolean networks under arbitrary switching signals[END_REF] investigated theoretically the stability, the controllability and the reachability of Switched Boolean Networks, but their method remains computationally expensive; Saadatpour et al. [START_REF] Saadatpour | Attractor analysis of asynchronous boolean models of signal transduction networks[END_REF] researched only the reachability of fixed points.

To tackle the complexity issue, symbolic model checking [START_REF] Burch | Symbolic model checking: 10 20 states and beyond[END_REF] based on ordered binary decision diagrams (OBDDs) and SAT-solvers (satisfiability) [START_REF] Parosh | Symbolic reachability analysis based on sat-solvers[END_REF] have been studied over years, but still fail to analyze big biological systems with more than 1000 variables. Bounded Model Checking (BMC) [START_REF] Clarke | Bounded model checking using satisfiability solving[END_REF] is an efficient approach but generally not complete as its searching depth is limited to a given integer k.

Beside these approaches, abstraction is an efficient strategy to deal with such models of big scale. It aims at approximating the model while keeping the most important parts influencing the reachability. Abstract methods often have better time-memory performance but with a loss of information. They solve usually a simplified version of the original model, i.e. the results from these approaches are not necessarily compatible with all the properties of the original model. While studying reachability problems, the system dynamics is abstracted to static causalities between states and transitions.

We have designed a new discrete modeling framework for a concurrent system [START_REF] Chai | A Heuristic for Reachability Problem in Asynchronous Binary Automata Networks[END_REF]: Asynchronous Binary Automata Network (ABAN). In ABAN, we applied the approach developed by Paulevé et al. [START_REF] Paulevé | Reduction of qualitative models of biological networks for transient dynamics analysis[END_REF][START_REF] Folschette | Sufficient conditions for reachability in automata networks with priorities[END_REF][START_REF] Paulevé | Refining dynamics of gene regulatory networks in a stochastic π-calculus framework[END_REF] to address reachability problem. This approach refers to a static abstraction of the reachability (with an over-approximation and an under-approximation of the real dynamics). It is based on an abstract interpretation: Local Causality Graph (LCG). This interpretation drastically reduces the searching state-space thus avoids costly global search [START_REF] Paulevé | Static analysis of biological regulatory networks dynamics using abstract interpretation[END_REF]. However, this pure static analysis is not complete as there are inconclusive cases which can not be decided reachable or not.

Many biological networks are encoded in Boolean style, e.g. [START_REF] Akutsu | Control of boolean networks: Hardness results and algorithms for tree structured networks[END_REF][START_REF] Kauffman | Homeostasis and differentiation in random genetic control networks[END_REF], because BN is a simple formalism but with strong applicability: discretization in BN is a way to handle the imprecision of a priori knowledge on the model. However BN may be not expressive enough. As to the modeling of the dynamic behavior "a ← 1 at moment t + 1 if b = 1 at moment t", one has a(t + 1) = b(t) in BN. a always follows the evolution of b but with a potentially unwanted behavior "a ← 0 when b = 0 at moment t". ABAN models this dynamics as via {b 1 } → a 1 without this redundancy. Besides, BNs are transformable to ABANs, and this property makes our approach applicable to a wider domain (Appendix A).

Our work shares similar concerns but we combine static analysis and bounded model checking. We have developed a heuristic approach PermReach to attack reachability problem [START_REF] Chai | A Heuristic for Reachability Problem in Asynchronous Binary Automata Networks[END_REF] which is more conclusive than pure static analysis but time-consuming and still not able to solve reachability problems under certain conditions. In this paper, we propose a hybrid approach ASPReach based on the former LCG reasoning and a non-exhaustive search in the LCG to obtain a more conclusive solution of reachability problems. ASPReach allows one to solve the cases where other static methods fail. Furthermore, it can also solve the reachability of a set of states which to our knowledge has never been done in static way. We assess the value of our contribution using benchmarks on biological examples from the literature. This paper is organized as follows: Section 2 introduces the formal background and the formalization necessary to the understanding of the work; Section 3 presents the concrete methods and algorithms composing the whole approach; Section 4 shows the benchmarks evaluating our approach and other alternatives; Section 5 concludes this paper.

Formalization

Notations:

a i means automaton a is taking value i;

x :: y is the sequential connector of entities x and y, where x appears just before y; a.next is the immediate successor of a; a.pred is the immediate predecessor of a.

Boolean Network is a traditional and efficient modeling framework, with many biological networks encoded in BN [START_REF] Kauffman | Homeostasis and differentiation in random genetic control networks[END_REF]. To describe the dynamical properties more precisely, Automata Network is introduced [START_REF] Chai | A Heuristic for Reachability Problem in Asynchronous Binary Automata Networks[END_REF][START_REF] Plateau | Stochastic automata network of modeling parallel systems[END_REF]. It can be considered as a subset of communicating finite state machines or safe Petri Nets.

Asynchronous Binary Automata Network (ABAN) is a special case of Automata Network. "Asynchronous" implies the update scheme that no more than one automaton can change its value at a time. Asynchronous update scheme makes a trade-off between the biological reality and complexity. From a given state, biological systems may evolve to multiple future states (e.g. cell differentiation), but it is costly to simulate generalized ANs where arbitrary subsets of automata can be updated simultaneously. "Binary" implies that every automaton has exactly two possible states (0, 1). In Appendix A, we show how to transform a Boolean Network into ABAN.

Definition 1 (ABAN). An ABAN is a triplet AB = (Σ, L, T ), where: To describe the evolution in an ABAN, we use the notion of trajectory. Definition 3 (Trajectory). Given an ABAN AB = (Σ, L, T ) and a global initial state α ∈ L, a trajectory t from α is a sequence of transitions t = tr 1 :: As to reachability problem, given an ABAN AB = (Σ, L, T ), the joint reachability REACH(α, Ω) is formalized as: joint state Ω is reachable iff there exists a trajectory t s.t. α • t = Ω. Partial reachability reach(α, ω) is defined analogously: local state ω = a i is reachable iff there exists a trajectory t s.t. (α • t)[a] = a i . REACH(α, Ω) and reach(α, ω) take Boolean values True, False or Inconclusive if it cannot be decided. In Fig. 1, Ω = a 1 , b 1 , c 1 , d 1 , e 0 or ω = a 1 is reachable from initial state α via trajectory t, i.e. reach(a 1 ) = True and REACH(α, Ω) = True. In fact, the reachability of a joint state even a global state is equivalent to that of one local state. Given a set of target states {a i , . . . z j }, by adding a new automaton x to Σ, setting its initial state to x 0 and adding transition {a i , . . . z j } → x 1 to T , REACH(α, {a i , . . . z j }) is thus equivalent to reach(α, x 1 ). For convenience, we study partial reachabilities in this paper.

• Σ {a, b, . . .}
Paulevé et al. [START_REF] Paulevé | Refining dynamics of gene regulatory networks in a stochastic π-calculus framework[END_REF] have proposed Local Causality Graph (LCG) to analyze reachability problems statically. LCG abstracts the original problem through an over-approximation (necessary condition) and an under-approximation (sufficient condition). It is a very efficient tool as there is no global search and all the operations are bounded in polynomial complexity. However LCG does not guarantee to obtain a result, i.e. some inconclusive instances satisfy the necessary condition and fail sufficient conditions. In this paper, we make use of LCG by removing some elements needed only in multivalued networks, then we try to analyze it more deeply to solve inconclusive cases of binary valued systems. In fact, Didier et al. have shown a technique to transform a multivalued network to Boolean network [START_REF] Didier | Mapping multivalued onto boolean dynamics[END_REF], which provides us the applicability to multivalued situations.

Definition 5 (LCG). Given an ABAN AB

= (Σ, L, T ), an initial state α and a target state ω, LCG l = (V state , V solution , E) is the smallest recursive structure with E ⊆ (V state × V solution ) ∩ (V solution × V state ) which satisfies: ω ∈ V state a i ∈ V state ⇔ {(a i , sol ai )|a i ∈ α} ⊆ E sol ai ∈ V solution ⇔ {(sol ai , V a (sol ai )} ⊆ E
where V state ⊆ LS is the set of local states, V solution ⊆ T is the set of solutions and V a is the set of required local states of sol ai .

Example 2. Fig. 2 shows the LCG for analyzing reach(a 1 ) in Fig. 1.

a 1 e 1 AND OR b 1 d 0 ∅ c 1 d 1
Figure 2: Visualization of LCG, with the squares representing local states and small circles representing solution nodes. ∅ signifies that there is no need to link any transitions, i.e the former state d 0 is in the initial state.

Algorithm 1 describes how to construct an LCG from an ABAN AB = (Σ, L, T ). Starting from a given target state Ls = ω, one can find all the transitions T s ⊆ T reaching ω and add edges ω → T s . Then we find all the heads A of T s and add edges T s → A and replace Ls with A (recursion). Finally, we update the structure until Ls ⊆ α or there is no transition with body in Ls.

Algorithm 1 Construction of LCG

Input: an ABAN AB = (Σ, L, T ), an initial state α, a target state ω

Output: LCG l = (V state , V solution , E) Initialization: Ls ← {ω}, V state ← ∅, V solution ← ∅, E ← ∅ while Ls = ∅ do Ls = Ls\V state for a i ∈ Ls do Ls ← Ls\{a i } if a i ∈ α then E ← E ∪ {(a i , ∅)} else
// Choose the transitions reaching a i , i.e. with body a 1-i

for sol = A → a 1-i ∈ T do V solution ← V solution ∪ {sol} E ← E ∪ {(a i , sol)} V state ← V state ∪ A for b j ∈ A do E ← E ∪ {(sol, b j )} Ls ← Ls ∪ A V state ← V state ∪ Ls V solution ← V solution ∪ a i .next return (V state , V solution , E)
Intuitively, when the recursive construction is complete, SLCG is in fact a digraph with state nodes V state and solution nodes V solution . E consists of the edges between local state nodes and solution nodes. To access certain local states, at least one of its successor solutions (corresponding transitions from solution nodes) need to be fired; to make one solution node firable, all of its successor local states need to be satisfied. A recursive reasoning of reachability begins with a state node representing target local state, goes through a i → sol ai → b j • • • and ends with initial state (possibly reachable) or a local state without solution successor (unreachable).

With LCG, it is easy to verify whether their are potential pathways from the target state ω to the initial state α. If there does not exist such a pathway, one can ensure that ω is not reachable from α. [START_REF] Paulevé | Static analysis of biological regulatory networks dynamics using abstract interpretation[END_REF] explains this local reasoning.

Definition 6 (Pseudo-reachability). Given an LCG l = (V state , V solution , E) with global initial state α, the pseudo-reachability of node v ∈ V state is defined as

reach (α, v) =      True if v ∈ α False if v ∈ α and ∃(s, sol) ∈ E (s,sol)∈E ( (sol,s)∈E reach (α, s)) otherwise
However, pseudo-reachability is named "pseudo" because it is only an overapproximation of reachability, i.e. it reveals a necessary condition of reachability.

Example 3. In Fig. 3, reach (l, c 1 ) = reach (l, a 1 )∧reach (l, b 1 ) = reach (l, a 0 )∧ reach (l, b 0 ) = True.
Both a 1 and b 1 are reachable, but they can not be reached simultaneously. In such LCG, there are two branches, a 1 → b 0 and b 1 → a 0 , the automata a and b involve themselves in different branches, the reachability of a 1 impedes the reachability of b 1 and vice versa.

a 0 1 b 0 1 c 0 1 {a1, b1} {b0} {a0} c 1 a 1 b 0 ∅ b 1 a 0 ∅ AND Figure 3: Σ = {a, b, c}, T = {{b 0 } → a 1 , {a 0 } → b 1 , {a 1 , b 1 } → c 1 }, ω = c 1
Also, the recursive reasoning does not terminate if there exists cycles in LCG. While computing the pseudo-reachability, self-dependent form reach (l, a i ) = . . . = reach (l, a i ) will appear. Dealing with cycles becomes inevitable. Definition 7 (Cycle). In an LCG, a cycle is formed by a sequence of nodes linked as follows:

a i → • → • • • → • → a i
To identify cycles, we search instead Strongly Connected Components (SCC) of size greater than one. Because cycles may interlace and there is no such problem in SCC. In other words, a SCC may contain several nested cycles which connect to each other. [START_REF] Tarjan | Depth-first search and linear graph algorithms[END_REF] shows that the detection of SCCs can be done in O(|V | + |E|) time, with |V | the cardinality of the vertices and |E| the cardinality of the edges. LCG is usually a sparse graph, as in biological systems, the components mostly interact with only a part of the system, hence the outdegree can be considered of O(1) and the detection of SCCs 1 can be done in O(|V |), i.e. linear time.

Theorem 1. Given a cycle x → • → • • • → • →
x in an LCG, if there is at most one incoming edge to the cycle, the cycle can be removed.

Proof. If there is no incoming edge, the target state y must be in the cycle. The edge y.pred → • → y can be removed, because the reachability of y.pred requires y, but y is the target state, which is never reached before the other local states in the LCG are reached. Thus the transition corresponding to this edge is never fired and the edge can be removed. Similarly, if there is an outside incoming edge a → • → x, a must be the successor of target state y or the target itself, x.pred → • → x can hence be removed. To reach x, we need to reach z, but z cannot depend on x as x is already to be reached. Self-dependence appears: x is reachable if x is reachable. Thus edge z → • → x is deleted (dashed line).

Unfortunately, not all cycles are removable via Theorem 1. Example 5 explains the issue.

Proof. Suppose an arbitrary cycle

C = a i → • • • b j → • • • → a i , with → an edge in the LCG. Note that reach (α, a i ) =⇒ reach (α, b j ) =⇒ reach (α, b j .next) =⇒ • • • =⇒ reach (α,

ASPReach Algorithm

In this section, we present the main contribution of this paper, our analyzer ASPReach: an algorithm for checking the reachability of a target local state ω from a global initial state α (which can also be partial) in a given ABAN. However, exhaustive search leads to heavy computation and huge need of memory. The algorithm proposed below tries to overcome those shortcomings by combining static analysis and stochastic search into the following hybrid approach.

ASPReach:

• Input: An ABAN AB, an initial state α, a target state ω and a max number of iterations k -Back to step (iv)

• Generate all trajectory that starts with α in l using ASP -If a trajectory t ending with ω is found, return True

return Inconclusive

LCG illustrates the causality between necessary transitions to be fired to reach the target state; the tentative of removing cycles simplifies the LCG and keep the reachability unchanged; pseudo-reachability allows one to filter some unreachable cases based on the topology of LCG. The heuristic approach is the core of our algorithm. Stochastic choices avoid combinatorial explosion on different OR Gates. The ASP part searches thoroughly the result but does not traverse the whole state space (ASP solver starts from constraints, finds one consistent order and terminate the search).

Algorithm 2 ASPReach

1: Input: LCG l = (V state , V solution , E), an integer k 2: Output: reachability r and a trajectory t 3: Compute SCCs, classify them into SCC1(l) with at most 1 incoming edge and SCC2(l) otherwise 4: // 1) Break all cycles and prune useless branches Algorithm 2 provides the detailed pseudocode of the algorithm taking an LCG l as input whose detailed construction is given in Algorithm 1. Lines 4-8 delete all cycles with at most one incoming edge. After removing cycles, the LCG may contain nodes without successor. Such nodes can be pruned since they do not lead to initial state (Line 9-18). This preprocessing reduces the search space of the stochastic search performed in step 4. Now l is pruned and might be cycle-free. Static analysis of l can then be used as heuristics to check pseudoreachability (Definition 6) in order to detect some unreachability cases (Lines 19-21) which may conclude before searching. LCG shows the dependencies between local states and transitions. A pathway in LCG suggests a possible trajectory of reaching the target state. If pseudoReach(l, ω) = False, we can ensure that ω is unreachable, as pseudo-reachability checks a necessary condition of reachability. If reach (l, ω) = True, static analysis is not sufficient for reachability analysis. When static analysis fails, a stochastic search is performed at most k times (line 22-29) to find a state sequence from the initial state α to target state ω. If there remain cycles with multiple incoming edges, according to Theorem 2, ω is unreachable. The value of k will be discussed later in the evaluation section. Random choices are made to fix a value for each OR gate of the LCG allowing to perform a reachability check by generating all possible variable assignment order using ASP. Keep in mind that every state node is an OR gate, we have to choose one of its successor solution nodes to access the state. A set of OR gate choices is called an assignment.

5: for each (V state ⊆ V state , V solution ⊆ V solution ) ∈ SCC1(l) do 6: for each v ∈ V state do 7: if ∃(v, v ) ∈ E, v ∈ (V solution \ V solution ) then 8: E ← E\{(v, v )|v ∈ V solution , (v, v ) ∈ E}
for v ∈ V state do 14: if ∃(v, v ) ∈ E then 15: V state ← V state \{v}; E ← E\{(v , v) ∈ E} 16: E ← E\{(v , sol) ∈ E|sol ∈ {sol = (A → a) ∈ V solution |v ∈ A}} 17: V solution ← V solution \{sol = (A → a) ∈ V solution |v ∈ A}
l = (V state , V solution , E ) ← (V state , V solution , E) 25: for v ∈ V state do // Treat each OR gates 26: pick a random element (v, v ) ∈ E 27: E ← E \{(v, v ) ∈ E |v = v } with i ∈ SCC2(l) and i ∈ E

Stochastic Search

As every OR gate has multiple choices, to avoid combinatorial explosion, we use a simple heuristic: choose randomly one assignment for each trial. Then we can construct a new LCG without OR gate, every state node has exactly one successor solution node.

After applying heuristics to delete OR gates, we use ASP (Answer Set Programming) [START_REF] Baral | Knowledge representation, reasoning and declarative problem solving[END_REF] to analyze the newly obtained LCG with only AND gates. ASP is a prolog-like declarative programming paradigm. It uses description and constraints of the problem (called rule) instead of imperative orders. ASP solver tackles problems by generating all the possibilities respecting the constraints. We use Clingo [START_REF] Gebser | Theory solving made easy with clingo 5[END_REF] which is a combination of grounder Gringo and solver Clasp. Given an input program with first-order variables, grounder computes an equivalent ground (variable-free) program for an ASP program, while solver selects admissible solutions (answer sets) in the ground.

A rule is in the following form:

a 0 ← a 1 , . . . , a m , not a m+1 , . . . , not a n .
where the element on the left of the arrow is called head and the ones on the right called body. a 0 is True if a 1 , . . . , a m are True and a m+1 , . . . , a n are False. Some special rules are noteworthy. A rule where m = n = 0 is called a fact and is useful to represent data because the left-hand atom a 0 is thus always True. It is often written without the central arrow. On the other hand, a rule where n > 0 and a 0 =⊥ is called a constraint. As ⊥ can never become True, if the right-hand side of a constraint is True, this invalidates the whole solution.

Constraints are thus useful to filter out unwanted solutions. The symbol ⊥ is usually omitted in a constraint. Programs can yield no answer set, one answer set, or several answer sets. For example, the program b:-not c. c:-not b. produces two answer sets: {b} and {c}. Indeed, the absence of c makes b true, and conversely absence of b makes c true. Cardinality constraints are another way to obtain multiple answer sets. The most usual way of using a cardinality is in place of a 0 : l{q 1 , . . . , q k }u ← a 1 , . . . , a m , not a m+1 , . . . , not a n .

where k ≥ 0, l is an integer and u is an integer or ∞. Such cardinality means that under the condition that the body is satisfied, the answer set X must contain at least l and at most u atoms from the set {q 1 , . . . , q m }, or, in other words: l ≤ |{q 1 , . . . , q m } ∩ X| ≤ u.

ASP Encoding

After deleting OR gates, to encode the reachability problem in ASP, we first describe the facts: Predicate init(a,i) shows the automaton a is at initial state i. Predicate node(a,i,n) shows the node a i in the LCG is numbered n, while parent(n1,n2) expresses node No.n 1 is the predecessor of No.n 2 . The LCG in Fig. 3 is encoded as follows: init(a,0). init(b,0). init(c,0). node(a,1,1). node(b,1,2). node(c,1,3). node(b,0,4). node(c,0,5).

parent [START_REF] Parosh | Symbolic reachability analysis based on sat-solvers[END_REF][START_REF] Akutsu | Control of boolean networks: Hardness results and algorithms for tree structured networks[END_REF]. parent [START_REF] Parosh | Symbolic reachability analysis based on sat-solvers[END_REF][START_REF] Baral | Knowledge representation, reasoning and declarative problem solving[END_REF]. parent [START_REF] Akutsu | Control of boolean networks: Hardness results and algorithms for tree structured networks[END_REF][START_REF] Bortolussi | Modeling biological systems in stochastic concurrent constraint programming[END_REF]. parent [START_REF] Baral | Knowledge representation, reasoning and declarative problem solving[END_REF][START_REF] Bockmayr | Using hybrid concurrent constraint programming to model dynamic biological systems[END_REF].

After the facts, we want the nodes to appear in an order by which we can fire all the transitions sequentially from initial state to target state.

The rough idea is: If different states of one automaton a appear, e.g. a 0 and a 1 . One of them must be in initial state (suppose a 0 ). The transitions with head a 0 have to be fired before a 0 flipping to a 1 , otherwise there is no solution node in the LCG which allows a 1 return to a 0 . In other words, the predecessor of a 0 must appear before a 1 . Core rule describes this constraint.

Predicate prior(N1,N2) signifies node N 1 appears earlier than N 2 in the resulting state sequence. seq(O,a,i) shows that state node a i appears in the O-th place in a trajectory. reachable/unreachable is the final result of the program.

%Rule 1, a node appears always earlier than its predecessor prior(N1,N2) :-parent(N2,N1). %Rule 2, transitivity prior(N1,N3) :-prior(N1,N2), prior(N2,N3). %Rule 3, Core rule prior(N1,N2) :-node(P1,S1,N1), node(P2,S2,N2), node(P2,S3,N3), parent(N1,N3), init(P2,S3), S2!=S3, P1!=P2. %target is unreachable if there is a conflict in order unreachable :-prior(N1,N2), prior(N2,N1), N1<N2. %One node appears once and at least once in a sequence 1seq(1..O,P,S)1 :-O=node(P1,S1,N1):node(P1,S1,N1), node(P,S,N), not unreachable. %Nodes in the sequence are consistent with the order :-prior(N1,N2), node(P1,S1,N1), node(P2,S2,N2), seq(O1,P1,S1), seq(O2,P2,S2), O1>O2. %One place in the sequence cannot be taken by multiple nodes :-seq(O1,P1,S1), seq(O2,P2,S2), P1!=P2, O1=O2. :-seq(O1,P1,S1), seq(O2,P2,S2), S1!=S2, O1=O2. %--------output formatting, displaying initial states first :-seq(O1,P1,S1), seq(O2,P2,S2), init(P1,S1), not init(P2,S2), O1>O2. :-seq(O1,P1,S1), seq(O2,P2,S2), init(P1,S1), init(P2,S2), P1<P2, O1>O2. reachable :-not unreachable.

Notation: a b means a appears before b. When analyzing the LCG in Fig. 3, Rule 1 gives b 0 a 1 , a 1 c 1 , a 0 b 1 , b 1 c 1 ; Rule 2 gives a 0 c 1 and b 0 c 1 ; Rule 3 gives a 1 b 1 and b 1 a 1 which is impossible, therefore there does not exist a state sequence to reach c 1 from initial state. c 1 is unreachable.

If we find a state sequence consistent with all the order constraints, we can obtain its corresponding trajectory, thus we are sure that the target state is reachable. Still, ASPReach is not complete. A counter-example is shown in Fig. 6, when there are multiple branches of one OR gate leading to unreachability, the result can be inconclusive. There is a tricky way to deal with this issue when |ORgates| is not big: we set a limit n, if |ORgates| < n, we shift the heuristics on the assignment of OR gates to the enumeration of all possible assignments. This "hacking" can deal with some inconclusiveness. In the benchmarks in the next section, inconclusive instances appear neither in biological examples nor in random generated tests.

We also show some algorithmic properties of ASPReach:

Theorem 3 (ASPReach termination and correctness). Let l = (V state , V solution , E) be an LCG with initial state α and target local state ω and k > 0 be an integer.

The call ASPReach(l,k) terminates. ASP Reach(l, k)=(False, ∅) only if t a trajectory in l from α to ω. ASP Reach(l, k)=(True, t) only if ∃t a trajectory in l from α to ω. The proof is given in appendix B.

Theorem 4 (ASPReach complexity). Let l = (V state , V solution , E) be an LCG with initial state α and k > 0 be an integer. Let s = |V solution | be the number of target state of l. Let v = |V state | be the number of vertices of l. Let e = |E| be the number of edges of l. The complexity of

ASPReach(l,k) is O(v + e + v/2 × v × e × s + v 2 × e + v × e + k × (v × e 2 + 2 v )) which is bound by O(k × 2 v ).
Proof is given in Appendix B.

Evaluation

In this section we evaluate our approach through various experiments. All tests were run on a Intel Core i7-3770 CPU, 3.4GHz with 8Gb of RAM computer. To validate our approach, we first tested a small model, λ-phage model [START_REF] Thieffry | Dynamical behaviour of biological regulatory networksii. immunity control in bacteriophage lambda[END_REF] to compare with an alternative reachability analyzer Pint [START_REF] Paulevé | Static analysis of biological regulatory networks dynamics using abstract interpretation[END_REF] implementing solely analysis using LCG [START_REF] Paulevé | Reduction of qualitative models of biological networks for transient dynamics analysis[END_REF][START_REF] Folschette | Sufficient conditions for reachability in automata networks with priorities[END_REF][START_REF] Paulevé | Refining dynamics of gene regulatory networks in a stochastic π-calculus framework[END_REF]. In this model with 4 automata and 12 transitions (without taking consideration of the self-regulations), our result shows complete conclusiveness while Pint cannot (Fig. 7). PermReach [START_REF] Chai | A Heuristic for Reachability Problem in Asynchronous Binary Automata Networks[END_REF] is not able either to handle some of the special cases where multiple states of one automaton appear in different branches (Fig. 8). Theses cases are solvable by ASPReach. The whole approach is implemented in Python32 . The call of ASP is done by package pyasp3 .

To evaluate the scalability in in silico networks, we take T-cell Receptor model (TCR) [START_REF] Saez-Rodriguez | A logical model provides insights into T cell receptor signaling[END_REF] and epidermal growth factor receptor model (EGFR) [START_REF] Samaga | The logic of EGFR/ErbB signaling: theoretical properties and analysis of high-throughput data[END_REF] as examples, with the former one containing 95 automata and 206 transitions and the latter one containing 104 automata and 389 transitions respectively.

These models are originally Boolean networks. According to the approach in Appendix A, BNs are transformed into ABANs. Here, we ran the same test as in [START_REF] Folschette | Sufficient conditions for reachability in automata networks with priorities[END_REF]. In the TCR model, we take 3 automata as input (cd4 cd28 tcrlig), varying exhaustively their initial states combinations (2 3 ), take the reachability of states of 5 automata (sre ap1 nfkb nfat sigmab) as output. Similarly we carried a bigger test on EGFR model with 13 automata as input and 12 automata as output. We first tested the performance of traditional model checkers, Mole4 and NuSMV5 , in which Mole turns out to be memory-out for 6 in 12 outputs, and all memory-out for NuSMV in model EGFR. Due to the big state space, traditional model checkers are not applicable. In the TCR tests, our approach gives exactly the same result as Pint did. As for EGFR tests, ASPReach returned no inconclusive output. As seen in Table 1, our approach can be more conclusive than Pint for ABANs. In the configuration of heuristics, we set a threshold for OR gates. If there are less than 10 OR gates after preprocessing, the computation will be shifted from heuristic to the enumeration of all combinations of OR gates. Here is the case for these three benchmarks. The experiments show the ability of ASPReach is already more conclusive than Pint in "simple" cases.

To test the global applicability of ASPReach, we ran two sets of tests on random models generated as follows: Given the number of transitions, for every transition tr, the head a h is randomly chosen from LS, the first element of the body A 1 is chosen from LS 1 = LS \{a h , a 1-h }. For i > 1, if A i-1 = b x exists, we generate A i with an 80% probability, choosing from

LS i = LS i-1 \ {b x , b 1-x }.
Fig. 9 ABAN. In this experiment, we fixed the number of transitions to |Σ| × 3 and each transition has a random number of heads from 1 to |Σ -1|, with Σ the number of automata of the generated ABAN. For traditional model checkers like NuSMV and Mole, memory-out cases begin to appear when |Σ| > 50, so the runtime results are not displayable in Fig. 9(a). Also, Pint is not able to check all test sets of any |Σ|. Even though the curve of average runtime of our approach shows an exponential-like tendency, ASPReach is very fast for |Σ| < 100 (less than 0.03s) and can also perform on models with 1000 automata. Moreover, the longest runtime among the test sets is less than 20s. Because we stop the computation if one reachability check takes more than 20s and we note it as timeout. This result suggests that the heuristic approach of ASPReach has a good performance in conclusiveness in general as the testing models are generated randomly.

Another test is on the different density with the same number of automata. In Fig. 9(b), we fixed |Σ| = 20 and vary the number of the transitions per automaton (density) from 1 to 12. The runtime peak is at density 8, a possible explanation is that even the topology of the network is more complex with the growth of density, more available transitions lead to more pathways from the initial state to the target state, thus the heuristics may end with less trials.

These experimental evaluations show that ASPReach has a better scalability on reachabilty analysis than traditional exhaustive model checkers and a better performance regarding conclusiveness than existing static analyzer.

Conclusion

In this paper, we present the ABAN modeling framework and its model-related reachability analyzer ASPReach. Facing the two critical challenges: complexity and conclusiveness, we combine static analysis and Answer Set Programming for a good performance on both criteria. ASPReach performs normally on models We are now considering an incorporation of heuristics in the preprocessing phase, like the work of [START_REF] Kolčák | Parameter Space Abstraction and Unfolding Semantics of Discrete Regulatory Networks[END_REF], to improve the performance of our approach regarding runtime. The development of dedicated heuristic for the orientation of the random search in ASPReach remains also a future work.

As to biological application, one path of research is the interfacing of our approach with model inference approaches. Biological models can be inferred from experimental data through machine learning techniques (e.g. [START_REF] Inoue | Learning from interpretation transition[END_REF]), but obtained models are not necessarily consistent with the a priori background knowledge about the dynamics of system. Thus it is important to check these obtained models and revise them to make them consistent with such background knowledge. In the case of reachability properties, our approach could be use as a model checker and a mean to enforce such properties in a model under construction.

The call ASPReach(l,k) terminates. ASP Reach(l, k)=(False, ∅) only if t a trajectory in l from α to ω. ASP Reach(l, k)=(True, t) only if ∃t a trajectory in l from α to ω.

Proof. 1: The algorithm starts by breaking all cycles of the LCG and according to Theorem 1 it terminates and does not affect the reachability of α in l.

2: Then all nodes of V state and (resp. solution) with no (resp. missing) outgoing edges are removed. Such nodes cannot be part of a trajectory leading to initial state α and thus this operation does not affect the reachability of α in l. The internal for loop of this operation iterates over V state which is finite. To continue looping, it requires one state deletion thus this operation will terminate atleast when V state becomes ∅.

Conclusion 1: ASP Reach(l, k) = False only if t a trajectory in l from α to v ∈ V solution .

Conclusion 2: the call ASPReach(l,k) terminates. Conclusion 3: After this pre-processing, pseudo reachability is checked and according to [START_REF] Paulevé | Static analysis of biological regulatory networks dynamics using abstract interpretation[END_REF], it terminates and is correct. It is the only possibility for ASPReach to output False.

Conclusion 4: Stochastic search follows by randomly reducing each OR gate of l to one of its edges to form l . This operation is run a finite time k and iterates over V state which is finite and thus it terminates. This operation does not create new edges, i.e. E ⊆ E. ASP solve(l ) generates all possible trajectories of l leading to α. The number of possible trajectory is finite and thus ASP solve(l ) terminates.

Furthermore when ASP solve(l ) = (True, t), t is a trajectory of l proving reachability of α in l and it is the only possibility for ASPReach to output True. 2: To remove useless nodes ASPReach iterates over all states and checking if one state has no successor in l requires to iterates over all edges. In worst case all states will be removed one by one and thus the complexity of this operation is op2 = O(v × (v + s) × e).

3: Computing pseudo reachability over l which have no loop correspond to perform a depth first search on all branch of a tree and thus bound to op3 = O((v + s) × e).

Example 1 .Figure 1 :

 11 Figure 1: An example of ABAN

Figure 4 :

 4 Figure 4: LCG l containing cycle x → • → y → • → z → • → x

Figure 5 :

 5 Figure 5: x, y, z all have external links, thus none of the links can be discarded

  a i ). According to the definition of reach , reach (α, a) = True only if ∃c k ∈ C and c k ∈ α. If there exists such c k , C should not exist as the reasoning stops at c k and does not form a cycle, contradiction. reach (α, a i ) = reach (α, b j ) = • • • = False.

• 3 .•

 3 Output: reach(ω) ∈ {False, True, Inconclusive} 1. Construct the LCG l = LCG(AB, α, ω) 2. Try to remove all cycles and prune useless edges from l Try to prove unreachability of ω in l using pseudo-reachability reach (l, ω) and return False if reach (l, ω) = False 4. Try at most k times • l ← l Simplify each OR gate such that l is a LCG with only AND gates • If there remain cycles:

9 :

 9 // 2) remove useless nodes/edges 10: pruned = True 11: while pruned do 12: pruned = False 13:

  True 19: // 3) Check pseudo-reachability 20: if pseudoReach(l) = False then 21: return (False, ∅) 22: // 4) main search loop 23: for each i in 1 . . . k do 24:

  if r = True then return (True, t) 30: return (Inconclusive, ∅)

Figure 6 :

 6 Figure 6: If an LCG contains such structure, the result could be inconclusive. However the inconclusiveness requires a 1 does not possess other reachable branches.

Figure 7 :Figure 8 :

 78 Figure7: One LCG in λ-phage model, automaton cl appears in both branches of the AND gate. Static analyzer Pint cannot decide whether in this case cro 1 is reachable or not, because it does not consider the order in the state sequence even though there exists a solution is of length 3: cll 0 :: cl 0 :: cro 1 corresponding to the trajectory {cl 1 } → cll 0 :: {cll 0 } → cl 0 :: {cll 0 , cl 0 } → cro 1 .

  .t. ABAN density (b) ABAN size fixed to 20

Figure 9 :

 9 Figure 9: Runtime of ASPReach on random ABANs

Conclusion 3 :

 3 ASP Reach(l, k) = (True, t) only if ∃t a trajectory in l from α to v ∈ V solution .Theorem 6 (ASPReach complexity). Let l = (V state , V solution , E) be an LCG with initial state α and k > 0 be an integer. Let s = |V solution | be the number of target state of l. Let v = |V state | be the number of vertices of l. Let e = |E| be the number of edges of l. The complexity of ASP Reach(l, k) is O(v + s + e + (v + s)/2 × v × e × s + v 2 × e + v × e + k × (v × e 2 + 2 v )) which is bound by O(k × 2 v ).Proof. 1: computing SCC(l) has a complexity of O(v + s + e). In worst case |SCC(l)| = (v + s)/2 and breaking one cycle of SCC(l) is O(v × e × s), thus complexity of removing cycle is op1 = O(v + e + s + (v + s)/2 × v × e × s)

4 :

 4 the stochastic search iterates atmost k times. Treating each OR gate to form l have a cost of O(v × e × e) ASP solve(l ) generates trajectories that can prove reachability of α in l . Each trajectory is a sequence where each element of |V state | appears exactly once. It correspond to the number of total order of |V state | which is 2 v . Thus ASP solve(l ) is bound by O(2 v ) and the whole stochastic search by op4 = O(k × (v × e 2 + 2 v )). Conclusion 1: The complexity of ASP Reach(l, k) is O(op1 + op2 + op3 + op4) = O(v +e+s+(v +s)/2×v ×e×s+v ×(v +s)×e+v ×e+k ×(v ×e 2 +2 v )). Conclusion 2: The complexity of ASP Reach(l, k) is bounded by O(k × 2 v ) C Pseudo Code of Pseudo-reachability

  L} is the set of transitions, where A (called head) is the set of required state(s) for transition tr = A → b i , which allows to flip b 1-i to b i (called body). In other words, transition tr is said fireable iff A ⊆ s, where s is the current global state. A local state represents the state of one automaton, e.g. a 1 , while a global state represents the joint state of all the automata in the network, e.g. a 0 , b 1 , c 0 where L = {a, b, c}.

		is the finite set of automata with every component having
	a Boolean state;
	• LS	∪ a∈Σ {a 0 , a 1 } is the set of all local states, L	× a∈Σ	{a 0 , a 1 } is the set
	of joint states where Σ ⊆ Σ. Particularly, if Σ = Σ, L is the set of global
	states.
	• T {A → b Definition 2 (Asynchronous Dynamics). From current global state s, the global
	state after firing transition tr = A → b i is denoted s•tr = s\{b 1-i }∪{b i }, where
	b 1-i ∈ s. The state of a certain automaton a is noted (s • tr)[a]. Particularly,
	if there is no fireable transition, the next state remains the same as the current
	state.

i | b ∈ Σ ∧ A ∈

  • • • :: tr i :: • • • :: tr n with tr i ∈ T and each tr i is fireable in (s • tr 1 • . . . • tr i-1 ). From α, the global state after firing all transitions of t is (s • tr 1 • . . . • tr n ), denoted s • t. Definition 4 (State sequence). Given an ABAN AB = (Σ, L, T ) and a global initial state α ∈ L and trajectory t, the state sequence seq = s 1 :: • • • :: s i :: • • • :: s n with s i ∈ LS is formed by the updated local states during the trajectory t.

  Example 5. Cycle x → • → y → • → z → • → x is unbreakable according to Theorem 1, is possesses 3 incoming edges. But with Theorem 2, if OR gates are removed, the cycle can be dealt with. The removal of OR gates is stated in the next section.Theorem 2. Given a cycle, if it contains no OR gate, all the local states in the cycles are unreachable.

	1 Python3	implementation	at	https://github.com/alviano/python/blob/master/
	rewrite_aggregates/scc.py

Table 1 :

 1 (a) shows the average run time of ASPReach on randomly generated AR=ASPReach. Results of the tests on small (λ-phage) and large (TCR,EGFR) examples from biology literature. Results of model-checkers using global search are memory-out so are not listed in the table. "Reachable", "Inconclusive" and "Unreachable" give respectively the number of different results of reachability. It is worth noticing that the inconclusive cases in Pint are caused by time-out[START_REF] Folschette | Sufficient conditions for reachability in automata networks with priorities[END_REF].

	Model	λ-phage	TCR	EGFR
	Inputs	4		3		13
	Outputs	4		5		12
	Total tests	2 4 × 4 = 64	2 3 × 5 = 40	2 13 × 12 = 98, 304
	Analyzer	Pint	AR	Pint	AR	Pint	AR
	Reachable	36(56%) 38(59%)	16(40%)	64,282(65.4%)	74,268(75.5%)
	Inconclusive	2(3%)	0(0%)	0(0%)	9,986(10.1%)	0(0%)
	Unreachable	26(41%)	24(60%)	24,036(24.5%)	24,036(24.5%)
	Total time	< 1s	< 1s	7s	40s	9h50min	3h46min

  Algorithm 3 Pseudo-reachability reach Input: an LCG l = (V state , V solution , E), an initial state α, a target state ω Output: a Boolean reach procedure pseudoReach(s) // If ω is in initial state, it is already reached if ω ∈ α then return True // The reachability of s depends on its successor solution nodes if ∃(s, sol) ∈ E then return False for each (s, sol) ∈ E do if fireable(sol) then return True return False procedure fireable(sol) for each (sol, s ) ∈ E do if pseudoReach(s ) then return False return True

Code and testing data available at https://github.com/XinweiChai/LCG-in-ASP

https://pypi.python.org/pypi/pyasp

http://www.lsv.fr/ ~schwoon/tools/mole

http://nusmv.fbk.eu

A Transformation from BNs to ABANs

Given Boolean functions v i (t + 1) = f i (V i ), with V i the set of participating variables among v 1 (t), • • • , v n (t). Boolean functions could be transformed to equivalent CNF (conjunctive normal form) and DNF (disjunctive normal form) if the length of Boolean functions is limited to O(1) [START_REF] Bro Miltersen | On converting CNF to DNF[END_REF] which is often the case.

Proposition 1 (Transformation from BN to ABAN). Given a BN G

B Proofs

Theorem 5 (ASPReach termination and correctness). Let l = (V state , V solution , E) be an LCG with initial state α and target local state ω and k > 0 be an integer.


Arnaud Guertin

Charlotte Duchemin

Arnaud Fardin

Corentin Guigot

Etienne Nigron

Charlotte Remy

Ferid Haddad

Nathalie Michel

Vincent Métivier

Férid Haddad

How nuclear data collected for medical radionuclides production could constrain nuclear codes

à 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.

Introduction

Our research activities are focused on radionuclide production mainly for medical applications, either for therapy or diagnosis. This work is carried out in close collaboration with the GIP ARRONAX that possesses a high energy and high intensity multi-particle cyclotron [1].

In this frame, production cross sections and thick target yields were measured for alpha emitters, such as the U-230/Th-226, Th-227/Ra-223 and Ac-225/ Bi-213 generators [2]; for photon, Tc-99m [2], and positon, Sc-44g [3], emitters for diagnosis; for electron emitters Re-186g [4], Tb-155 [5] and Sn-117m [6] for therapeutic applications. From the irradiated materials, new experimental production cross section data of interest for medical applications and monitor reactions have been extracted which allow to expand our knowledge on these excitation functions, to confirm the existing trends and to give additional values on a wider energy range. These experiments were conducted at the ARRONAX facility using the stacked-foil technique over the four last years.

This data set only represents a small part of the data needed in the field studying the production of innovative radionuclides for medical applications. In order to get answers quickly without the need of new experiments, it is interesting to use theoretical models. The TALYS code gather together several of such models for each step of a nuclear interaction.

A systematic comparison of our results with the output of the TENDL-2015 nuclear data library and TALYS code has been done. In this latter case, several combinations of a e-mail: Arnaud.Guertin@subatech.in2p3.fr models have been tested in order to better reproduce the available data.

Material and methods

Experimental set-up

The production cross section data are obtained using the stacked-foil method [2,7], which consists of the irradiation of a set of thin foils, grouped as patterns. Each pattern contains a target to produce the isotopes of interest. Each target is followed by a monitor foil to have information on the beam intensity thanks to the use of a reference reaction recommended by the International Atomic Energy Agency [8]. A degrader foil is placed after each monitor foil to change the incident beam energy from one target foil to the next one. Each foil in the stack is weighed before irradiation using an accurate scale (±10 -5 g) and scanned to precisely determine its area. The thickness is deduced from these values, assuming that it is homogeneous over the whole surface. In this work, we used thin target and monitor foils of a few tens of micrometres thick and degrader foils of few hundreds of micrometres thick. All foils were purchased from the Good Fellow c company.

These foils were irradiated by proton (up to 70 MeV), deuteron (up to 34 MeV) or alpha particle (68 MeV) beams provided by the ARRONAX cyclotron. It delivers these beams with an energy uncertainty of ±0.50 MeV, ±0.25 MeV and ± 0.61 MeV, respectively for p, d or alpha particle, as specified by the cyclotron provider using simulations. The beam line is under vacuum and closed using a 75 µm thick kapton foil. The stacks were located about 6.8 cm downstream in air.

ND2016

The energy through each foil has been determined in the middle of the thickness of the foil using the SRIM software [9]. Energy losses in the kapton foil and air have been taken into account in our analysis. All along the stack, depending on the number of foils, the energy uncertainty increases up to ±2.0 MeV due to the energy straggling. Irradiations were usually carried out for half an hour, with beams of mean intensity between 100 and 150 nA for proton, between 50 and 140 nA for deuteron and between 140 and 200 nA particles for alpha beam.

The recommended production cross section values [7] of the Ti-nat(d, x)V-48 (all energies), Cu-nat(p,x) Co-56, Zn-62 (> 50 MeV), Ti-nat(p,x)V-48 (< 20 MeV), Ni-nat(p,x)Ni-57 (20-50 MeV), Cu-nat(α,x)Ga-67 (up to 50 MeV) and Al-27(α,x)Na-22 (from 50 MeV to 70 MeV) reactions were used to get information on the beam intensity, depending on the investigated energy range.

The activity measurements in each foil were performed using a high purity germanium detector from Canberra (France) with low-background lead and copper shielding. The first measurements started the day after the irradiation (after a minimum of 15 hours cooling time) during one hour, for all target and monitor foils. Our data are then limited to γ emitter radionuclides with a half-life higher than few hours. A second series of measurements was performed one week after End Of Beam, during a minimum of 24 hours (one day) and up to 60 hours. Third measurements were devoted to long half-life radionuclides and also waiting for the decay of some radionuclides.

Gamma spectra were recorded in a suitable geometry calibrated in energy and efficiency with standard Co-57, Co-60 and Eu-152 gamma sources from LEA-CERCA (France). The full widths at half maxima were 1.04 keV at 122 keV (Co-57 γ ray) and 1.97 keV at 1332 keV (Co-60 γ ray). The samples were placed at a distance of 19 cm from the detector, which is suited to reduce the dead time and the effect of sum peaks. The dead time during the counting was always kept below 10%.

Production cross section calculation

The production cross section values are calculated using the well-known activation formula, defined as a relative equation in which the knowledge of the beam current is no longer necessary thanks to the recommended reactions. The uncertainty is expressed as a propagation error calculation (see [6] for more details).

The TALYS code and TENDL data library

In this work, all the experimental production cross section values are compared with the version 1.6 of the TALYS code released in December, 2013 [10]. TALYS is a nuclear reaction program which simulates reactions induced by light particles on nuclei ranging from Li to Ds. It incorporates theoretical models to predict observables including production cross section values as a function of the incident particle energy (up to 1 GeV). A combination of models that best describes the whole set of available data for all projectiles, targets and incident energies have been defined by the authors and put as default in the code. In this way, a calculation can be performed with minimum information in the input file: the type of projectile and its incident energy, the target type and its mass. The results of this combination of models are referenced in Figs. 1 to 4 as TALYS 1.6 Default.

TENDL is a TALYS-based Evaluated Nuclear Data Library using both default and adjusted TALYS calculations and data from other sources [10,11]. Our experimental production cross section values have been compared with the 8 th TENDL release: TENDL-2015. It provides evaluated data for seven types of incident particles (n, p, d, t, He-3, alpha-particle, gamma) and for all the isotopes living more than 1 second (∼2800 isotopes), up to 200 MeV.

Since there are some differences between experimental data and the results of the TALYS code using default models and the TENDL-2015 library, we have defined a combination of models, already included in the TALYS code, which better describes the production cross sections, for a variety of projectiles, incident energies and target masses.

The description of the optical, pre-equilibrium and level density models have been found to have a great influence on the calculated production cross section values. When proton, deuteron and alpha particle are used as projectile, better results are in general obtained using the optical model described, respectively by [12], [13] and [14]. And for these three projectiles, when a preequilibrium model based on the exciton model including numerical transition rates with optical model for collision probabilities [15] and a model for the microscopic level density from Hilaire's combinatorial tables [16] is used. The results of this combination of models are plotted in Figs. 1 to 4 as TALYS 1.6 Adj. Parametrizations of the models have not been changed. They are used as they are implemented in TALYS. Level density or optical model parameters haven't been tuned.

Results and discussion

TALYS version 1.6 results have been obtained using the default models (labelled Default) and using the combination of models (labelled Adj.) described in the section 2.3. In Figs. 1 to 4, the full circles correspond to our experimental values, the other geometrical symbols to literature data, the full lines to TALYS Default and the dash lines to TALYS Adj. calculations. Four reactions, extracted from previous published articles, are presented. They cover proton, deuteron, alpha particle beams and a wide range of masses.

The Th-226 production from the Th-232(p,3n)Pa-230 reaction

Th-226 is a promising radionuclide for α RIT in leukaemia treatment. Th-226, with a half-life of 31 minutes, can be produced from the α decay of U-230, which has a half-life of 21 days, thanks to a generator system U-230/Th-226. We studied the U-230 production using the Th-232(p,3n)Pa-230(β -)U-230 reaction.

Our set of data for the Th-232(p,3n)Pa-230 reaction is plotted in Fig. 1. There is a good agreement with our experimental data and those available in the literature.

The TALYS code version 1.6 with default models and the TENDL-2015 library values, plotted with a dashedpointed line, are not able to reproduce the production cross section amplitude even if the shape is in good agreement with our experimental data. The TALYS 1.6 Adj. results are in perfect agreement with these experimental data.

The Re-186g production from natural tungsten

Re-186g, with a half-life of 3.7 days, is a β -emitter used in clinical trials for the palliation of painful bone metastases resulting from prostate and breast cancer. It emits a 137 keV gamma ray, suited for SPECT imaging.

The Fig. 2 shows that our experimental production cross section values are in agreement with the literature. The TALYS 1.6 Default calculation is far from the cross section amplitude. The TALYS 1.6 Adj. values are in better agreement but still underestimate the maximum of the cross section by 20%.

The Sn-117m production from natural cadmium

Sn-117m is a conversion and Auger electron emitter, with a half-life of 13.6 days, useful for the palliation of painful bone metastases. Moreover, imaging can be achieved with its 158 keV gamma line.

Our experimental production cross section values are in agreement with the values of Hermanne et al., published in 2010 (see Fig. 3).

However, the Qaim and Döhler values show a shift in energy, roughly 10 MeV higher, for the location of the cross section's maximum. The points have been degraded rather far from the initial beam energy, which could explain the shift by energy straggling through the foils. The TALYS 1.6 Default models give the good shape but results are slightly underestimated and shifted in energy (3 MeV lower). TALYS 1.6 Adj. models are in good agreement with the general behaviour of the experimental data. This combination allows reproducing the shape, the position of the maximum and the amplitude below 45 MeV.

The Sc-47 production from natural titanium

Sc-47, with a half-life of 3.35 days, is a β -emitter of interest for a theranostic approach with Sc-44.

Our experimental production cross section values are plotted in Fig. 4. They are in agreement with the literature data sets even if the amplitude of the first maximum around 40 MeV is lower.

The TALYS 1.6 Default calculation reproduces the right shape but with an overestimation of the amplitude below 50 MeV. TALYS 1.6 Adj. calculation is in good agreement with the general behaviour of the experimental data of A. Hermanne et al. reproducing well the shape and the amplitude.

Conclusions

A large set of production cross section of medical radionuclides have been collected using the stacked-foil technique with proton, deuteron and alpha beams delivered by the ARRONAX cyclotron and for various materials on a wide range of masses. Some comparisons have been systematically performed with the TALYS code, version 1.6. This code has been chosen because it includes a large number of theoretical models, the possibility to combine these models to better describe the experimental data and because the authors' team is very reactive and helpful. Three main mechanisms have been investigated, which have a great impact on the studied observables in our work. They are the optical potential, the level density description and the pre-equilibrium model. A set of models have been found allowing a good description of all our collected data, which is different from the suggested default combination of models. They can be used to get high quality data when no data are available in databases. Further investigations on other mechanisms that can affect our observable are ongoing.

Figure 1 .

 1 Figure 1. Th-232(p,3n)Pa-230 cross section.

Figure 2 .

 2 Figure 2. W-nat(d,x)Re-186g cross section.

Figure 3 .

 3 Figure 3. Cd-nat(α,x)Sn-117m cross section.

Figure 4 .

 4 Figure 4. Ti-nat(α,x)Sc-47 cross section.

c The Authors, published by EDP Sciences. This is an Open Access article distributed under the terms of the Creative Commons Attribution License 4.0 (http://creativecommons.org/licenses/by/4.0/).

The ARRONAX cyclotron is a project promoted by the Regional Council of Pays de la Loire financed by local authorities, the French government and the European Union. This work has been, in part, supported by a grant from the French National Agency for Research called "Investissementsd Avenir", Equipex Arronax-Plus No. ANR-11-EQPX-0004 and Labex No. ANR-11-LABX-0018-01.


Vincent Métivier

Charlotte Duchemin

Arnaud Guertin

Nathalie Michel

Ferid Haddad

Férid Haddad

Thorium-232 fission induced by light charged particles up to 70 MeV

come   L'archive ouverte pluridisciplinaire

Introduction

The irradiation of thorium by light charged particles like protons and deuterons leads to the production of several radionuclides among which radium-223 [1], bismuth-213 [2] and thorium-226 [3], alpha emitters having a great potential in oncologic therapy. The GIP ARRONAX is focused on the production of medically relevant radionuclides and possesses a multi-particle cyclotron [4]. This accelerator has been used for our study on the Th-232(p,x) and Th-232(d,x) reactions. The main motivation was to study the production cross section of Pa-230 that decays to Th-226 [5] via U-230. After the irradiations, the activity values were determined by gamma spectrometry, and the associated spectra give information on the production of several fission products inside the thorium-232 target. The activity of each detectable and quantifiable fission product (FP) has been determined and their associated production cross section has been extracted. From these data, we determined the mass distribution of the FP and the sum of the fission product cross section values. A systematic comparison with the results of the TALYS code has been done.

Material and methods

Experimental set-up

The cross section data are obtained using the stacked-foils method [5,6], which consists of the irradiation of a set of thin foils, grouped as patterns. Each pattern contains a target to produce the isotopes of interest. Each target is followed by a monitor foil to have information on the beam intensity thanks to the use of a reference reaction recommended by the International Atomic Energy Agency [7]. In our experiment, the monitor foil acts also as a catcher to stop the recoil nuclei produced in the target foil. A degrader foil is placed after each monitor foil to a e-mail: Vincent.Metivier@subatech.in2p3.fr change the incident beam energy from one target foil to the next one. Each foil in the stack has been weighed before irradiation using an accurate scale (± 10 -5 g) and scanned to precisely determine its area. The thickness is deduced from these measurements, assuming that it is homogeneous over the whole surface. In this work, we used 10 and 40 µm thick thorium foils, 10 to 25 µm thick titanium, copper and nickel monitor foils depending on the incident particle and its energy, and 100 to 1000 µm thick aluminium or copper degrader foils. These foils were irradiated by the proton (up to 70 MeV) and deuteron (up to 33 MeV) beams provided by the ARRONAX cyclotron. Proton and deuteron beams have, respectively, an energy uncertainty of ± 0.50 MeV and ± 0.25 MeV, as specified by the cyclotron provider using simulations.

The beam line is under vacuum and closed using a 75 µm thick kapton foil. The stacks were located about 6.8 cm downstream in air. The energy through each target and monitor foils has been determined in the middle of the thickness of the foil using the SRIM software [8]. Energy losses in the kapton foil and air have been taken into account. Five stacks were irradiated with protons and five with deuterons, covering respectively, the energy range from 70 MeV down to 11 MeV and from 33 MeV down to 8 MeV. The use of several stacks allows us to minimize the energy uncertainty in our experiments. All along the stack, depending on the number of foils, the energy uncertainty increases up to ± 1.8 MeV due to the energy straggling. Irradiations were carried out for half an hour, with a mean intensity between 100 and 150 nA particles for proton beams and between 50 and 140 nA for deuteron beams. The recommended cross section values [7] of the Ti-nat(d,x)V-48 (all energies), Cu-nat(p,x)Co-56,Zn-62 (> 50 MeV), Ti-nat(p,x)V-48 (< 20 MeV) and Ni-nat(p,x)Ni-57 (20-50 MeV) reactions were used to get information on the beam intensity.

The activity measurements in each foil were performed using a high purity germanium detector with low-background lead and copper shielding. Gamma spectra were recorded in a suitable geometry calibrated in energy and efficiency with standard Co-57, Co-60 and Eu-152 gamma sources. The full widths at half maxima were 1.04 keV at 122 keV (Co-57 γ ray) and 1.97 keV at 1332 keV (Co-60 γ ray). The samples were placed at a distance of 19 cm from the detector which is suitable to reduce the dead time and the effect of sum peaks. The dead time during the counting was always kept below 10%. The first measurements started the day after the irradiation (after a minimum of 15 hours cooling time) during one hour, for all target and monitor/catcher foils. A second series of measurements was performed one week after the end of irradiation, during a minimum of 24 hours (one day) and up to 60 hours. A third series of measurements was done for long half-life radionuclides and waiting for the decay of some radionuclides. Our data are then limited to γ emitter radionuclides with a half-life higher than few hours. FP have been detected and quantified from Zn-72 to Pm-151. The majority of them have been measured after the decay of their parents (filiation).

Cross section calculation

The cross section values are calculated using the wellknown activation formula, defined as a relative equation in which the knowledge of the beam current is no longer necessary thanks to the recommended reactions. The uncertainty is expressed as a propagation error calculation (see [9] for more details).

Catcher correction

In our experiments devoted to medical isotope production studies, only one catcher was placed and was necessary to collect the recoil nuclei after the target. However, in the case of fission, the FP are emitted both forward and backward direction from the target. The use of a catcher before the target should then be mandatory to collect all the nuclei. In order to correct it, the proportion of FP emitted backward the target has been estimated by kinematic and Monte Carlo calculations.

In our calculation, the projectile is sent on a thorium target followed by a catcher. The FP emission is computed in the centre of mass of the compound nucleus (Th+p or Th+d): the total kinetic energy is calculated according to the Viola et al. formula [10], both fragments are supposed to be emitted with an energy distributed according to the linear momentum conservation and distributed isotropically in the centre of mass of the reaction. The FP energy is then calculated in the laboratory system and its emission point randomly computed in the depth of the thorium target. Thanks to SRIM calculation [8], one can then determine if the FP remains in the target or exits it, forward or backward. Experimentally we know that FP from Zn-72 to Pm-151 have been detected. The proportion of FP in the target as well as backward and forward the target is determined in the case of Zn-72, I-131 and Pm-151. Fig. reffig1 shows these proportions as a function of the incident energy for proton beams.

Based on the results shown on Fig. 1, one can consider that approximatively the same amount of FP is collected in both catchers. This has been validated with a dedicated experiment with two catchers, forward and backward (see Fig. 4). Thus, the amount of FP activity measured in the (forward) catcher is simply multiplied by a factor 2 and added to the activity obtained in the target to determine the cross section values discussed hereafter. The small difference due to kinematics between forward and backward visible on Fig. 1 is neglected compared to the uncertainty on the experimental determination of the activity on the catcher.

The TALYS code

In this work, all the experimental cross section values are compared with the version 1.6 of the TALYS code released in December, 2013 [11]. TALYS is a nuclear reaction program which simulates reactions induced by light particles on target nuclei heavier than carbon. It incorporates theoretical models to predict observables including cross section values as a function of the incident particle energy (from 1 keV to 1 GeV). A combination of models that best describes the whole set of available data for all projectiles, targets and incident energies have been defined by the authors and put as default in the code. In this way, a calculation can be performed with minimum information in the input file: the type of projectile and its incident energy, the target type and its mass.

Since there are some differences between experimental data and the results of the TALYS code using default models, we have defined a new combination of models, already included in the TALYS code. The description of the optical, preequilibrium and level density models have been found to have a great influence on the calculated production cross section values. Better results are, in general, obtained when proton and deuterons are used as projectile using the optical model described, respectively by [12] and [13]. And for both projectiles, but also with alpha particles [9], when a preequilibrium model based on the exciton model including numerical transition rates with optical model for collision probabilities [14] and a model for the microscopic level density from Hilaire's combinatorial tables [15] is used. The results referenced as TALYS 1.6 Adj in Figs. 4 and5 correspond to TALYS calculations performed with this new combination of models.

Results and discussion

The mass distributions

The data depicted in Fig. 2 and Fig. 3 show the cumulative cross section values of the detected FP, respectively with protons and deuterons as projectiles, plotted as a function of the FP mass, A. The points correspond to our measurements and the dash lines are spline fits on these experimental values. The full lines to TALYS results with default models, selecting only FP that have been detected in our experimental conditions (presented at the end of the part 2.1). The apparent depletion of the heavier FP peak is due to these experimental conditions (γ spectrometry detectability).

Proton results (Fig. 2) show that the symmetric fission (highlighted around A = 110) becomes more probable than the asymmetric fission with the increase of the incident energy (above 30 MeV). This has already been observed with protons [16] and neutrons [17]. Figure 2 shows saturation, even decrease of the amplitude of the symmetric fission above 57 MeV. This trend has to be confirmed and explained. The cross section values are similar in the case of proton and deuteron irradiations (see Fig. 3) considering the same incident beam energy. In all cases, the TALYS code with default models is not able to reproduce our experimental data leading to a too low production of fission fragments.

The Th-232 total fission cross sections

Figure 4 and Fig. 5 show the sum of the FP production cross sections measured in our experiments, for protons and deuterons as projectiles, respectively. These plots also show the TALYS code results (limited to the detected FP) and the TALYS Th-232 total fission cross section (including all FP). For both projectiles, protons (Fig. 4) and deuterons (Fig. 5), our experimental total fission value reaches 800 mb for the incident energy of 30 MeV. Both curves show the same trend.

Figure 4 presents also the results from an experiment made with two catchers; one backward and one forward the target. Two energy points (69.9 MeV and 57.1 MeV) are in agreement with values obtained when only one catcher is placed forward relative to the beam direction and its contribution is doubled. This confirms our approach for the catcher correction.

In Fig. 4 and Fig. 5 TALYS results using the default models are referred as 'Default' whereas those obtained using the combination of models described in the part 2.4 are refered as 'Adj.'. We found that TALYS Adj. is more able to reproduce the shape and the amplitude of the FP cross sections than TALYS with default models. In addition, the TALYS Adj. results allow estimating the proportion of FP that have been produced during the irradiations but not detected in our experimental conditions (i.e. cooling time, γ spectrometry) by just subtracting the TALYS result without and with our experimental constrains.

Conclusion

Th-232(p,x) and Th-232(d,x) reactions have been studied using the stacked-foils method. These reactions have led to the production of several fission products. The absence of catcher backward the targets in our stacked-foils experiments has been corrected after kinematic and Monte Carlo calculations. These corrections have been validated by a dedicated experiment using two catchers.

Even if our work come from undedicated experiments and give partial measurements, the increase of symmetric fission of Th-232 with the proton and deuteron incident energy has been observed. Total fission cross sections have then been estimated thanks to TALYS calculations. A special effort has been done on the determination of which models included in the TALYS code can better reproduce the experimental data than the default ones, considering our experimental conditions. Our new experimental results could help to contribute to the improvement of the theoretical models. Studies are still in progress on this part.

Figure 1 .

 1 Figure 1. Fraction of fission products exiting the thorium target (40 µm thick), backward and forward, as a function of the incident proton energy.

Figure 2 .

 2 Figure 2. Mass distribution of the fission products for protons as projectiles.

Figure 3 .

 3 Figure 3. Mass distribution of the fission products for deuterons as projectiles.

Figure 4 .

 4 Figure 4. Th-232 total fission cross section for protons as projectiles.

Figure 5 .

 5 Figure 5. Th-232 total fission cross section for deuterons as projectiles.

c The Authors, published by EDP Sciences. This is an Open Access article distributed under the terms of the Creative Commons Attribution License 4.0 (http://creativecommons.org/licenses/by/4.0/).

The ARRONAX cyclotron is a project promoted by the Regional Council of Pays de la Loire financed by local authorities, the French government and the European Union. This work has been, in part, supported by a grant from the French National Agency for Research called "Investissements d'Avenir", Equipex Arronax-Plus No. ANR-11-EQPX-0004 and Labex No. ANR-11-LABX-0018-01.


Laurent Pizzagalli
email: laurent.pizzagalli@univ-poitiers.fr

Pierre Beauchamp

Hannes Jónsson

H J Ónsson

Calculations of dislocation mobility using Nudged Elastic Band method and first principles DFT calculations

come

Introduction

Dislocations modelling has first been developed in the framework of linear elasticity theory, not only because of historical reasons, but also because this approach has been able to capture most of their physical properties [START_REF] Hirth | Theory of dislocations[END_REF]. However, this approach is inadequate to deal with the severe lattice distortion near the core of a dislocation. In consequence, it is not suited for describing important properties governed by the structure of the core, such as dislocation mobility or dislocation interactions at short distances. Major improvements have been obtained thanks to the Peierls-Nabarro (PN) model and its recent developments [START_REF] Lu | I: Methods and Models[END_REF][START_REF] Joós | [END_REF][4], since it combines elasticity theory with the discrete nature of the lattice. This model removes the artificial divergence at the core originating from the continuum description, allowing to determine several key quantities such as the core width, the Peierls energy and the Peierls stress. Nevertheless, while this approach has be proved to be successful in many cases, it is not accurate enough for several classes of materials such as covalent or ionic materials, in which dislocation cores are very narrow or reconstructed. Although the PN model has been recently improved [4], it intrinsically includes some limitations by construction [5].

A full atomistic treatment of the dislocation is expected to yield better results. Hence, accurate values can be determined in case of f.c.c. metals, for which reliable classical interatomic potentials are available. However, for other materials such as semiconductors or ceramics, energetic properties of dislocation cores are poorly described by potentials. This point has been brought to light in the case of silicon for example [6]. Another option is to use first principles density functional theory (DFT) calculations which give an accurate description of most materials even those including unusual and unfavorable bonding configurations. The problem with these methods, however, is that the simulated systems are limited in size, and as a consequence boundary conditions require careful and special treatments [7][8][9], or the use of multiscale modelling methods [10,11]. This is especially true for Peierls stress determination. Recently, Pizzagalli

Position

Peierls energy

E p E f 2 E f 1 E p E f 1 E p

Periodic conditions

Fixed conditions

Figure 1. (color online) Effect of boundary conditions on Peierls energy variation as a function of the dislocation core position. For fixed conditions, the interaction between the single dislocation and the surface cluster ranges from E 1

f to E 2 f (dash-dotted line). By combining periodic conditions and NEB, and by displacing the dislocations coherently, the interaction Ep between dislocations remains constant (dashed line) and does not affect the deduced Peierls stress and activation energy.

and Beauchamp have calculated the Peierls stress of screw dislocation in silicon using reduced systems and both periodic and fixed conditions [12]. This work shows that Peierls stress determination using first principles accuracy is possible, but that corrections are required to compensate for boundary effects.

In this paper, we propose a method allowing an accurate determination of the Peierls stress, without the need for boundary corrections. Also, the calculations provide the full energy curve for the movement of the dislocation through lattice, i.e. the whole Peierls barrier, as well as the associated core structure configurations. The method is based on the combination of the nudged elastic band (NEB) method for finding a minimum energy path (MEP) [START_REF] Jónsson | Classical and Quantum Dynamics in Condensed Phase Simulations[END_REF] of a transition and first principles density functional theory (DFT) calculations of the atomic forces subject to periodic boundary conditions. Previously, the NEB method has been used in a few studies of dislocation mobility but only in combination with empirical potentials and fixed atom boundaries in finite systems [START_REF] Rasmussen | [END_REF][15][16][17]. By using periodic boundary conditions, as we do here, smaller systems can be used in the calculation and it becomes easier to make use of first principles DFT methods. In the following we first describe errors coming from fixed boundaries in small system, and we show how this effect is greatly reduced using NEB and periodic boundaries. Then, the Peierls energy and stress are determined with first principles accuracy in the case of a screw dislocation in silicon, moving along two different directions. Finally we discuss the advantages and the limitations of this method, compared to previous models.

The straightforward approach to computational studies of dislocations is to construct a system that contains one dislocation and fixed atoms at the boundaries of the system. In that case, boundary conditions will introduce errors in the energy variation as a function of dislocation position, especially in the case of small systems. This is schematically depicted in figure 1. The initial and final configuration states include the dislocation-surface interaction energy E 1 f , positive (negative) for a fixed (free) surface. For the transition state, ideally located in the center of the system, the interaction energy E 2 f is different, i.e. E 2 f = E 1 f . As a result, the calculated curve will be different from the real one (Fig. 1). For a screw, the energy change ∆E due to fixed cylindrical boundaries is proportional to the square of the ratio between the dislocation distance from the cell center and the cell size, with a factor depending on the shear modulus and the anisotropy of the material [18]. For silicon and cluster sizes considered in first principles DFT calculations, ∆E is about 50 meV. This effect could be reduced by considering larger cluster systems [17].

Here, we propose to use periodic boundary conditions, which is well suited for standard plane-waves first principles DFT calculations, and as such have been widely used for investigating dislocations [19][20][21][22][23][24]. The initial configuration is a system including an even number of dislocations, with a zero Burgers vectors sum (Figure 1). The final configuration is then obtained by a simple translation equal to the investigated dislocation displacement. Finally, a standard NEB calculation between these two configurations is performed. The important point here is that distance between dislocations remain constant along the whole path. This condition is fullfilled at the beginning, if the initial MEP is built from linear interpolation between the end point configurations. We have found that if the springs linking adjacent replicas in the NEB are chosen to be stiff enough and a direct minimization method is used to relax the system to the MEP closest to the initial guess, then the the distance remains constant during the full relaxation until convergence has been reached. The advantage of the method is illustrated in figure 1. Since distances between all dislocations remain constant throughout the NEB calculation, the energy variation as a function of position is simply shifted by a constant term, provided that the total interaction between dislocations E p is elastic. As a result, the shape of the curve is the same as the real one.

Model and methods

We apply this method to a non-dissociated screw dislocation in silicon, with an initial stable shuffle configuration [6]. This choice has two motivations. First, in a previous study, the Peierls stress has been determined with first principles DFT and corrections due to periodic and fixed boundary conditions [12]. Second, in the same work, it has been shown that interatomic potentials for silicon are not reliable enough for this particular problem. We have then an interesting system, for which first principles accuracy is required, and with available data for comparison. We considered a system similar to that used in [12], including a dipole of dislocations, a cell geometry of 12 × 12 × 1 (144 atoms) corresponding to an infinite quadrupolar arrangement [12]. First principles DFT calculations have been performed with the VASP distribution [25,26], with 2 special k-points along the dislocation line direction, an ultrasoft pseudopotential for Si [27], a plane waves energy cutoff of 140 eV, and the GGA functional PW91 [START_REF] Perdew | Electronic Structure of Solids '91[END_REF]. System relaxation is stopped when all forces are lower than 3×10 -3 eV Å-1 . With these parameters, a lattice constant of 5.475 Å and a bulk modulus of 99.7 GPa have been obtained, in good agreement with experiments [START_REF] Hull | Properties of crystalline silicon[END_REF]. The NEB calculations have been performed with both improved tangent and climbing images algorithms [START_REF] Henkelman | [END_REF]31]. Figure 2 shows the zinc-blende structure, and two dislocation displacements, investigated in this work. D' corresponds to the easiest path for displacing a non-dissociated screw, while D" is the direct path allowing the screw to cross-slip from shuffle to glide {111} planes.

In order to obtain the Peierls stress from the energy variation, an accurate determination of the dislocation position relative to the lattice is required. Distances between images computed in configuration space, as output by NEB, are useless in that case because the dislocation position does not linearly follow the configuration space coordinates. As stated in [32], dislocation centers can not be derived uniequivocally from the atomic positions. Dislocation centers are then determined according to a well-known procedure. Relative displacements of atoms on both sides of the plane perpendicular to { 101} and including the dislocation displacement vector are first calculated. Figure 2 shows the bonds linking these atoms for D' and D". Then, the dislocation centers are determined by fitting relative displacements using the elastic expression b/2π arctan((x -x 0 )/∆), ∆ allowing to define a dislocation core width [START_REF] Hirth | Theory of dislocations[END_REF]. This approach assumes that the whole dislocation can still be described as a Peierls dislocation throughout its displacement in the lattice. As a check, the distance between the two dislocations in the dipole has been determined for all relaxed NEB images. The maximum deviation is 3 × 10 -2 Å, i.e. less than 1% of the total displacement. This shows that the two dislocations have been coherently displaced, as assumed.

Results

The Peierls energy, i.e. the energy variation associated with the displacement of a dislocation in a crystalline lattice, is represented in figure 3 for two paths, D' and D". Note that a multiplying factor of 1/2 has been applied since two dislocations are included in our cell. For both paths, starting from the stable configuration, the energy is smoothly increasing until a halfway maximum, then decreasing to the next stable configuration. In the initial configuration A, the bond linking atoms 1 and 2 is markedly stretched (Fig. 2 and Fig. 4). When the dislocation moves along D', this stretching intensifies until the bond breaks at the saddle configuration (Fig. 4), which is a high symmetry position at the intersection of shuffle and glide planes. A simple force relaxation of this configuration indicates it is metastable. In previous, DFT/LDA calculations, we have found that this configuration was weakly stable [6], with an energy difference of 0.32 eV per Burgers vector, in good agreement with the value of 0.41 eV/b obtained in this work. Finally, atoms 1 and 2 form bonds with neighbor atoms to displace the dislocation core in A'. The D' migration is then done by the breaking and reforming of a stretched bond, explaining why this mechanism is the easiest. For D", we found that the migration mechanism is totally different. Initially, the bond between atoms 2 and 3 is slightly stretched and out of the (101) plane (Fig. 2 and Fig. 4). This bond is progressively twisted when the dislocation moves along D", until it belongs to the (101) plane at the saddle point (Fig. 4). This corresponds to a screw in the glide set, with a sp 2 hybridized core [6]. Actually, it is rather a weakly stable configuration, not seen here due to our limited NEB images number. We have then two totally different mechanisms, in one case, the dislocation is displaced via bond breaking and reforming in one case (D'), and via severe bond twisting in the other (D"). It is also interesting to compare calculated Peierls energies with classical models. In fact, it is common to use a phenomenological expression W sin 2 (x) for approximating the Peierls energy [START_REF] Hirth | Theory of dislocations[END_REF]. Introducing our calculated maximum Peierls energy in this simple model, we obtain a correct estimation of our computed energy curves, in particular for the D' path.

The Peierls stresses can be accurately obtained as the extrema of the derivative of the spline-fitted energy variations (Figure 3), divided by b 2 , b being the Burgers vector. For the D' path, the Peierls stress is 2.47 eV Å-3 , in excellent agreement with the value 2.6 ± 0.2 eV Å-3 obtained from DFT calculations of sheared systems combined with boundary corrections [12]. For the D" path, the calculated value is 4.05 eV Å-3 . Obviously, comparing energy barriers and Peierls stresses shows that dislocation displacement would be easiest along the D' path, so that the dislocation remains in the shuffle plane. Shuffle-glide transition by dislocation cross-slip along the D" path is then unlikely. This is also in agreement with previous calculations showing that even if a shear stress favoring the D" displacement is applied on the system, the dislocation is always displaced along the D' path [12]. Figure 3 also represents variation of the core width ∆ as a function of the dislocation displacement in the lattice. In both cases D' and D", the core is widened, with a maximum increase of about 20%, corresponding approximately to the maximum Peierls energy. The D" curve shows a local minimum compared to the D' curve. A possible explanation is related to the weak stability of the transition configuration along the D" path, for which the dislocation core is slightly narrowed, whereas for D', the transition configuration is a true saddle point.

sible paths, such as the D" path corresponding to a shuffle-glide cross-slip. We have shown recently that the latter could not be obtained by simply stressing the system [12]. Also, it has been shown recently that the more stable structure of a non-dissociated screw dislocation in the glide set exhibits a 2 × 1 reconstruction along the dislocation line [33]. Using our approach, it should be possible to investigate a transition between the single-period shuffle to the double-period glide screw dislocation, provided that the initial system dimension along the dislocation line is 2b, and that the reconstruction path is simple enough to be obtained with a small number of images in the NEB band.

In our approach, what is really calculated is the coherent displacement of a periodic arrangement of dislocations in a crystalline lattice. One may wonder whether this situation is equivalent to the displacement of a single dislocation. Obviously, even if the distance between dislocations remains constant, the core structure, and then the dislocation core field, is varying during the displacement. So the interaction between dislocations is not constant during the displacement. However, this effect is expected to be negligible and lower than surface-dislocation interactions present in fixed boundary calculations, even in small cell. This is the same assumption done in standard calculations when comparing the stability of several dislocation core configurations, where a constant elastic interactions is assumed between dislocations, independently of the core. In order to get an estimation of this effect, we have performed dislocation calculations with classical potentials and periodic boundary conditions, comparing energy differences between two core configurations as a function of cell size. In all cases, we have found that for the cell size considered here, i.e. 12 × 12, the core field interaction energy is lower than 1% of core energy differences. The error is then in the range of standard first principles accuracy, and can be safely neglected.

Another question concerns the fact that during the NEB calculations the coherency of the collective displacement is retained. In fact, it is likely that it is energetically more favorable to first displace one dislocation in the dipole system, then the other. This is especially true in hard materials with covalent bonds for which a strong interaction is expected between close dislocations. However, we usually found in our calculations for silicon that the dislocations were coherently displaced. The key explanation lies in the the fact that a relaxation algorithm such as steepest descent or the velocity Verlet algorithm with damped dynamics [START_REF] Jónsson | Classical and Quantum Dynamics in Condensed Phase Simulations[END_REF] brings the NEB to the minimum energy path that is nearest to the initial guess even though this may not be the global minimum energy path. While such convergence to the nearest local minimum rather than the global minimum is often a problem when dealing with systems represented by a complex energy surface, it is an advantage here because one can select out the right minimum energy path by choosing the initial path correctly. In order to keep the coherency, it is important that initial and final configurations in the band are well relaxed, with symmetric dislocation cores. We found that if this condition is not fulfilled, a non-coherent displacement could occur. It is also better to use systems including two dislocations instead of four, since it is easier to keep the coherency in the former case. In NEB calculations, we used standard spring constants (i.e. 5.0 eV Å-2 ). It is possible that using slightly stiffer springs would also help for keeping coherent displacements.

In conclusion, we propose a new method for investigating the mobility of dislocations using periodic boundary conditions and first principles accuracy. With this method, the Peierls energy and stress, and the dislocation structure during the displacement could be investigated, as well as high energy barriers, out of reach of standard calculations. In principle, the approach could be used for all kind of dislocations and for all materials, as long as dislocations in the computational cell are stable towards each other.

Figure 2 .

 2 Figure 2. (color online) Ball-and-stick representation of the zinc-blende structure (( 101) plane). The path D' (D") is defined by a dislocation displacement from A to A' (A to A"). Dashed lines show bonds linking atoms used for calculating the disregistry due to the dislocation core. In particular, atoms located in dislocation cores for saddle configurations are marked by numbers (1 and 2 for D', and 2 and 3 for D").

Figure 3 .

 3 Figure 3. (color online) Energy (in eV/Burgers vector, top graph), stress (derivative of the plotted energy curves, middle graph), and core width variation (bottom graph) as a function of the dislocation core position in the lattice for paths D' (full lines) and D" (dashed lines).

Dislocation mobility using NEB and DFT calculations

Acknowledgements

The authors are grateful to Dr. Ladislas Kubin for a critical reading of the manuscript. One of us (L.P.) thanks the Poitou-Charentes region and the University of Iceland Research Fund for financial support during his stay in Reykjavik.

Concluding discussion

Our approach combining NEB and first principles DFT calculations allows to determine the energy of a dislocation translated in a crystalline lattice, as well as the associated Peierls stress, without the need for any corrections. In particular, the calculated Peierls stress for a screw dislocation in silicon is in excellent agreement with published results [12]. In addition to the Peierls energy determination, our approach has several other advantages compared to previous models. In particular, we have access to the dislocation structure during its displacement. Another advantage is the possibility of easily investigating several pos-