Vincent Danos

Tobias Heindel

Ilias Garnier

Jakob Grue Simonsen

Continuous-time Markov chains as transformers of unbounded observables

This paper provides broad sufficient conditions for the computability of time-dependent averages of stochastic processes of the form f (Xt) where Xt is a continuous-time Markov chain (CTMC), and f is a real-valued function (aka an observable). We consider chains with values in a countable state space S, and possibly unbounded f s. Observables are seen as generalised predicates on S and chains are interpreted as transformers of such generalised predicates, mapping each observable f to a new observable Ptf defined as (Ptf )(x) = Ex(f (Xt)), which represents the mean value of f at time time t as a function of the initial state x. We obtain three results. First, the well-definedness of this operator interpretation is obtained for a large class of chains and observables by restricting Pt to judiciously chosen rescalings of the basic Banach space C0(S) of S-indexed sequences which vanish at infinity. We prove, under appropriate assumptions, that the restricted family Pt forms a strongly continuous operator semigroup (equivalently the time evolution map t → Pt is continuous w.r.t. the usual topology on bounded operators). The computability of the time evolution map follows by generic arguments of constructive analysis. A key point here is that the assumptions are flexible enough to accommodate unbounded observables, and we give explicit examples of such using stochastic Petri nets and stochastic string rewriting. Thirdly, we show that if the rate matrix (aka the q-matrix) of the CTMC is locally algebraic on a subspace containing f , the time evolution of projections t → (Ptf )(x) is PTIME computable for each x. These results provide a functional analytic alternative to Monte Carlo simulation as test bed for mean-field approximations, moment closure, and similar techniques that are fast, but lack absolute error guarantees.

Introduction

The study of properties of stochastic processes is a currently highly active field of research in computer science, but has been avidly studied in mathematics even prior to Kolmogorov's development of the axiomatic approach to probability. For the particular case of continuous-time Markov chains (CTMCs), the crucial property of interest is how it acts on a function from its underlying state space to the reals -a "generalised predicate" in analogy to the usual notion of predicate transformers in computer science [Koz83]. Here, we are interested in the problem of computing time-dependent expectations of real-valued functions of the state space (so-called observables) on a continuous-time Markov chain X t . Ultimately, one wants provably correct, practical algorithms that given some continuous-time Markov chain X t and some observable f computes the transient conditional mean t → (x → E x (f (X t ))), i.e., the average of f conditioned on the initial state as a function of time. However, it is not a priori clear that this function is even computable let alone practically so; hence, it must be established that there is an algorithm that computes the transient conditional mean to arbitrary desired precision. As the basic objects are typically (computable) real numbers, we will employ the type 2 theory of effectivity (TTE) for computability.

It is useful to recall the nature of the objects that we manipulate and want to compute in the finite case. A continuous-time Markov chain (CTMC) on a finite state space S is entirely captured by its q-matrix, which is an S × Sindexed real matrix in which every row sums to zero and all negative entries lie on the diagonal. Any finite q-matrix Q induces a semigroup t → P t = e tQ which describes the time evolution of the CTMC. For any π on S (viewed as a row vector), the map t → πP t corresponds to the probability of being at any given state at time t for the initial probabilistic state π. Therefore, if X 0 is distributed according to π, the associated stochastic process X t is distributed according to πP t . Dually, one can interpret a CTMC as a transformer of observables, using the discrete time setting [Koz83] as an analogy. For any f : S → R (seen as a column vector), P t f = (x → y∈S p t,xy f (y)) is the vector of conditional means3 of f at time t as a function of the initial state. The function we seek to compute is precisely t → P t f . Formally, P t f is the (column) vector of conditional means (i.e., (P t f )(x) = E x (f (X t )) for x ∈ S) where X t is the CTMC with transition function P t . For the finite state case, t → P t f = e tQ f is the unique solution to the initial value problem

d dt u t = Qu t u 0 = f (1)
and this can be solved by any method for computing matrix exponentials or numerical solvers for finite ODEs. In the countably infinite state case, observables are possibly unbounded-this makes computing transient conditional means fundamentally harder than in the finite case and calls for more sophisticated mathematics, as we shall describe in § 3.2.

For motivation, we give two paradigmatic examples of transient means. The classic example of a stochastic process of the form f (X t ), i.e., a pair of a CTMC X t and an observable f , where X t models a set of chemical reactions with the molecule count of a certain chemical species as observable: states of X t are multisets over a finite set of species; then, we are interested in computing the evolution of the mean count of a certain species.

Another natural example native to computer science is a stochastic interpretation of any string rewriting system as a CTMC X t . An obvious class of observables for string rewriting are functions f that count the occurrence of a certain word as sub-string in each state of the CTMC X t ; note that this is different from counting "molecules", as there is always only a singe word! For example, consider the string rewriting system with the single rule a aba and initial state a: the mean occurrence count of the letter a grows as the exponential function e t while the mean occurrence count of the word aa is zero at all times; adding the rule ba ab does not change the mean a-count but renders the mean count of the word aa non-trivial.

Note that these two classes of models are only meant as simple examples and the results of this paper are not geared towards any particular modelling language for CTMCs; the results apply equally well to Kappa models [DFF + 10] or stochastic graph transformation. 4Our main result is a sufficient condition for E x (f (X t )) to be computable, powerful enough to encompass many interesting unbounded observables. We proceed by construction of a suitable computable Banach space where to restrict our operators P t , in such a way that they form a strongly continuous semigroup (SCSG). It follows that the initial value problems corresponding to observables in the domain of the generator of the semigroup admit computable solutions. The construction of our Banach space rests on functional analytic techniques recently developed in Ref. [Spi12,Spi15], which we combine with the work on SCSGs by Weihrauch & Zhong [WZ07].

We do not assume that the reader has some acquaintance with continuoustime Markov chains on a countable state space and recall the basic concepts of transition functions and q-matrices (Def. 5 and 6) in the preliminaries. Having said that, and both in order to save space and to keep the logic of the paper clearly apparent, most of the mathematical material used is postponed to a series of appendices.

Two motivating examples of CTMCs with observables

We illustrate our constructions with: (i) chemical reaction networks (CRN), aka stochastic Petri nets, and (ii) stochastic string rewriting as a simple example of (rule-based) modelling. In both cases, the construction of the q-matrix implied by a model is readily done, and so is the definition of a natural set of unbounded observables with clear relevance to the dynamics of a model: word occurrence counts for stochastic string rewriting (Def. 2) and multiset inclusions for Petri nets (Def. 3).

Stochastic string rewriting and word occurrences

Stochastic string rewriting can be thought of as never ending, fair competition between all redexes of rules, "racing" for reduction; the formal definition is as follows, in perfect analogy to Ref.

[HLM06] which covers the case of graphs.

Definition 1 (Stochastic string rewriting). For each rule ρ = l r ∈ Σ + × Σ + we define the q-matrix of ρ, denoted by Q ρ , as the q-matrix Q ρ = (q ρ uv ) u,v∈Σ + on the state space of words Σ + with off-diagonal entries

q ρ uv = (w, w ) ∈ Σ * × Σ * u = wlw , v = wrw
for each pair of words u, v ∈ Σ + such that u = v, and diagonal entries,

q ρ uu = -v =u q ρ uv for all u ∈ Σ + . For a finite set of rules R ⊆ Σ + × Σ + , we define Q R = ρ∈R Q ρ ,

and with additional choices of rate constants

k : R Q + , we define Q R,k = ρ∈R k ρ Q ρ .
For a given rule set R, each entry q R uv of the q-matrix corresponds to the propensity to rewrite: it is just the number of ways in which u can be rewritten to v. We shall usually work without rate constants for the sake of readability. Note that the use of Σ + for the left and right hand side of rules is convenient to get string rewriting as a special case of graph transformation in a straightforward manner.

The occurrence counting function of a word as sub-string in the state of the CTMC of R is as follows.

Definition 2 (Word counting functions). Let w ∈ Σ + be a word. The wcounting function, denoted by w :

Σ + → R ≥0 , maps each word x ∈ Σ + to w (x) = |{(u, v) ∈ Σ * × Σ * | x = uwv}|.

Stochastic Petri nets and sub-multiset occurrences

We recall the definition of stochastic Petri nets and occurrence counting of a multisets. Note that for the purposes of the present paper, places and species are synonymous.

Definition 3 (Multisets and multiset occurrences).

A multiset over a finite set P of places is a function x : P → N that maps each place to the number of tokens in that place. Given a multiset, x ∈ N P , the x-occurrence counting function x : N P → N is defined by

x (y) = y! (y-x)! x ≤ y 0 otherwise
where z! = p∈P z(p)! is the multiset factorial for all z ∈ N P .

Definition 4 (Stochastic Petri net). Let P be a finite set of places. A stochastic Petri net over P is a set

T ⊆ N P × R >0 × N P
where N P is the set of multisets over P, which are called markings of the net; elements of the set T are called transitions. The q-matrix Q l,k,r on the set of markings for a transition (l, k, r) ≡ l → k r ∈ T has off-diagonal entries

q l,k,r xy = k • l (x) l (x) > 0, y = x -l + r 0 otherwise
where addition and subtraction is extended pointwise to

N P . The q-matrix of T is Q T = (l,k,r)∈T Q l,k,r .

Preliminaries

For the remainder of the paper, we fix an at most countable set S as state space.

Transition functions and q-matrices

We first recall the basic definitions of transition functions and q-matrices. We make the usual assumptions [And91] one needs to work comfortably: namely that q-matrices are stable and conservative and that transition functions are standard and also minimal as described at the end of § 3.1 (and detailed in Appendix A). With these assumptions in place, transition functions and q-matrices determine each other, and one can freely work with one or the other as is most convenient. z∈S p t,xz ≤ 1 for all x ∈ S and t ∈ R ≥0 . Thus, each row of a transition function corresponds to a sub-probability measure, and transition functions converge entry-wise to the identity matrix at time zero. Unless stated otherwise, all transition functions in this paper are standard.

Taking entry-wise derivatives at time 0 gives a q-matrix.

Definition 6 (q-matrix). A q-matrix on S is an S × S-matrix Q = (q xy ) x,y∈S with real entries q xy such that q xy ≥ 0 (if x = y), q xx ≤ 0, and z∈S q xz = 0 for all x, y ∈ S.

Conversely, for each q-matrix, there exists a unique entry-wise minimal transition function that solves Equation (2) [And91, Theorem 2.2],

d dt P t = QP t , P 0 = I (2)
which is called the transition function of Q. From now on, we assume that all transition functions are minimal solutions to Equation (2) for some q-matrix Q (cf. Definition 15, see also [Nor98,p. 69]).

The Abstract Cauchy problem for P t f

Abstract Cauchy problems (ACPs) in Banach spaces [Ein52], are the classic generalisation of finite-dimensional initial value problems. (See Appendices C and D for details.) Concretely, we want to obtain P t f as unique solution u t of the following generalisation of our earlier equation (1):

d dt u t = Qu t (t ≥ 0) u 0 = f (acp)
where f is an observable and Q is a linear operator which plays the role of the q-matrix.

There are a few points worth noting. First, the topological vector space of all observables R S cannot be equipped with a suitable complete norm. Therefore, one has to look for a Banach subspace B ⊂ R S wherein to interpret the above equation. There are several ways to do this, and they are not equally interesting. Second, as P t f = u t is the desired solution, and P 0 = I, it follows that d dt P t f | t=0 = Qf . If this derivative does not exist, Qf is simply not defined. In fact, as is clear from the examples, we can only expect Q to be partially defined and an unbounded operator.5 On the positive side, if P t is a strongly continuous semigroup (Def. 20) on B, meaning lim h 0 P h f = f for all f in B, and we take Q as its generator defined on a subspace of B (Def. 21) by the above formula, then P t f is the unique solution of (acp), provided f is in the domain of Q [EN00, Proposition II.6.2].

Even better, in this case, not only does (acp) have P t f as unique solution, but we get an explicit approximation scheme:

P t f = lim n→∞ e tAn f ( 3
)
where θ is a constant of the SCSG such that nI -Q is invertible for n > θ and

A n = nQ(nI -Q) -

Banach space wanted!

Table 1 gives an overview of initial value problems for transient distributions (first row) and transient conditional means (second row). Transient distributions are summable sequences, and transition functions form SCSGs [Reu57] and therefore allow for a well-posed corresponding ACP. But the classic example of a Banach space to reason about conditional means [RR72] is the space C 0 (S) of functions vanishing at infinity, i.e., functions f : S → R such that for all > 0, the set {x ∈ S | f (x) ≥ } is finite, equipped with the supremum norm (Def. 17). The corresponding processes are called Feller transition functions [And91, § 1.5] and verify a principle of finite velocity of information flow (for all t, y, p t,xy vanishes as x goes to infinity).

Spieksma's theorem

A solution is provided by a result of Spieksma [Spi12, Theorem 6.3], which constructs a Banach space for a given q-matrix Q and an observable f of interest. The formulation will be in terms of so-called drift functions.

Definition 7 (Drift function).

Let Q be a q-matrix on S, and let c ∈ R.

A function W : S → R >0 is called a c-drift function for Q if for all x ∈ S (QW )(x) := y∈S q xy W (y) ≤ cW (x).
We will say that W is a drift function for Q if there exists c ∈ R such that it is a c-drift function for Q. One can show that P t W ≤ e ct W in this case. Thus, drift functions control their own growth under the transition function.

We also need to work with weighted variants of C 0 (S).

Definition 8 (Weighted C 0 (S)-spaces). Let S be a set and let W : S → R >0 be a positive real-valued function, referred to as a weight. The Banach space C 0 (S, W ) consists of functions f : S → R such that f /W vanishes at infinity, where

( f /W )(x) = f (x) /W (x). The norm _ W on such functions f is f W = sup x∈S | f (x) /W (x)|.
As C 0 (S, W ) is isometric to C 0 (S) it is indeed a Banach space. It is also a closed subspace of L ∞ (S, W ), the set of functions such that f /W is bounded. We will use later the fact that: Lemma 1. Finite linear combinations of indicator functions,7 form a dense subset of C 0 (S) Indeed, suppose S = N for convenient notations, and define g n = f truncated after n. It is easy to see that g n converges to f in the _ W norm. In fact,

g n -f W = sup x>n |f (x)|/W (x) → 0 iff f /W vanishes at infinity. Now the idea is to find positive drift functions V, W for Q such that V ∈ C 0 (S, W ), i.e.
, such that V /W vanishes at infinity. Thus, both functions can grow at most exponentially in mean and V is negligible compared to W at infinity. Intuitively, functions on the order of V are as good as functions vanishing at infinity, in the case of Feller processes [RR72], i.e., CTMCs such that transition functions induce SCSGs on C 0 (S).

Theorem 1. Let P t be a transition function on S with q-matrix Q, let V, W : S → R >0 be drift functions for Q such that V ∈ C 0 (S, W ). Then:

1. P t induces an SCSG on C 0 (S, W ). 2. For all f ∈ C 0 (S, W ), P t f is given by Equation (3) in C 0 (S, W ) where Q is the generator of P t .
The first part of the theorem is proved in [Spi12, Theorem 6.3]; the second part follows from the general theory of ACPs. Note that f does not need to be in the domain of Q, in which case we only obtain a 'mild' solution to the ACP, i.e., a solution to its integral form which might not be everywhere differentiable. If on the other hand, f is in the domain of Q, then P t f is continuously differentiable.

We have now covered the mathematical ground needed to make sense of conditional means of not necessarily bounded observables. This, however, does not immediately yield an algorithm for computing transient means, as it does not yet ensure that transient means are computable even in presence of a complete specification of the underlying system. Even transient conditional probabilities can fail to be computable [AFR11]! Before we proceed to computability questions, let us return to our two classes of examples.

Applications: string rewriting and Petri nets

We now give examples of drift functions for stochastic string rewriting and Petri nets. The former case is well-behaved since the mean letter count grows at most exponentially. The case of Petri nets will be more subtle and we shall give an example of an explosive Petri net such that we can nevertheless reason about conditional means of unbounded observables.

For string rewriting, we have canonical drift functions.

Lemma 2 (Powers of length are drift functions). Let R ⊆ Σ + × Σ + be a finite string rewriting system and let n ∈ N + be a positive natural number. There exists a constant

c n ∈ R >0 such that |_| n : Σ + → R ≥0 is a c n -drift function.
Proof. The proof is deferred to Appendix E. Now, we can apply Spieksma's method to get a Banach space for reasoning about conditional means and moments of word counting functions.

Corollary 1 (Stochastic string rewriting).

Let R be a finite string rewriting system, let n ∈ N \ {0}, and let |_| : Σ + N be the word length function. The transition function

P t of q-matrix Q R is a SCSG on C 0 (Σ + , |_| n ).
Thus, all higher conditional moments of word counting functions can be accommodated in a suitable Banach space. The case of Petri nets is more subtle, since, in general, the (weighted) token count is not a drift function.

Example 1. Consider the Petri net with single transition 2A → 1 3A and with one place A. The token count A is not a drift function. In fact, the corresponding CTMC is explosive (by [Spi15, Theorem 2.1]).

Our final example is an extension of the previous explosive CTMC with a new species whose count can nevertheless be treated using Theorem 1.

Example 2 (Unobserved explosion). Consider the Petri net with transitions

{2A → 1 3A, B → 1 2B}.
The underlying CTMC is explosive, and we cannot apply Theorem 1 to compute the transient conditional mean of the A-count for the exact same reason as in Example 1. However, we can do so for the B-count, using the weight function W = B 2 and observable f = B . Putting V = f allows to apply Spieksma's recipe (ruling out states with B-count 0 for convenience). The conditional mean E 2A+B ( B (X t )) can be best understood by adding a coffin state, on which both the A-and B-count are zero and in which the Markov chain resides after (the first and only) explosion.

Computability

We follow the school of type-2 theory of effectivity, as reviewed in Appendix G following [BHW08]. A real number x is computable iff there is a Turing machine that on input d ∈ N (the desired precision), outputs a rational number r with |r -x| < 2 -d . Next, a function g : R → R is computable if there is a Turing machine that, for each x ∈ R, takes an arbitrary Cauchy sequence with limit x as input and generates a Cauchy sequence that converges to g(x)-where convergence has to be sufficiently rapid, e.g., by using the dyadic representation of the reals.

Computability extends naturally to any Banach space B other than R. We only need a recursively enumerable dense subset on which the norm, addition and scalar multiplication are computable, thus making B a computable Banach space (Definition 30); usually, the dense subset is induced by a basis of a dense subspace. For weighted C 0 -spaces (with computable weight functions) and their duals, we fix an arbitrary enumeration of all rational linear combinations of indicator functions 1 x ; for the Banach space of bounded linear operators on weighted C 0 -spaces we use the standard construction for continuous function spaces [WZ07, Lemma 3.1]. A SCSG P t is computable if the function t → P t from the reals to the Banach space of bounded linear operators is computable. The computable SCSGs correspond to those obtained from CTMCs through Theorem 1.

For simplicity, motivated by the observation that we do not lose any of the intended applications to rule-based modeling, we restrict to row-and columnfinite q-matrices with rational entries in our main result.

Theorem 2 (Computability of CTMCs as observation transformers).

Let Q be a q-matrix on S, let V, W : S → R ≥0 be positive drift functions for Q such that V ∈ C 0 (S, W ). If the q-matrix Q is row-and column-finite, consists of rational entries, and is computable as a function Q : S 2 → Q, and -W : S → Q is computable, the following hold.

1. The SCSG P t , i.e., the function t → P t , is computable.

The evolution of conditional means

(t, f ) → P t f is computable. 3. The evolution of means (π, t, f ) → πP t f is computable for all π : S → R such that x∈S W (x)|π(x)| < ∞.
Proof. We will apply a result by Weihrauch & Zhong on the computability of SC-SGs [WZ07, Theorem 5.4]. Applying this result requires some extra information:

1. the SCSG P t must be bounded in norm by e θt for some positive constant θ; 2. we must have a recursive enumeration of a dense subset of the graph of the infinitesimal generator of the SCSG P t .

We first show that the constant θ, featuring in Theorem 1, i.e., the witness that W is a θ-drift function for Q, satisfies P t ≤ e θt (using the first part of the proof of Theorem 6.3 of Ref. [Spi12]). Next, we obtain a recursive enumeration of a dense subset A ⊆ Q of the domain of the generator Q by applying Q to all rational linear combinations of indicator functions 1 x . Note that for the latter, we use that indicator functions belong to the domain of the generator and

Q1 x = Q1 x (see Proposition 2).
Now, by Theorem 5.4.2 of Ref.

[WZ07], we obtain the first two computability results, as θ, A, 1 is a so-called piece of type IG-information [WZ07, p. 513]. Finally, the third point amounts to showing computability of the duality pairing

_, _ : L 1 (S, W ) × C 0 (S, W ) → R where L 1 (S, W ) is (isomorphic to) the dual of C 0 (S, W ) (Lemma 3).
This theorem immediately gives computability of the CTMCs and (conditional) means for stochastic string rewriting and Petri nets discussed in Sect. 4. Note that the theorem does not assume that V itself is computable; its role is to establish that the transition functions is a SCSG, but V plays no role in the actual computation of the solution. Note also that the algorithms that compute the functions t → P t , (t, f ) → P t f , and (π, t, f ) → πP t f push the responsibility to give arbitrarily good approximations of the respective input parameters π, t and f to the user. This however is no problem for any of our examples or rulebased models in general: t is typically rational, f is computable and even to the natural numbers, and π is often finitely supported or a Gaussian.

Remark 1 (Domain of the generator). The difficulty in working with SCSGs is to find a tractable description of the domain of their generator. Such a description is important for two reasons: first, as said before the theorem, for computability, we need a recursively enumerable dense subset of the graph of the domain of the generator; second, as the unique solution to the ACP is also differentiable, it is useful for numerical integration.

Computability ensures existence of algorithms computing transient means, but yields no guarantees of the efficiency of such algorithms. We now proceed to a special case that (i) encompasses a number of well-known examples, including context-free string rewriting, and (ii) leads to PTIME computability, by reducing the problem of transient conditional means to solving finite linear ODEs.

The finite dimensional case and PTIME via ODEs

We now turn to the special case where we can restrict to finite dimensional subspaces B ⊆ R S . 8 The prime example will be word counting functions and context-free string rewriting systems. Hyperedge replacement systems [DKH97], the context-free systems of graph transformation, can be handled mutatis mutandis. The main result is PTIME computability of conditional means. For convenience, we extend the usage of the term locally algebraic as follows.

Definition 9 (Locally algebraic). We call a q-matrix Q on S locally algebraic for an observable f ∈ R S if the set {Q n f | n ∈ N}, containing all multiple applications of the q-matrix Q to the observable f , is linearly dependent, i.e., if there exists a number N ∈ N such that the application Q N f of the N -th power is a linear combination

N -1 i=0 α i Q i f of lower powers of Q applied to f .
Using local algebraicity of a q-matrix Q for an observable f , one can generate a finite ODE with one variable for each conditional mean

E[Q n f (X t ) | X 0 = x]
(as detailed in the proof of Theorem 3); then, recent results from computable analysis [PG16] entail PTIME complexity.

Theorem 3 (PTIME complexity of conditional means).

Let Q be a qmatrix on S, let x ∈ S, let f : R → S be a function such that f (x) is a computable real and

Q N f = N -1 i=0 α i Q i f for some N ∈ N and computable coefficients α i .
The time evolution of the conditional mean P t f (x), i.e., the function t → P t f (x), is computable in polynomial time.

Proof. Consider the N -dimensional ODE with one variable E n for each n ∈ {0, . . . , N -1} with time derivative

d dt E n (t) = E n+1 (t) if n < N -1 N -1 i=0 α i E i (t) if n = N -1 and initial condition E n (0) = Q n f (x). Solving this ODE is in PTIME [PG16]
(even over all of ≥0 , using the "length" of the solution curve as implicit input).

Finally,

E i (t) = P t Q i f (x) is the unique solution.
Note that the linear ODE that we construct has a companion matrix as evolution operator, which allows to use special techniques for matrix exponentiation [TR03,BCO + 07].

Proposition 1 (Local algebraicity of context-free string rewriting). Let R be a string rewriting system, let w ∈ Σ + , let m ∈ N. The q-matrix Q R of the string rewriting system R is locally algebraic for the m-th power of w-occurrence

counting w m if R ⊆ Σ × Σ + .
Proof. For every product of word counting functions m i=1 wi , applying the q-matrix Q R to this product yields the observable Q R m i=1 wi . Using previous work on graph transformation [DHHZS15], restricted to acyclic, finite, edge labelled graphs that have a unique maximal path (with at least one edge), the observable Q R m i=1 wi is a linear combination k j=1 α j kj l=1 w j,l of word counting functions w j,l (with all

k j ≤ m). Moreover, if R is context-free (R ⊆ Σ × Σ + ), we have kj l=1 |w j,l | ≤ m i=1 |w i |
for all j ∈ {1, . . . , k}. Thus, we stay in a subspace that is spanned by a finite number of products of word counting functions.

Corollary 2. For context-free string rewriting, conditional means and moments of word occurrence counts are computable in polynomial time.

We conclude with a remark on lower bounds for the complexity.

Remark 2. The complexity of computing transient means, even for context-free string rewriting, is at least as hard as computing the exponential function. This becomes clear if we consider the rule a → aa, and the observable of a-counts a . Now, the time evolution of the a -mean conditioned on the initial state to be a, i.e., the function t → E a ( a (X t )), is exactly the exponential function e t . Tight lower complexity bounds for the exponential function are a longstanding open problem [Ahr99].

The main result is computability of transient (conditional) means of Markov chains X t "observed" by a function f , i.e., stochastic processes of the form f (X t ). For this, we have described conditions under which a CTMC, specified by its q-matrix, induces a continuous-time transformer P t that acts on observation functions. In analogy to predicate transformer semantics for programs, this should be called observation transformer semantics CTMCs; formally, P t is a strongly continuous semigroup on a suitable function space. Finally, motivated by important examples of context-free systems -be it the well-known class from Chomsky's hierarchy or the popular preferential attachment process (covered in previous work [DHHZS15]) -we have considered the special case of locally finite q-matrices. For this special case, we obtain a first complexity result, namely PTIME computability of transient conditional means.

The obvious next step is to implement our theoretical results since one cannot expect that the general algorithms of Weihrauch & Zhong [WZ07] perform well for every SCSG on a computable Banach space. For example, the Gauss-Jordan algorithm for infinite Matrices [Par12] should already be more practicable for inverting the operator nI -Q from Equation (3) compared to the brute force approach used by Weihrauch & Zhong [WZ07]. Computability ensures existence of algorithms for computing transient means, but yields no guarantees of the efficiency of such algorithms.

Even if it should turn out that efficient algorithms are a pipe dream -after all, transient probabilities p t,xy are a special case of transient conditional meanswe expect that already implementations that are slow but to arbitrary desired precision will be useful for gauging the quality of approximations of the "meanfield" of a Markov process, especially in the area of social networks [Gle13], but possibly also for chemical systems [SSG15]. Theoretically, they are a valid alternative to Monte Carlo simulation, or even preferable. A From Transition to q-matrices and back The first two properties in Definition 10 make transition functions monoids of real-valued S × S-matrices equipped with matrix multiplication; the last property requires that rows induce sub-probability measures. 10 Thus, an honest transition function P t consists of proper probability matrices, i.e., the entries of every row sum up to 1. Unless explicitly allowed to be non-standard, we shall assume transition functions to be standard, i.e., entry-wise continuous at time 0.

9 The reader familiar with CTMCs might object that in shifting focus from Markov chains to pairs of initial distributions and transition functions, we lose information. However, "[a]ll the probabilistic information about the process, insofar as it concerns only countably many time instants, is contained in the transition function and initial distribution. One could almost say that the transition function is the Markov chain." [And91, p. 4]. 10 A sub-probability measure on an at most countable set S (equipped with the σalgebra of all subsets) can be safely identified with a function to the non-negative reals f : S → R ≥0 that sums up to at most 1, i.e., such that x∈S f (x) ≤ 1 holds.

Example 3 (Yule-Furry process). The Yule-Furry process of parameter β > 0, the archetypal example of a branching process, can be thought of as a simplistic model of population growth; for n 0 ≤ n, its transition probabilities have the following closed form.

p t,n0 n = n -1 n -n 0 (e -βt ) n0 (1 -e -βt ) n-n0
Using the metaphor of population growth, if we start with a population of size 1, the probability at time t to have grown to exactly size n ≥ 1 is e -βt (1 -e -βt ) n-1 .

Typically, transition functions of CTMCs do not have simple closed expressions for their entries; the Yule-Furry process is one of the few exceptions. In contrast, taking (time) derivatives of transition functions (at time zero) typically yields natural functions, intuitively corresponding to a flow of probability mass; for a Yule process of parameter β with initial population n 0 , we think of n 0 independent individuals that can replicate at any given moment in time, and thus the outflow of probability mass is -βn 0 .

Transition functions are always entry-wise differentiable [Kol51,Aus55] and the derivatives at time zero p 0,xy form the so-called q-matrix of the transition function [And91,p. 14]. Provided that all states are stable, i.e., p 0,xx = -∞ for all x ∈ S, q-matrices of transition functions can safely be identified with stable q-matrices according to the following definition [And91, Theorem 2.2].

Definition 11 (q-Matrix [And91, p. 64]). A q-matrix is an S × S-matrix Q = (q xy ) x,y∈S with entries real numbers or -∞ such that

∞ > q xy ≥ 0 (if x = y),
-∞ ≤ q xx ≤ 0, and z∈S q xz ≤ 0 hold for all x, y ∈ S; it is stable (resp. conservative) if q x := -q xx < ∞ (resp.

z∈S q xz = 0) holds for all x ∈ S. In this paper, we assume q-matrices to be stable and conservative. For each q-matrix Q = (q xy ) x,y∈S on an at most countable state space S, the Kolmogorov backward equation

d dt P t = QP t , P 0 = I (4)
has a unique minimal non-negative solution P t = (p t,xy ) x,y∈S (see, e.g., Theorem 2.2 of Ref.

[And91]), which is a possibly dishonest transition function.

Definition 12 (Transition functions of q-matrices). For a stable q-matrix Q on S, the transition function of q-matrix Q, is the unique minimal transition function P t solving the Kolmogorov backward equation, i.e., any other transition function Pt = ( pt,xy ) x,y∈S that solves Equation (4) satisfies pt,xy ≥ p t,xy for all (t, x, y) ∈ R ≥0 × S × S.

In case of a finite state space, all the above is trivial: the transition function P t is just the matrix exponential e Qt , and similarly, if the entries of the q-matrix are bounded. In general, the hard part is to find for a given q-matrix a corresponding transition function, while the opposite direction, passing from a transition function to its q-matrix, is always via differentiation [Kol51,Aus55]; however, with the usual assumptions on CTMCs (as in Norris's standard textbook [Nor98]), qmatrices are stable and conservative, which allows to pass back and forth between transition functions and q-matrices without complications [And91, Theorem 2.2].

B Markov chains

For completeness sake, we recall the classic definitions of continuous-time Markov chains in terms of a stochastic process satisfying the Markov property. However, note that all proofs rely on the functional analytic perspective described in Appendix A, which does not hinge on any non-trivial measure theory. The following definitions are based on [And91].

Definition 13 (Continuous-time Markov chain (CTMC)).

A stochastic process {X t } t∈R ≥0 , defined on a probability space (Ω, F, Pr), with values in a countable set S (to be called the state space of the process), is called a continuoustime Markov chain (CTMC) if for any finite set

0 ≤ t 1 ≤ t 2 ≤ • • • ≤ t n ≤ t n+1
of "times", and corresponding set i 1 , i 2 , . . . , i n-1 , i, j of states in S we have

Pr X tn+1 = j X tn = i, X tn-i = i n-1 , . . . X t1 = i 1 = Pr X tn+1 = j X tn = i whenever Pr X tn = i, X tn-i = i n-1 , . . . X t1 = i 1 > 0.
If for all 0 ≤ s ≤ t and all i, j ∈ S the equation

Pr X t = j X s = i = Pr X t-s = j X 0 = i
holds, the stochastic process {X t } t∈R ≥0 is homogeneous.

As described by Anderson [And91,p. 4], all probabilistic information about a CTMC is given by their transition function.

Definition 14 (Transition function of a CTMC).

Let {X t } t∈R ≥0 be a homogeneous CTMC. The t-dependent matrix P t = (p t,ij ) i,j∈S with

p t,ij = Pr X t = j | X 0 = i for all i, j ∈ S and t ∈ R ≥0 is the transition function of the CTMC.
We only consider homogeneous CTMCs that are stable, standard11 and minimal (cf. Assumption 2.1 of Ref. [Spi12]).

Definition 15 (Stable, standard, minimal). Let {X t } t∈R ≥0 be a homogeneous CTMC. It is stable if the q-matrix of its transition function is stable, it is standard if its transition function is standard, and it is minimal if its transition function is the minimal transition function that solves the Kolmogorov (backward) equation (4).

C Linear operators on Banach spaces

In the present paper, we only consider vector spaces over the field R. We recall some core definitions (see also [EN00, Appendices A and B]).

Definition 16 (Banach space). A Banach space is a complete normed vector space.

We now give examples of Banach sequence spaces, i.e., Banach spaces consisting of functions f : S → R from a fixed set S, which will always be the state space of some CTMC is this paper.

Definition 17 (C 0 -spaces). A function f : S → R vanishes at infinity if for all > 0, f (x) < holds for almost all x ∈ S, i.e., the set {x ∈ S | f (x) ≥ } is finite. The Banach space C 0 (S) is the Banach space that consists of all functions f : S → R that vanish at infinity and its norm _ is given by f

= sup x∈S |f (x)|.
The duals of C 0 -spaces exists and are isomorphic to L 1 -spaces; however, C 0 -spaces are in general not isomorphic to their double dual (which is isomorphic to the corresponding L ∞ -space).

Definition 18 (L 1 -spaces). A function

f : S → R is summable if the sum x∈S |f (x)
| converges (for some order of summation). The Banach space L 1 (S) consists of summable functions and its norm _ is given by f

= x∈S |f (x)|.
Definition 19 (Linear operators, bounded linear operators). Let B be a Banach space. A linear operator on B is a pair (A, D(A)) where D(A) ⊆ B is a linear subspace, called domain, and A is a linear map A : D(A) → B. A linear operator on B is bounded if there exists a constant K ∈ R ≥0 such that A(u) ≤ K u holds for all u ∈ D(A). We denote by L(B) the set of all linear operators (A, D(A)) on B that are bounded and defined on the whole space, i.e., satisfy D(A) = B.

D Abstract Cauchy problems of SCSGs

Every transition function P t = (p t,xy ) x,y∈S (of some CTMC X t with state space S) induces a partial function on the vector space R S that maps a function f : S → R to the function P t f := (x → y∈S p t,xy f (y)) (of the same "type" S → R) whenever the respective sums converge absolutely for all states x ∈ S, and it is undefined on f otherwise. This function is indeed always a linear operator on the vector space of bounded functions [And91, Sect. 1, Lemma 4.4.]. However, we are interested in a stronger property (Theorem 1).

Definition 20 (Strongly continuous semigroup).

Let B be a Banach space, and let L(B) be the set of bounded linear operators on B with B as domain of definition. A strongly continuous semi-group (SCSG) on B is a family {P t } t∈R ≥0 of bounded linear operators P t ∈ L(B) satisfying 1. P 0 = I B , the identity on B; 2. P t+s = P t P s for all s, t ∈ R ≥0 ; and 3. lim h 0 P h f -f = 0 holds for all f ∈ B. 12In the finite state case, each transition functions P t can be obtained from their q-matrix Q -simply by taking the matrix exponential P t = e tQ . In general, a SCSG is only the limit of a sequence of matrix exponentials of bounded operators, the so-called Yosida-approximants (see, e.g. Ref.

[EN00]). Part of the intricacy is explained by the fact that the time derivative d dt P t f at 0, i.e., the limits lim h 0 1 /h(P h f -f ), do not exists for all elements f of the Banach space. This brings us to a new object, which is typically related but in general not identical to the q-matrix (cf. [Reu57, Theorem 2]).

Definition 21 (Infinitesimal generator). The infinitesimal generator Q of a strongly continuous semigroup P t on a Banach space B is defined by

Qf = lim h 0 1 /h(P h f -f ) for all f ∈ dom(Q) = f ∈ B The limit lim h 0 1 /h(P h f -f ) exists. .
Infinitesimal generators determine the semigroup uniquely [EN00, Theorem 1.4]. To understand why, note first that it is sufficient to describe each P t on a dense subset of the Banach space; the domain of definition of the infinitesimal generator is simply the canonical such subset. Finally, for each f in the domain of the domain of the generator, the function t → P t f is characterised as the unique differentiable solution to the abstract Cauchy problem [Ein52] in Equation (5).

d dt u t = Qu t (t ≥ 0) (5) u 0 = f

E Drift functions for observables

Powers of length are drift functions (Lemma 2)

Let R ⊆ Σ + × Σ + be a finite string rewriting system and let n ∈ N + be a positive natural number. There exists a constant c n ∈ R >0 such that

|_| n : Σ + → R ≥0 is a c n -drift function.
Proof. Let K = max l r∈R |r|; it is a (rough) bound for the maximal length "growth" of rules in R. Define c n = (2K) n |R|. For all words v ∈ Σ + and all rules l r ∈ R,

(Q l r |n|)(v) = w∈Σ + q l r vw |w| n (by definition of Q l r |n|) = w∈Σ + \{v} q l r vw |w| n -   w∈Σ + \{v} q l r vw   q l r v |v| n (by conservativity) = w∈Σ + \{v} q l r vw (|w| n -|v| n ) ≤ w∈Σ + \{v} q l r vw ((|v| + K) n -|v| n ) (by definition of K) ≤ ((|v| + K) n -|v| n ) w∈Σ + \{v} q l r vw ≤ |v|((|v| + K) n -|v| n ) (≤ |v| redexes of l r in v) = |v| n k=1 n k |v| n-k K k (by Binomial Theorem and cancelling -|v| n ) ≤ |v| n k=1 n k |v| n-k K k ≤ |v| n k=1 n k |v| n-k K k |v| k-1 K n-k = |v| n k=1 n k |v| n-1 K n ≤ 2 n K n |v| n
The desired follows since Q R = l r∈R Q l r .

F Spieksma, corrigendum

The statement of Theorem 6.3 of Ref.

[Spi12] contains an error (corrected in our version) that can be spotted if one tries to fill in the details that are left out at the end of Spieksma's proof in op. cit. A counterexample is given in Appendix I.

The following proposition gives a useful and correct weakening of her statement, on a non-empty subset of the domain of the generator of the relevant SCSGs.

Proposition 2. Let P t be a transition function on S with q-matrix Q such that there exists a positive function W : S → R >0 and a constant θ such that QW ≤ θW , (i.e., (QW )(x) ≤ θW (x) for all x ∈ S). Let the family {P t } t∈R ≥0 be a SCSG on C 0 (S, W ), witnessed by V ∈ C 0 (S, W ) that is positive and a drift function for Q. For all f ∈ C 0 (S, W ) that satisfy f V < ∞ and Qf ∈ D(Q), we have

Qf = Qf = lim h 0 1 /h(P t f -f ), (6)
i.e., the latter limit exists in C 0 (S, W ), and in particular f ∈ D(Q).

G Computable analysis

We collect the main definitions and terminology from [BHW08], as required in [WZ07].

G.1 Naming schemes for sets

The central point of computable analysis in the tradition of the type-2 theory of effectivity (aka tte) are so-called representations of elements of sets (that later will always come equipped with a topology).

Definition 22 (Notations, representations, naming systems). Let X be a set, and let Σ be finite alphabet. A notation of X is a surjective partial function ν :⊆ Σ * → X, and, similarly, a representation is a surjective partial function δ :⊆ Σ ω → X; finally, a naming system for X is a notation or a representation of X. Given an element x ∈ X and a naming system of X, a γ-name is an element of γ -1 ({x}) (which thus is a finite or infinite string over Σ that is mapped to x).

G.2 Computable reals

There are several equivalent definitions of computable real numbers [BHW08, Theorem 3.2].

Definition 23 (Computable real [BHW08, Definition 3.1]). A real number

x ∈ R is computable if {q ∈ Q | q < x}
is a decidable set (of rational numbers).

Definition 24 (Computable sequence of reals [BHW08, Definition 3.8]).

A sequence (x n ) n∈N of real numbers is called computable if there exists a computable sequence (q k ) k∈N of rational numbers such that |x n -q n,i | < 2 -i for all n, i ∈ N.

Definition 25 ([BHW08, Definition 3.11.]).

1. If (r n ) n∈N is a convergent sequence of real numbers with limit x, and m :

N → N is a function such that i ≥ m(n) implies |x -r i | < 2 -n (for all i, n ∈ N),
we call the function m a modulus of convergence of the sequence (r n ) n∈N . 2. A sequence (r n ) n∈N of real numbers converges computably if it has a computable modulus of convergence.

G.3 Computable functions on infinite sets

A function f :⊆ R n → R is computable if there is an oracle Turing machine that, given any k ∈ N, may ask for arbitrarily good rational approximations of the input x ∈ dom(f );13 after finitely many steps, it writes a rational number q on the output tape with |f (x) -q| < 2 -k [BHW08, Definition 4.1.]. However, for the purposes of computability, instead of using oracle Turing machines, it suffices to consider a simpler model, namely Turing machines with infinite I/O, i.e., Turing machines equipped with one infinite input and one infinite output tape that are read only and one-way write only, respectively. This computation model allows for a natural notion of computability on infinite strings.

Definition 26 (Computable stream functions [BHW08, Definition 4.7]).

Let Σ be an alphabet. A function F :⊆ Σ ω → Σ ω is computable if there exists a Turing machine that, given a stream p ∈ dom(F ) on the input tape, writes the output stream F (p) on the one-way output tape, and that, given a stream p ∈ Σ ω dom(F ) does not write infinitely many symbols on the output tape (and instead moves to the halting state at some point).

If X is a space, a (computable) representation of X is a (computable) map δ X :⊆ Σ ω → X.

A particular representation that we shall invariably use for R with the usual (Euclidean metric) is ρ, given as follows: Let, for each n

∈ N, Q n = {m2 --n : m ∈ Z}, and let Q D = n Q n be the set of dyadic numbers. Let ν D : N → Q D be any computable bijective numbering such that if •, •, • : N 3 → N is a computable pairing function, we have ν D ( i, j, k ) = (i-j)2 -k . The representation ρ : {0, 1} * → R is now defined by dom(ρ) = {p ∈ {0, 1} * : ∀k.(ν D (p(k) ∈ Q k ∧|ν D (p(k))-ν D (p(k +1)| < -(k+1) )))} and ρ(p) = lim n ν D (p(n)
Note that ρ is merely a representation that uses rapidly converging Cauchy sequences and that the choice of {0, 1} * is immaterial: we could have chosen a larger alphabet, or even a unary alphabet if need be.

The above definitions carry over with ease to computation between arbitrary sets w.r.t. suitable representations.

Definition 27 (Computable function [BHW08, Definition 4.7]). Let Σ be an alphabet, and let

f :⊆ X → Y be a function; moreover, let δ X :⊆ Σ ω → X and δ Y :⊆ Σ ω → Y be representations. A function F :⊆ Σ ω → Σ ω is called a (δ X , δ Y )-realizer of f if δ Y F (p) = f δ X (p) holds for all p ∈ dom(f δ X ), i.e., such that for any δ X -name of some x ∈ dom(f ), the value F (p) is a δ Y -name of f (x). The function f is (δ X , δ Y )-computable if a computable (δ X , δ Y )-realizer of f exists.

G.4 Computable metric spaces

Definition 28 (Computable metric space [BHW08, Definition 7

.1]). A triple B = (X, d, α) is called a computable metric space, if 1. d : X × X → R is a metric on X, 2. α : N → X is a sequence such that the set {α(n) | n ∈ N} is dense in X, 3. d • (α × α) : N 2 → R is a computable map (with R equipped with the repre- sentation ρ).
Observe in the above definition that countable metric spaces are separable spaces and they come with a choice of a dense subset via an enumeration.

H Computability of duality pairing

Weighted C 0 -spaces C 0 (S, W ) and their duals L 1 (S, W ) with recursively enumerable S will be equipped with a recursive enumeration α of all finite rational linear combinations of indicator functions 1 x of states x ∈ S (observe that the set of such finite linear combinations is countable and can obviously be enumerated by an algorithm). As the indicator functions themselves are computable, and the map W → Q + is computable, it follows easily that both spaces are computable Banach spaces. Moreover, we equip the product space L 1 (S, W ) × C 0 (S, W ) with the sum of the norms of the underlying spaces, i.e., (π, f

) = π * W + f W where f W = sup x∈S |f (x)| /W (x) and π * W = x∈S W (x)|π(x)|.
Finally, we choose a recursive enumeration of all linear combinations of the set x∈S {(1 x , 0), (0, 1 x )} to make the product space a computable Banach space.

Lemma 3 (Computability of duality pairing).

Let S be a recursively enumerable set, and let W : S → Q + be a computable function. The dual pairing , : L 1 (S, W ) × C 0 (S, W ) → R, mapping (π, f ) to x∈S π(x)f (x), is computable.

Proof. For any two elements (π, f ), (π , f ) ∈ L 1 (S, W )×C 0 (S, W ) of the product space, if π is finitely supported, i.e., the set of states K = {x ∈ S | π (x) = 0} is finite, we have the following estimate for the distance of π , f and π, f in R.

Next, we describe an algorithm that for each desired precision = 2 -d as input, first, makes three successive requests for approximations (π , f ) of (π, f ) of ascending precision based on the estimate from (8) and Equation ( 9), followed by computation of a rational approximation r of π, f , before finally, giving the result r as output by writing up to d+1 digits of precision in binary representation.

1. Request (π 1 , f 1 ) as approximation of (π, f ) up to precision 1. Store F = 2 f 1 W + 1; now we have F > f W as a save over-approximation of the norm f W of f . 2. Request (π 2 , f 2 ) as approximation of (π, f ) up to precision min{ /3F , 1}). This ensures f 2 W π -π 2 * W < /3 (cf. ( 8)). Moreover, with ψ x = π(x) -π 2 (x),

x∈S |ψ x f 2 (x)| = x∈S W (x)|ψ x | |f 2 (x)|/W (x) ≤ f 2 W x∈S |W (x)ψ x | = f 2 W π -π 2 * W
< /3 [START_REF] Vincent Danos | Abstracting the differential semantics of rule-based models: Exact and automated model reduction[END_REF] 3. Compute the support of π 2 , i.e., K 2 = {x ∈ S | π 2 (x) = 0} using the representation of (π 2 , f 2 ). Compute M = max x∈K2 W -1 (x). 4. Request (π 3 , f 3 ) up to precision min{ /3M2 π2 * W , /3F , 1}. thus, the second approximation of π 2 to keep the support K 2 , and the best approximation of f 3 . 6. Return enough digits of r in the dyadic representation, d + 1 after the dyadic point.

To show that the algorithm is correct, we derive as in ( 8)

|r -π, f | = | π 2 , f 3 -π, f | ≤ x∈K |π 2 (x)f 3 (x) -π(x)f (x)| + f W π -π 2 * W
and then, as in 9, we obtain

≤ x∈K π 2 (x)f 3 (x)- π(x)(f 3 (x) -f (x)) + ψ x f 3 (x) + π 2 (x)f 3 (x) + f W π -π 2 * W = x∈K |-π(x)(f 3 (x) -f (x)) -ψ x f 3 (x)| + f W π -π 2 * W ≤ x∈K |π(x)(f 3 (x) -f (x))| + |ψ x f 3 (x)| + f W π -π 2 * W ≤ x∈K |π(x)(f 3 (x) -f (x))| + |ψ x f 3 (x)| + f W π -π 2 * W
and then as in ( 10) and ( 11)

≤ /3M2 π2 * W M 2 π 2 * W + f 3 W π -π 2 * W + f W π -π 2 * W ≤ /3 + f 3 W /3F + f W /3F < .

I Counterexample [Spi12, Theorem 6.3]

Example 4 (Counterexample to Theorem 6.3 of Ref. [Spi12]). Consider the positive natural numbers N \ {0} as state space of a CTMC with q-matrix Q on S = N \ {0} given by

q xy =      x 2 y = x -1 -x 2 y = x = 1 0 y = x = 1 (12)
and consider the function W (x) = x 2 , which is positive as x > 0. Now, W is a drift function since, for x > 1, we have

QW (x) = y∈S q xy W (y) = x y=x-1
q xy W (y) = q x(x-1) W (x -1) -W (x))

= q x(x-1) (x -1) 2 -x 2 )

< 0 ≤ W (x)
and QW (1) = 0 ≤ W (1). Thus, Q is a drift function.

Next, the identity function V (x) = x on S belongs to C 0 (S, W ) and is positive. Moreover, QV (x) = q x(x-1) V (x -1) -V x)

= q x(x-1) x -1 -x)

= -q x(x-1)

= -x 2

and QV (1) = 0, whence

QV (x) = -(x) 2 • 1 S\{1} (x).
This implies that V is also a drift function and hence the minimal transition function P t with q-matrix Q is a SCSG on C 0 (S, W ), applying the first part of Theorem 6.3 of Ref. [Spi12]. Finally, we have V V = 1 < ∞ an also QV W = -W W = 1 < ∞. Thus, by the second part of Theorem 6.3 of Ref.

[Spi12], we conclude that V belongs to the domain D(Q) of the generator Q. This means by definition that the limit

|

  π , f -π, f | = | x∈S π (x)f (x) -π(x)f (x)| (7) ≤ x∈S |π (x)f (x) -π(x)f (x)| = x∈K |π (x)f (x) -π(x)f (xof f W ) ≤ x∈K |π (x)f (x) -π(x)f (x)| + f W x∈S\K |W (x)(π(x) -π (x) e x = f (x) -f (x) and φ x = π(x) -π (x), π(x)f (x) = π(x)(f (x) + e x ) = π(x)e x + π(x)f (x) = π(x)e x + (φ x + π (x))f (x)= π(x)e x + φ x f (x) + π (x)f (x).

  x∈K2 |π(x)(f (x) -f 3 (x))| ≤ f -

lim h 0 P

 0 h V -V h exists in C 0 (S, W ), thus V ∈ D(Q), and QV ∈ C 0 (S, W ). Now QV = QV in analogy to Sect. 1.4 of Ref [And91] which leads to QV ∈ C 0 (S, W ), which is false: QV = -(•) 2 • 1 S\{1} and -(•) 2 • 1 S\{1} / ∈ C 0 (S, (•) 2 )!

Definition 5 (Standard transition function [And91, p. 5f.]). A transition function on S is a family {P t } t∈R ≥0 of S × S-matrices P t

	= (p t,xy ) x,y∈S with
	non-negative, real entries p t,xy such that
	1. lim t 0 p t,xx = 1 for all x ∈ S;
	2. lim t 0 p t,xy = 0 for all x, y ∈ S such that y = x;
	3. P t+s = P t P s = ( z∈S p s,xz p t,zy ) x,y∈S for all s, t ∈ R ≥0 ; and
	4.

Table 1 .

 1 Transition functions acting on Banach spaces: state of the art

	1 is called a Yosida approximant. 6 The Yosida approximants
	are the cornerstone of the generation theorems [EN00, Corollary 3.6] that allow
	one to pass from the generator Q to the corresponding SCSG. The constant θ
	also bounds the growth of the SCSG in norm, that is to say P t ≤ M e θt for
	some M . This should already make clear that Equation (3) is crucial to obtain
	error bounds for results on the computability of SCSGs. In fact it is the starting
	point of the proof of the main result on the computability of SCSGs [WZ07,
	Theorem 5.4.2,p. 521].
	It remains to see whether we can find a Banach space to build ACPs which
	accomodate interesting (specifically unbounded) observables.

By definition, Ex(f (Xt)) = ω∈Ω f (Xt(ω)) dpx for some suitable probability space (Ω, p), where px is the probability measure p on Ω conditioned by X0 = x. Performing a change of variables, Ex(f (Xt)) = S f dX * t (px), where X * t (px) is the image measure of px through Xt. Since pt,x-= X * t (px)(-), we deduce Ptf (x) = Ex(f (Xt)).

In fact, string rewriting is simply the restriction of stochastic graph transformation [HLM06] to directed, connected, acyclic, edge labelled graphs with in and out degree of all nodes bounded by one, i.e., to graphs consisting of a unique maximal path.

Even when Qf is defined, one has to check Qf = Qf , that is to say: 1 /h(P h f -f ) converges to Qf in the Banach space norm. But this will turn out to be easy compared to finding sufficient conditions for Qf to be defined.

In case f does not belong to the domain of the infinitesimal generator Q, (3) will generate a 'mild' solution that is a solution to the integral form of (acp) which might not be differentiable everywhere.

The indicator function 1x is defined as usual as 1x(y) = δxy.

Note that Anderson only considers standard transition functions [And91, p. 6].

This last condition, in presence of the first two, is equivalent to continuity of the function t → Pt, from R ≥0 to L(B), if we consider L(B) with the strong operator topology.

Thus, the oracle Turing machine may ask finitely many questions of the kind "Give me a vector p ∈ Q n of rational numbers with d(x, p) < 1/2 i ," where the exponent i may depend on the answers to the previous questions.


G Machado

G Chagnon

D Favier

Analysis of the isotropic models of the Mullins effect based on filled silicone rubber experimental results

Keywords: filled silicone rubber, Mullins effect, stress-softening, deformation criteria, hyperelasticity

The Mullins effect of rubber like material is classically defined as the stress softening during initial loading cycles. This effect is not accounted when the mechanical properties of material are modeled by a simple hyperelastic strain-energy function. In order to capture the stress softening it is necessary to define a set of supplementary variables as well a dissipation function, which evolves with the deformation history. In this paper we first describe experimental results that illustrate stress softening in particle-reinforced silicone rubber for uniaxial, planar and equibiaxial traction. The results allow to analyze the stress softening for the three different load cases. First, with respect to the choice of a stress-softening measure, the energy loss was evaluated by comparing the stored elastic energy for the first and the second loadings. The results point out that the virgin energy and the first invariant parameters are the best choice. Nevertheless, the maximum principal elongation, classically used in Mullins effect modeling, is not able to describe the different load cases. Furthermore, the ability of different class of models to describe filled silicone rubber was studied. The results show that models with a non-proportional and non-homothetical second load paths seem to be more efficient.

Introduction

Like many filled rubber-like materials, as well as soft biological tissues and other biomaterials, filled silicone rubber exhibit highly non-linear mechanical behavior. This non-linearity is commonly characterized by large strain levels under static conditions, a nonlinear stress-strain response and strain rate dependence. Finally, under cyclic loading conditions, hysteresis and stress-softening behaviors are observed on both filled and unfilled elastomers. For the last half century all these cyclic phenomena were evidenced and modelled to support the most widely engineering applications. However, their microscopic mechanisms (breaking of weak chains, breaking of links, desabsorption of chains, etc.) explanation remains non-unanimous.

Hysteresis is rather related to the dissipative nature of material, i.e., related with viscoelasticity [START_REF] Bergstrom | Large strain time-dependant behavior of filled elastomers[END_REF] and viscoplasticity [START_REF] Miehe | Superimposed finite elastic-viscoelastic-plastoelastic stress response withdamage in filled rubbery polymers. experiments, modelling and algorithmicimplementation[END_REF] behaviors. It is characterized as the difference between the loading and unloading paths during a cycle. On the other hand the stress-softening phenomenon, also called Mullins effect [START_REF] Mullins | Softening of rubber by deformation[END_REF], can be described as a stress-softening phenomenon of the material after a first loading. It can be idealized as an instantaneous and irreversible softening of the stress-strain curve, due to rearrangements in the microstructure of the material that occurs whenever the load increases beyond its prior all-time maximum value. At times when the load is less than a prior maximum, nonlinear elastic behavior prevails. One time the previous maximum stretch is reached the loading path turns up and follows the primary curve again up to a new maximum.

Through the years, the Mullins effect has experimentally been observed in different deformation states and numerous models have been proposed, but there is non-unanimous opinion on the choice of the variable to describe the evolution of the Mullins effect. Many different criteria have been chosen during the last fifty years, the most used remains the strain energy density or the maximum principal elongation. Within this context, under the assumption of isotropy, the aims of this work is to analyze which are the most suitable modelling for Mullins effect for different loading cases. In this way, a filled silicone rubber is studied under various mechanical loadings, where the Mullins effect becomes evident. Then, the determination of the best criterion to describe the evolution of stress softening and the kind of model able to describe it for every loading cases are evaluated. Hence, the different type of modelling are summarized in 2. Material preparation procedures are described in Section 3, testing procedures and strain field measurements techniques based on optical methods are described in Section 4. Experimental results are then presented, in Section 5. An analyze of the strain measure criteria and the best way to model Mullins effect is finally presented in Section 6.

Mullins effect modelling

There exist a large number of Mullins effect constitutive equations, with different kind of isotropic or anisotropic modelling hypothesis. More and more modelling tackles the anisotropy of the Mullins effect (see for example [START_REF] Ehret | Modeling of anisotropic softening phenomena: Aplication to soft biolgical tissues[END_REF], [START_REF] Diani | Observation and modeling of the anisotropic visco-hyperelastic behaviorof a rubberlike material[END_REF] and [START_REF] Miehe | A micro-macro approach to rubber-like materials. part iii: The micro-sphere model anisotropic mullins-type damage[END_REF]). In this paper, we focus on isotropic approaches of the Mullins effect to analyze what are the best ways of modelling. Whatever is the modelling, a deformation measure parameter is used to control Mullins effect evolution and its choice is critical.

Choice of the modelling

A classical approach to describe Mullins effect is to use damage mechanic theory. The strain energy density of a virgin material W 0 (F) submitted to the deformation gradient F can be modified incorporating continuum damage parameter d in order to take into account the stress softening. The classical form is given by

W(F) = (1 -d)W 0 (F).
(1)

All proposed models vary from each other by the definition of the damage criteria and the form of the damage evolution function. The definition can involve any microstructural chain damage and microvoid formations both associated with physical phenomena like breaking of chains and/or links, desabsorption of chains, etc. A common approach assumes that damage is a function of a discontinuous quantity like the maximum applied stretch or the maximal energy density with respect to the primary load path. Nevertheless, as reported by [START_REF] Chagnon | On the relevance of continuum damage mechanics as applied to the mullinseffect: theory, experiments and numerical implementation[END_REF], these models assume that the stress ratio between two different secondary load curves are related by a constant multiplicative factor:

P j s P i s = 1 -d(x j s ) 1 -d(x i s ) (2)
where s denotes the secondary load curves numbered i = 1, ..., 3 and j = 2, ..., 4 (i < j) for the presented experimental data. Pseudo elasticity models are based on a different approach. [START_REF] Ogden | A pseudo-elastic model for the mullins effect in filled rubber[END_REF] present an additive function to describe the second loading curves, this permits to have two distinguished functions to describe the first and the second loading curves. In this case some conditions must be imposed to ensure the continuity of the strain-stress behavior. Now, no more proportionality is imposed between the second loading curves and the evolution function for the second loading curves is described thanks to W 0 measure. This permits to have a model able to characterize the different loading cases. This approach has been very employed in the past few years (see for example [START_REF] Dorfmann | A constitutive model for the mullins effect with permanent set in particulereinforcedrubber[END_REF] and [START_REF] Horgan | A theory of stress softening of elastomers based on finite chain extensibility[END_REF]).

The double network theory is different from damage mechanics as the multiplicative function used for describing Mullins effect is not constant on the second loading curves. It can no longer be considered as a 2

Table 1: Definition of damage criteria for some representative models Variable Description Reference 2W 0 (F) [START_REF] Simo | On a fully three-dimensional finite-strain viscoelastic damage model : formulationand computational aspects[END_REF], [START_REF] Miehe | Superimposed finite elastic-viscoelastic-plastoelastic stress response withdamage in filled rubbery polymers. experiments, modelling and algorithmicimplementation[END_REF], [START_REF] Septanika | Application of the network alteration theory for modeling the time dependent constitutive behavior of rubber. Part I -general theory[END_REF] W 0 (F) [START_REF] Ogden | A pseudo-elastic model for the mullins effect in filled rubber[END_REF], [START_REF] Laraba-Abbes | A new taylor-made methodology for the mechanical behaviour analysisof rubber like materials: II. application of the hyperelastic behaviourcharacterization of a carbon-black filled natural rubber vulcanizate[END_REF], [START_REF] Miehe | Discontinuous and continuous damage evolution in ogden type large strainelastic materials[END_REF] λ max Maximum largest principal stretch [START_REF] Laiarinandrasana | Visco-hyperelastic model with internal state variable coupled with discontinuousdamage concept under total lagrangian formulation[END_REF]) [START_REF] Chagnon | On the relevance of continuum damage mechanics as applied to the mullinseffect: theory, experiments and numerical implementation[END_REF], [START_REF] Qi | Constitutive model for stretch-induced softening of the stress-stretch behaviorof elastomeric materials[END_REF] . damage function as the value of the function is evolving on second loading curves and it comes back to a value of one on the first loading curve whereas the material has already be stretched. DeSouza [START_REF] Neto | A phenomenological three dimensional rate independent continuum damage modelfor highly filled polymers : formulation and computational aspects[END_REF], based on the idea of [START_REF] Gurtin | Simple rate-independent model for damage[END_REF] proposed a damage master curve, that describes the evolution of the network, that is magnified to the considered maximum previous strain level.

I 1 /3 -1 I 1 = trace(C) = trace(B)
ε: ε ε = 1 2 C -1 -I Lion (
A physical approach, proposed by [START_REF] Marckmann | A theory of network alteration for the mullins effect[END_REF] is based on the position of the strain hardening when the material reaches it maximum previous deformation level. They assume that the macro-molecular chain network is evolving with the maximum deformation, they assume that the number of active monomers is constant. This hypothesis was relaxed by [START_REF] Chagnon | Development of new constitutive equations for mullins effect in rubber usingthe network alteration theory[END_REF] who considered that some chains can become hanging, that means that the number of active monomers is a decreasing function of the deformation. In another way, [START_REF] Qi | Constitutive model for stretch-induced softening of the stress-stretch behaviorof elastomeric materials[END_REF] developed a double network model also based on evolution of chain elongation deformation measure. This permit to have independent second loading curves.

As shown, there exists a large kind of models (all existent models have not been presented here) but it seems that the choice of an energy deformation measure is quite important. This choice ensures an ability of the model to describe the different load cases, moreover the different modelling have different consequences on the strain-stress curves.

The history parameter choice

The choice of the Mullins effect evolution criterion is fundamental for the description of different loading cases. As pointed out by [START_REF] Diani | A review on the mullins effect[END_REF], many criteria are exhibited by different authors. The Table 1 summarizes different choices of the past years. Usually, the choice of the measure is not justified, but is able to represent the experimental curves. For an uniaxial representation, the choice of the parameter is not essential, as all can describe the stress softening. The Mullins effect evolution function is then built in function of this parameter. Nevertheless, in a tridimensional representation, the choice becomes very crucial as it rules the position of the second loading curves in the different deformation modes. In the next sections the experimental procedures and results are provided to enhance this discussion.

Preparation of the Silicone Specimens

The samples used in the experiments were made of a filled silicone rubber called Rhodorsil RTV3428 supplied as two liquid components: the uncured silicone and the curing agent. The final material is produced by a polyaddition, curing at room temperature. The liquid mixture is molded by injection to obtain a sheet with constant thickness. Although the experiments are performed separately and the strain states are different, data from all of the individual experiments is often used as a set. This means that the specimens used for each of the experiments must be made using the same protocol of elaboration, in order to obtain specimens with reproducible mechanical properties, as follow: (i) mixing the components with a 10/1 mass ratio, (ii) putting the mixture under vacuum for 20 minutes in order to eliminate undesirable entrapped air bubbles, (iii) sheet molding injection, (iv) putting mold inside an oven at 70 o C for 4 hours in order to accelerate the curing process and assure a sufficient cross-linking density; (v) sheet demolding after 1 hour exposed at ambient temperature. Finally, the external surface of the molded plate was coated with a stochastic silicone paint pattern [START_REF] Meunier | Mechanical experimental characterisation and numerical modelling of an unfilledsilicone rubber[END_REF]. Made of small speckles, the pattern is necessary for the digital image correlation (DIC) and stereo digital image correlation (SDIC) field measurements. Note that the quality of the coated pattern (size, density and grey contrast level) is a critical point to obtain a good estimation of the strain field measurement.

Testing procedures and strain field measurements techniques

In-plane tests

In plane quasi-static experiments were conducted on a MTS 4M universal testing machine with an Entran ELPM-T2±250N load cell. The images were recorded at 0.5 Hz with a Jai TM-4200GE CCD camera using a reduced scan of 2048 × 1000 pixels. were performed on rectangular samples having an initial gauge length l 0 = 40 mm, width w 0 = 13 mm and a thickness e 0 = 2 mm. Since the experiment was not intended to fail the specimen, there was no need to use a dumbbell shaped specimen commonly used to prevent specimen failure nearby the clamps. But knowing that the effects of clamps create an indeterminate state of stress and strain in the region surrounding the clamp, due to the process of gripping, the initial gauge length was adopted as being less than the real physical size of the samples. Figure 2(a) presents the tensile test specimen. Second, the pure shear strain state was approached by performing planar shear tensile test (see Figure 2(b)). The initial height l 0 , the constant width w 0 and the thickness e 0 of samples were 8 mm, 70 mm and 2 mm, respectively. These dimensions have as objective to create an experiment where the specimen is perfectly constrained in the lateral direction such that all specimen thinning occurs in the thickness direction.

The DIC technique, like a non-contact method, is often used to characterize rubber-like materials, see for example [START_REF] Meunier | Mechanical experimental characterisation and numerical modelling of an unfilledsilicone rubber[END_REF] and [START_REF] Sasso | Characterization of hyperelastic rubber-like materials by biaxial and uniaxialstretching tests based on optical methods[END_REF]. In each testing time step an image of the deformation of samples was recorded using a CCD camera. Using the DIC, it is possible to reconstruct the surface of the deformed samples and determinate the full-field surface displacements. An accuracy on the order of 10 -2 pixels or better for in-plane displacement components and point-to-point strain accuracy of ±10 -4 in-plane can be obtained. The DIC method is preferred because there is no change in stiffness resulting from the presence of a attached sensor (classical extensiometry) and for being insensitive to ambient vibrations and rigid body motions (speckle pattern interferometry). Moreover the method is able to deal with high strain levels, what is a very desirable feature given the high strain level experienced by the tested material. Also the DIC allows demonstrating and measuring an heterogeneous deformation field. See [START_REF] Sutton | Digital image correlation for shape and deformation measurements[END_REF] for further explanation about digital image correlation method.

Out-plane tests

A bulge test was also conducted in order to determinate a equibiaxial state. This test, also called "balloon" test, consists of a thin disk specimen, of initial diameter d 0 = 200 mm with e 0 = 2 mm of thickness, constrained between two clamping flanges. This test machine, performed at the same way in [START_REF] Meunier | Mechanical experimental characterisation and numerical modelling of an unfilledsilicone rubber[END_REF], consists of a air piston connected to the bottom circumferential clamp. The internal pressure is measured with a pressure sensor Foxboro 0-500 mbar. For these tests, out-of-plane displacements measurements were needed, therefore a Dantec Q-400 Stereo Digital Image Correlation (SDIC) system was used to acquire and correlate the images. Figure 3 presents the experimental setup. Since a large out plane displacement is experienced by the silicone sample, a good focal distance is absolutely necessary. Nevertheless, a high spatial resolution is needed to permit an accurate measurement. Using the stereocorrelation it is possible to reconstruct the deformed sample surface geometry as well as determinate the displacement and strain fields. Because rigid body motion has no effect on the measurements, a high level of strain has been achieved in this kind of experiment. The Figure 4 shows the strain level curves for the maximal pressure value of 18 KPa in the center of the membrane. The curvature on the top of the bulge is computed, using a post-treatment routine.

Experimental results

Different tests are realized to evaluate the Mullins effect in different loading cases. First, time dependence of the material is studied and then the stress-softening is analyzed.

Time dependent response

A creep tensile test was performed, where a constant monotonic load of 20 N was applied resulting in a initial stretch of λ i = 2.160. It can be observed in Figure 5 that the material response tends towards an equilibrium state, which cannot be reached within laboratory time scale. However the stretch increased by less than 1% within the observed time range. In a second time, a tensile specimen has been subjected to a load/unload sequence at different elongation rates λ = 2.5 × 10 -3 s -1 to 1.25 s -1 in order to verify the rate-dependence influence. No noticeable difference between stress-strain responses was observed at the considered strain rate range. Consequently the RTV3428 behavior can be assumed independent of the rate of deformation for the observed ranges.

Uniaxial tension test

During the test, using an elongation rate of λ = 1.60×10 -2 s -1 , the nominal stress tensor P (First Piola-Kirchhoff stress tensor) is assumed to be homogeneous within the gauge region as well as the deformation gradient tensor F. Since the current thickness is not measured, the material is assumed to be incompressible, i.e., det (F) = 1 for convenience. In the central zone, the deformation gradient, considering that the direction 1 is the loading one, is given by:

F = λ(e 1 ⊗ E 1 ) + λ -1/2 (e 2 ⊗ E 2 + e 3 ⊗ E 3 ) (3)
and the nominal stress tensor is:

P = P 11 (e 1 ⊗ E 1 ) (4)
where (E 1 , E 2 , E 3 ) and (e 1 , e 2 , e 3 ) are the initial and deformed orthonormal basis respectively, here identical. A cyclic loading test was realized, the results are presented in Figure 6. Different phenomena are highlighted, first a large Mullins effect appears by comparing the two first loading at each strain level, but with very few residual elongation. Moreover, a little hysteresis after the first cycle (difference between loading and unloading during the second cycle) is observed.

Plane strain tensile test

Also, the strain state during the test was supposed to be homogeneous, by assuming the direction 1 is the loading one and direction 2 is the breadth of the specimen direction, the deformation gradient tensor was expressed as:

F = λ (e 1 ⊗ E 1 ) + (e 2 ⊗ E 2 ) + λ -1 (e 3 ⊗ E 3 ) (5)
where the stress tensor is

P = P 11 (e 1 ⊗ E 1 ) + P 22 (e 2 ⊗ E 2 ) (6)
A cyclic planar loading test was realized, the results are presented in Figure 7. Planar tensile response, likewise uniaxial traction, presents the same phenomena. Here, the residual stretch shows to be more evident. The maximum principal stretches experienced by the planar specimens are smaller if compared with uniaxial traction specimens. This limitation lies in the fact that the planar traction specimens must to be constrained in the lateral direction without slipping.

Equibiaxial tension state approached by the bulge test

Due to the axial-symmetry of the experimental configuration the equibiaxiality of the stress and strain is obtained on the top the inflated sample. The deformation gradient tensor in a local coordinate system is

F = λ t (e 1 ⊗ E 1 + e 2 ⊗ E 2 ) + λ -2 t (e 3 ⊗ E 3 ) (7)
where λ t is the elongation in the tangential direction of the specimen measured by SDIC. Finally, with all assumptions above, the stress tensor is given by Once thickness dimension is much smaller than the other two planar dimensions (d 0 e 0 ), the stress is also assumed to be negligible along the thickness. Another consideration comes from the fact that the curvature of the inflated sample was presumed to be the same along all directions at the disk center.The component P 11 can be calculated from knowledge of curvature radius r, the tangential elongation λ t and the pressure p recorded during the test:

P = P 11 (e 1 ⊗ E 1 + e 2 ⊗ E 2 ) + P 33 (e 3 ⊗ E 3 ) . ( 8
P 11 = P 22 = pr 2e 0 λ t ; P 33 ≈ 0. ( 9
)
The results for the central area are presented in Figure 8 for a cyclic loading. The response are qualitatively similar to uniaxial loadings with hysteresis and Mullins effect.

Discussion

On the one hand, by use of experimental data, we propose to analyze the consequence of the choice of a stress-softening measure. Like, often described, by considering that the main part of the Mullins effect is represented by the difference between the first and the second loading paths, the energy loss is evaluated by comparing the stored elastic energy for the first and the second loadings. The evolution of the ratio of these two quantities is plotted according to all the Mullins effect measure parameters of Table 1, for the three loading cases: uniaxial and equibiaxial extensions and planar tension. The results are presented in Figure 9. It can be noticed that none of the parameters present a perfect correlation for all loading cases. Although, for the tested strain range, it seems that the virgin energy (W 0 ) and the first invariant (I 1 ) parameters are the best choice. The maximal principal stress (P max ) can be acceptable for low strains levels. What clearly appears is that the maximum principal elongation (λ max ), usually used in Mullins effect modelling, is not able to describe the different loading cases. All the other parameters derived from the trace operator of a strain measure are not more efficient. All these measures become less and less efficient at large deformation. It emphasizes the fact that the Mullins effect should be described by W 0 or I 1 parameter. All this means that the use of a non-efficient parameter implies the construction of the second loading modelling on wrong bases. Whatever is the model, the evolution function used to describe the stress-softening will be fitted on a non-adequate parameter. Many Mullins effect modelling have a large number of parameters, this permits to have a good fit, but if the wrong deformation criterion is used, the physics of the phenomenon is lost.

On the other hand, some modelling hypothesis are discussed. A damage modelling approach imposes to these ratios to be constant whatever is the maximum deformation parameter. Figure 10 presents the ratios between the different second loading curves (plotted according to the principal elongation). The beginning of the curves varied very quickly, but this is due to experimental errors which are more significative at low strain levels. It clearly appears that no damage model can correctly represent the Mullins effect of a material that has a similar behavior than the RTV3428 silicone. A choice must be done, but the beginning of second loading curves and strain hardening cannot be described correctly at the same time using these models. The proportionality between the second loading curves prevents from describing large strain hardening, then it can only be applied to material that have a weak influence of Mullins effect. But this discards the imposition of proportionality between second loading curves, like proposed by [START_REF] Beatty | A theory of stress-softening in incompressible isotropic materials[END_REF] and [START_REF] Zuñiga | A new phenomenological model for stress-softening in elastomers[END_REF]. But their proposition is written thanks to a strain measure that does not seem to be very efficient to characterize all loading cases. This was improved by [START_REF] Zuñiga | A phenomenological energy-based model to characterize stress-softening effectin elastomers[END_REF] by changing the strain measure to an energy measure. An interesting analysis and criticism of this model is proposed by [START_REF] Kazakeviciute-Makovska | Experimentally determined properties of softening functions in pseudo-elasticmodels of the mullins effect[END_REF].

DeSouza [START_REF] Neto | A phenomenological three dimensional rate independent continuum damage modelfor highly filled polymers : formulation and computational aspects[END_REF] proposed a master second loading curve to describe Mullins effect, the curve is magnified to join the stress-free state to the maximum deformation. Figure 11 presents a normalization of the strain-stress curve for the three simple tests for the maximum strain and the normalized strain energy. It clearly appears that all the curves are not superimposed, that means that the idea of a single master curve is not adequate. The form of the second loading curve is varying a lot, according to the maximum deformation. The strain hardening is very different and is not appearing at the same in a second loading curve. This phenomenon is also more important for Natural Rubber because of strain induced crystallization.

All these remarks emphasize the fact that an efficient Mullins effect modelling, by an isotropic approach, must not use a formalism where the second loading curves are proportional or homothetical. The modelling should be able to describe different form of strain-hardening according to the maximum previous deformation level. Moreover governing Mullins effect parameter should be built thanks to different loading cases; a simple uniaxial test can drive to abnormal conclusion. We suggest to use the elastic energy or the first strain invariant. The physical or pseudo-elastic approaches seems to be the best way to tackle correctly Mullins effect evolution functions to be able to well describe all second loading curves.

Conclusion

In the present paper, limitation of isotropic models of the Mullins effect have been tackled thanks to experimental tests. In this way, filled silicone rubber has been characterized. The reproducibility and accuracy of measurements were evaluated through several successive tests. Results have shown that our protocol to manufacture the silicone specimens and our experimental methodology work well and provide a precise experimental characterization of softening phenomena. Uniaxial, equibiaxial and planar tensions have been realized, this permits to compare first and second loading curves in different strain modes. Energy losses are compared and analyzed. It clearly appears that the strain energy measure (W) is the best choice and in a lesser the first strain invariant (I 1 ).

The analysis of experimental data proved that proportionality between second loading curves is not a good way to model the strain-stress curves. There are two main ways to well describe the Mullins effect. The first one, by using a non-constant multiplicative function on second loading curves or by using different non-proportional functions, this respectively corresponds to double-network and physical approaches.

The governing parameters and description function forms have been analyzed under isotropic assumption. The main extension of this study is to, now, take into account the influence of the loading direction of the stress-softening.

Acknowledge

We would like to thank the French ANR for supporting this work through the project RAAMO ("Robot Anguille Autonome pour Milieux Opaques")

  , ε and e are respectively the Right Cauchy-Green,Left Cauchy-Green, Piola and Euler-Almansi strain tensors. I is the metric tensor.

  Figure 1 presents the experimental setup. First, simple tensile tests

Figure 1 :

 1 Figure 1: In-plane experimental setup. (a) silicone specimen; (b) clamps; (c) load cell; (d) omni light; and (e) CCD camera.

Figure 2 :

 2 Figure 2: Test specimen dimensions (a) uniaxial tensile specimen (b) planar tensile specimen.

Figure 3 :

 3 Figure 3: Bulging test set-up. (a) bulge specimen; (b) circular clamp; (c) CCD cameras; and (d) pressure sensor.

Figure 4 :Figure 5 :

 45 Figure 4: First principal strain (Green-Lagrange measure) in the center of the membrane for the maximal pressure value of 18 KPa.

Figure 6 :

 6 Figure6: Cyclic loading-unloading tensile test with increasing maximum stress: 0.2, 0.4, 0.6 and 0.8 MPa at λ = 1.60×10 -2 s -1 .

Figure 7 :

 7 Figure 7: Cyclic loading-unloading planar tensile test. Maximum stress: 0.18, 0.28, and 0.38 MPa at λ = 2.60 × 10 -2 s -1 .

Figure 8 :

 8 Figure8: Cyclic equibiaxial stress-strain behavior obtained using the tension bulge test.

Figure 9 :

 9 Figure 9: Energy loss between the primary path and the second load for the different mechanical loadings with respect to (a) the first strain invariant I 1 ; (b) the second strain invariant I 2 ; (c) maximal principal stretch λmax; (d) maximal principal stress Pmax; (e) strain energy density W 0 ; (f) Left Cauchy-Green strain measure; (g) Piola strain measure; (h) Euler-Almansi strain measure.

Figure 10 :

 10 Figure 10: Second loads ratio. Rows correspond, respectively, to uniaxial planar and equibiaxial traction. In the left column dot lines represent the primary path and continuous lines are the second loads. The right column shows the ratio between second loads.

Figure 11 :

 11 Figure 11: Normalized second stress loading curves with respect the maximum strain (left column) and the normalized energy of deformation (right column). Rows correspond, respectively, to uniaxial planar and equibiaxial traction.