Landy Rabehasaina

Jae-Kyung Woo

Multitype branching process with nonhomogeneous Poisson and contagious Poisson immigration

Keywords: AMS 2000 subject classications: Primary 60J80, 60J85, secondary 60K10, 60K25, 90B15 Multitype branching process with immigration, Nonhomogeneous Poisson process, Contagious Poisson process, Convergence in distribution

In a multitype branching process, it is assumed that immigrants arrive according to a nonhomogeneous Poisson or a contagious Poisson process (both processes are formulated as a nonhomogeneous birth process with an appropriate choice of transition intensities). We show that the normalized numbers of objects of the various types alive at time t for supercritical, critical, and subcritical cases jointly converge in distribution under those two dierent arrival processes. Furthermore, some transient expectation results when there are only two types of particles are provided.

Introduction

Single or multitype branching processes with dierent stochastic assumptions on the immigration process have been applied in diverse elds in applied probability such as biology, epidemiology, and demography. For example, [START_REF] Resing | Polling systems and multitype branching processes[END_REF] utilized the theory of multitype branching processes in discrete time with immigration to study the joint queue length process in the dierent queues of a polling system in queueing theory. More generally, a network of innite server queues may be seen as multitype Galton Watson processes with immigration, see e.g. [START_REF] Altman | On stochastic recursive equations and innite server queues[END_REF]. Some actuarial applications of branching processes such as a reinsurance chain were discussed in [START_REF] Rolski | Stochastic Processes for Insurance and Finance[END_REF]Section 7.5]. Regarding applications in biology, the recent paper [START_REF] Hyrien | Stochastic modeling of stress erythropoiesis using a two-type age-dependent branching process with immigration[END_REF] considered multitype branching processes with homogeneous Poisson immigration to study stress erythropoiesis, although the authors pointed out that a nonhomogeneous Poisson process (NHPP) might be more realistic in that situation. The reader is referred to [START_REF] Mitov | Multitype branching processes with inhomogeneous Poisson immigration[END_REF] for a detailed discussion about the relevant literature on various types of branching processes.

In this paper, we consider a multitype branching process in which there are dierent types of particles, and new particles that arrive according to a NHPP or a contagious Poisson process (CPP). For immigration processes, an alternative to homogeneous Poisson process, NHPP and CPP are chosen for the following reasons. NHPP and CPP are within the class of nonhomogeneous birth processes, which means the intensity of event occurrence possibly varies with time (e.g. seasonality of catastrophe incidence) (NHPP) or the past state of the process (e.g. the number of previous shocks or the number of accidents incurred in the past) (CPP).

As a result, these non-Poissonian processes have been widely used in various areas such as engineering, applied probability, biological science, and actuarial science. In particular, a CPP is an immigration process with a linear birth rate and is basically a homogeneous case of generalized Polya processes (GPP) studied in [START_REF] Willmot | Distributional analysis of a generalization of the Polya process[END_REF][START_REF] Cha | Characterization of the generalized Pólya process and its applications[END_REF][START_REF] Cha | New shock models based on the generalized Polya process[END_REF] for example. This process is also regarded as a contagious model which is suitable to describe the spread of an infection as proposed in [START_REF] Pólya | Sur quelques points de la théorie des probabilités[END_REF] and the accident proneness as considered by [START_REF] Greenwood | An inquiry into the nature of frequency distributions representative of multiple happenings with particular reference to the occurrence of multiple attacks of disease or of repeated accidents[END_REF], [START_REF] Bates | Joint distributions of time intervals for the occurrence of successive accidents in a generalized Polya scheme[END_REF], and [START_REF] Bühlmann | Mathematical methods in risk theory[END_REF]Section 2.2].

Since a branching process can be used to study the dynamic network of the spread of infectious diseases, a CPP is naturally a suitable choice for the immigration arrival process to model the occurrence of contagious events due to the explanation given previously. In addition, we refer to some papers which considered GPP in other elds. For example, see [START_REF] Konno | On the exact solution of a generalized Polya process[END_REF] for a more detailed discussion on the use of GPP in the framework of a non-stationary type master equation approach in mathematical physics. As a GPP models the decreasing sequence of interarrival times depending on the number of the events in the past (of which their positive dependence is proved in [START_REF] Cha | Justifying the Gompertz curve of mortality via the generalized Polya process of shocks[END_REF]), it becomes a more realistic choice in modeling mortality as well as the increasingly damaging impact on an aging system as studied in [START_REF] Cha | Justifying the Gompertz curve of mortality via the generalized Polya process of shocks[END_REF].

The focus of our paper is to study the joint asymptotic behavior of a process representing the numbers of dierent types of particles alive at time t when the immigration process is described by a NHPP or a CPP. These models may be interpreted dierently depending on whether we are in an epidemic, actuarial, queueing or reliability setting. In an epidemic setting, the particles represent contaminated cells, the types represent their locations, and those cells can move to other locations where they can possibly contaminate other cells. In an actuarial setting, a particle may represent a certain type of claims or tasks that needs to be processed in dierent branches of an insurance company before being settled, or are in dierent stages of a reinsurance contract as explained in [START_REF] Rolski | Stochastic Processes for Insurance and Finance[END_REF]Section 7.5]. In a queueing setting, a particle is a customer who arrives and gets served immediately in the setting of innite server queues and after leaving the queue, is replicated into several new customers who are sent to other queues for their subsequent service. In a reliability setting, particles are interconnected parts in a system which can be damaged upon external shock arrivals and need to be repaired or are dependent line outages in a power network which may cause cascading blackouts, see [START_REF] Qi | Estimating the propagation of interdependent cascading outages with multi-type branching processes[END_REF].

Most papers in the literature consider the critical case of underlying branching mechanisms, see for e.g. [START_REF] Weiner | On a multi-type critical age-dependent branching process[END_REF][START_REF] Durham | A problem concerning generalized age-dependent branching processes with immigration[END_REF][START_REF] Weiner | A multi-type critical age-dependent branching process with immigration[END_REF][START_REF] Mitov | Multitype branching processes with inhomogeneous Poisson immigration[END_REF]. In this paper, we consider all three underlying branching mechanisms, including supercritical, critical and subcritical cases. Indeed, it is well known that in the subcritical and critical cases for a continuous time multitype Galton-Watson process, i.e. when the eigenvalue of the mean matrix of ospring does not exceed 1, extinction is certain, whereas the survival probability in an innite horizon is positive in the supercritical case.

Depending on the modeling context, either a critical, subcritical or supercritical case may be more appropriate for modeling the dynamics of the quantity of interest. In particular, in a queueing or actuarial context where the clients or tasks will eventually exit the system, it may be more plausible to use a critical or subcritical case. For example, in the case of polling systems, the stable case corresponds to the subcritical branching process and the heavy trac limit is studied using the near critical branching process in [START_REF] Van Der Mei | Towards a unifying theory on branching-type polling systems in heavy trac[END_REF]. On the other hand, in the context of epidemiology, a supercritical case may be more appropriate to model the rapid expansion of a particular disease in the beginning of the outbreak.

The three cases of underlying branching mechanisms exhibit dierent behaviors of the branching process with immigration. Concerning the immigration arrival process, particular attention in the forthcoming results is given to the case where the immigration rate increases very fast. This is a property of the GPP, as the (stochastic) arrival rate is linear with respect to the number of arrivals at the current time, so that the interarrival time decreases with respect to the failure rate order, as explained in [START_REF] Cha | Justifying the Gompertz curve of mortality via the generalized Polya process of shocks[END_REF][START_REF] Badía | Univariate and multivariate stochastic comparisons and ageing properties of the Generalized Pólya process[END_REF]. In the NHPP, assumptions on the intensity function are such that an exponential behaviour for the latter is studied, yielding dierent kinds of asymptotic results. For a single type branching process with general lifetime distribution of the particles, the reader is referred to [START_REF] Hyrien | Supercritical Sevastyanov branching processes with non-homogeneous Poisson immigration[END_REF] and [START_REF] Hyrien | Subcritical Sevastyanov branching processes with nonhomogeneous Poisson immigration[END_REF] which proved asymptotic results and functional central limit theorems respectively in the supercritical and subcritical cases.

In this paper, we also investigate the transient moment of the process for two types of particles branching mechanism when the renewal function associated with the arrival process is explicit. See e.g. [START_REF] Hyrien | Stochastic modeling of stress erythropoiesis using a two-type age-dependent branching process with immigration[END_REF] for a similar study when the particular lifetime particles have a general distribution.

The rest of the paper is organized as follows. In Section 2, a multitype branching process without/with immigration and relevant assumptions are described. Some convergence results for the distribution of a number of dierent types of particles at time t (denoted in vector form as N (t)) under dierent immigration processes (NHPP, CPP) are obtained from its characteristic function (CF). In Section 3, NHPP is assumed for the arrival process of immigrants and convergence results for N (t) are given in Theorems 3 and 4, with a particular emphasis in the case when the intensity of the arrival process increases exponentially. A result in the critical case is given in Theorem 5, which agrees with the previous results in the same context in [START_REF] Weiner | A multi-type critical age-dependent branching process with immigration[END_REF][START_REF] Vatutin | A critical Bellman-Harris branching process with immigration and several types of particles[END_REF][START_REF] Mitov | Multitype branching processes with inhomogeneous Poisson immigration[END_REF]. For the critical case, some remarks for homogeneous Poisson immigration and one dimensional branching process with immigration are provided in Remark 6. In subsections 3.1, 3.2, and 3.3, detailed proofs of Theorems 3, 4 and 5 are given. Section 4 considers a CPP for the immigration process. Asymptotic behaviors of N (t) are studied in Theorem 10 in the function of the parameters of the arrival process with detailed proofs included in subsections 4.1, 4.2 and 4.3. In the proofs of Theorems 3, 4, 5 and 10, we show that, for a conveniently chosen normalizing function g(t), the process N (t)/g(t) converges in distribution to an identiable limit as t → +∞ by proving that the corresponding CF converges. Finally, some transient results in scenarios when there are two types of particles in the branching process are presented in Section 5.

Lastly, the following matrix notations will be used throughout the paper. For any matrix

M ∈ R m×n , M ∈ R n×m denotes its transpose. By < u, v >= k i=1 u i v i we denote the usual inner product between two vectors u = (u 1 , ..., u k ) and v = (v 1 , ..., v k ) , associated with the corresponding Euclidean norm ||u|| = k i=1 u 2
i . We let 1 = (1, ..., 1) , a vector with 1's of an appropriate dimension, R k + = [0, +∞) k and R * k + = (0, +∞) k . We also let L 1 (R + ) be the set of integrable measurable functions from R + to R. Finally, we denote N * = {1, 2, ...} and N = N * ∪ {0}.

The model

The baseline model, a classical multitype branching process (without immigration), is described as follows. We consider a set of particles of k possible types, with a type i particle having exponential lifetime with mean 1/µ i for i = 1, ..., k, denoted by E(µ i ) for µ i > 0. Upon its death, a type i particle produces Y (i) j copies of type j particles for all j = 1, ..., k, where (Y

(i) 1 , ..., Y (i) k ) is a random vector with corresponding probabilities p i (n) = p i (n 1 , ..., n k ) = P(Y (i) j = n j , j = 1, ..., k) for n = (n 1 , ..., n k ) ∈ N k ,
and generating functions dened by

h i (z) = h i (z 1 , ..., z k ) = n∈N k p i (n)z n 1 1 • • • z n k k , i = 1, ..., k, z ∈ [0, 1] k . (1)
In other words, p i (n) is the probability that type i particle produces n 1 , ..., n k copies of type 1, ..., k particles respectively. Then all copies evolve independently and have the same dynamics. Note that p j (0, ..., 0) is the probability that no replica is made, i.e. the probability that the particle does not produce any copies at the end of its lifetime. The mean numbers of copies from type i particle are denoted by

(m i,1 , ..., m i,k ) = (E[Y (i) 1 ], ..., E[Y (i) k ]). We suppose that (Y (i) 1 , ..., Y (i)
k ) is not a.s. a zero vector, which guarantees that the base process does not die out immediately. We let the vector process

N o (t) = (N o 1 (t), ..., N o k (t))
where N o j (t) represents the number of j type particles at time t. Also it is assumed that at time 0 the number of particles N o (0) is a random vector I = (I(1), . . . , I(k)) with distribution

I ∼ n∈N k p n δ n , (2)
for some probability vector (p n ) n∈N k , and where δ n is the Dirac distribution concentrated at n ∈ N k , and the second order moment holds

Assumption (A) : E(||I|| 2 ) = k j=1 E(I(j) 2 ) < +∞.
A central remark is that, given N o (0) = n 0 = (n 0 (1), ..., n 0 (k)) ∈ N k , N o (t) has the same distribution as the independent sum

k j=1 n 0 (j) l=1 N j,l (t), (3)
where, for each j = 1, ..., k, {N j,l (t), t ≥ 0} l∈N * are independent identically distributed (iid) copies of a generic multitype branching process {N j (t) = (N j 1 (t), ..., N j k (t)), t ≥ 0} with the branching mechanism given by (1) and N j i (t) represents the number of type i particles at time t produced from a type j particle generated at time 0, i.e. N j (0) is a vector of which jth entry is 1 and 0 elsewhere.

The CF of N o (t) is denoted by ϕ o t (s) := E[e <is,N o (t)> ] for s ∈ R k . According to [5, Chapter V],
{N o (t), t ≥ 0} is a continuous time multitype branching process (without immigration). Note that, in view of the representation (3), one has the following expression

ϕ o t (s) = E k j=1 ϕ j t (s) I(j) = n 0 ∈N k k j=1 ϕ j t (s) n 0 (j) p n 0 , (4)
where ϕ j t (s) := E(exp(< is, N j (t) >)).

We recall some useful results which will be used often in the subsequent study. First, it is convenient to introduce a k × k matrix A = (a ij ) i,j=1,...,k where the a ij 's are dened by

a ij = µ j (m ij -1 [i=j]
), i, j = 1, ..., k.

(

) 5
We suppose that A is regular i.e. all entries of the matrix exp(t 0 A) are positive for some t 0 > 0 (see [START_REF] Athreya | Branching processes[END_REF]Denition (10) p.202]). This entails that the largest eigenvalue ρ of A is with multiplicity 1. It is commonly known that we are in the subcritical case if ρ < 0, in the critical case if ρ = 0, and in the supercritical case if ρ > 0. We let u and v be the k × 1 right and left eigenvectors respectively, i.e. such that Au = ρu and v A = ρv , with positive entries, and normalized in such a way that < u, 1 >= 1 and < u, v >= 1. Then, in [5, Theorem 1 p.209] it was shown that {< u, N o (t)e -ρt >, t ≥ 0} is a martingale when N o (0) is a constant n 0 ∈ N k . In the case when N o (0) = I given by (2), it is clear that this property still holds because Assumption (A) in particular implies that I is integrable. Finally, we will suppose that the process is non singular, i.e. that for all i = 1, ..., k, p i (n) in ( 1) is not in the form

p i (n) = 1 [n i =1, n j =0, j =i] for all i, j = 1, ...,
ϕ W (x) = E[e ixW ] = n 0 ∈N k E[e ixW |N o (0) = n 0 ] p n 0 , x ∈ R.
Because of (3), one has

E[e ixW |N o (0) = n 0 ] = E[e ix k j=1 n 0 (j) l=1 W j,l ] = k j=1 ϕ j (x) n 0 (j) ,
where W j,l is the a.s. limit of N j,l (t)e -ρt as t → +∞. For each j = 1, ..., k, (W j,l ) l∈N is an iid sequence having a CF denoted by ϕ j (x) = E(e ixW j,l ), and the W j,l , j = 1, ..., k, l = 1, ..., n 0 (j) are independent. The vector ϕ(x) := (ϕ 1 (x), ..., ϕ k (x)) is in general not explicit but satises a particular integral equation (see [5, Eq.(28) p.206] for details). Finally, the CF of W is then given as a function of this vector ϕ(x)

ϕ W (x) = n 0 ∈N k k j=1 ϕ j (x) n 0 (j) p n 0 = E k j=1 ϕ j (x) I(j) , x ∈ R. (6)
We then move on to a multitype branching process with immigration which is the central stochastic process studied in this paper. Let us consider that new particles (immigrants)

arrive at time T i , i ≥ 1 according to a random vector I i = (I i (1), ..., I i (k)) having the same distribution (2) as I. Then it evolves according to the branching mechanism described at the beginning of this section. We assume that the immigration sequence (I i ) i∈N is iid. The vector process N (t) = (N 1 (t), ..., N k (t)) represents the number of each type of particles at time t dened as:

N (t) = S(t) i=1 N o,i (t -T i ), t ≥ 0, (7)
where {N o,i (t), t ≥ 0} i∈N are iid copies of {N o (t), t ≥ 0} with N o,i (0) = I i . The process {S(t), t ≥ 0} with S(0) = 0 is the arrival process for new particles associated with a nondecreasing sequence (T i ) i∈N with T 0 = 0 which represents the arrival time of the ith particle, and interarrival times (T i -T i-1 ) i∈N * . In other words, N o,i (t -T i ) is a vector of the number of particles in each system at time t generated from the I i (j) particles of type j, j = 1, ..., k, arrived at T i . Also, an underlying assumption is that N o,i j (t) = 0 when t < 0 for j = 1, ..., k, {S(t), t ≥ 0} and {N o,i (t), t ≥ 0}, i ∈ N, are processes independent of each other. Hence, N (t) is a continuous time multitype branching process with immigration given by the process {S(t), t ≥ 0}.

Immigration modelled by a nonhomogeneous Poisson process (NHPP)

We assume in this section that {S(t), t ≥ 0} is a NHPP with intensity λ(t) > 0 for t ≥ 0 being locally integrable, and set Λ(t) := t 0 λ(y)dy for t ≥ 0.

To study the asymptotic behavior of N (t) in (7) when t → +∞, we rst need the CF of N (t). The following result is an easy extension of [START_REF] Durham | A problem concerning generalized age-dependent branching processes with immigration[END_REF]Equation (2)]; see also [START_REF] Mitov | Multitype branching processes with inhomogeneous Poisson immigration[END_REF]Theorem 1] for a similar result that concerns the probability generating function of N (t).

Lemma 2. The CF of N (t) in ( 7) admits the following expression

ϕ t (s) = E[e <is,N (t)> ] = exp t 0 [ϕ o t-x (s)-1]λ(x)dx = exp t 0 [ϕ o x (s)-1]λ(t-x)dx , (8)
for all s ∈ R k .

Proof. Given S(t) = n, (T 1 , ..., T n ) are distributed as the ordered statistics (U (1) , ..., U (n) ) with (U 1 , ..., U n ) being independent with density λ(y)

Λ(t) 1 [0,t] (y), one nds ϕ t (s) = ∞ n=0 E exp < is, n l=1 N o,l (t -U (l) ) > e -Λ(t) (Λ(t)) n n! . Since n l=1 N o,l (t -U (l) ) = n l=1 N o,l (t -U l )
and by independence of (U 1 , ..., U n ) and the process {N o,l (t), t ≥ 0}, one obtains

ϕ t (s) = ∞ n=0 1 Λ(t) t 0 E[exp(< is, N o (t -y) >)]λ(y)dy n e -Λ(t) (Λ(t)) n n! .
A change of variable x := ty concludes the proof.

The following results show that the normalized process converges towards dierent limits depending on the assumptions regarding the intensity of the arrival process.

Theorem 3. Suppose that the intensity λ(t) of the NHPP {S(t), t ≥ 0} satises that e -ρt λ(t)

is integrable. Then e -ρt N (t) D -→ ∞ 0 e -ρz dY W z , t → +∞, (9)
where {Y W t , t ≥ 0} is a nonhomogeneous compound Poisson process with intensity λ(y) and jumps distributed as W v.

Proof. The proof is given in Section 3.1 Although Theorem 3 is valid when the eigenvalue ρ has any sign, it is especially interesting in the supercritical case ρ > 0, as e -ρt λ(t) ∈ L 1 (R + ) implies that the intensity λ(t) can, for example, grow exponentially as e δt for 0 ≤ δ < ρ. However Theorem 3 becomes less interesting in the critical case ρ = 0 or subcritical case ρ < 0, as the condition e -ρt λ(t) ∈ L 1 (R + ) roughly means that the intensity tends to 0 potentially very fast. Hence, the following supplementary results show how the renormalized process N (t) converges in distribution or in probability

• when the intensity grows exponentially in the critical or subcritical case,

• when the intensity grows exponentially as e δt with δ ≥ ρ in the supercritical case, complementing Theorem 3.

Theorem 4. Let us suppose that the intensity λ(t) of the NHPP {S(t), t ≥ 0} satises λ(t) ∼ λ ∞ e δt as t → +∞ for some δ ≥ 0 and λ ∞ > 0. Then the following convergences hold as t → +∞:

N (t) D -→ ν if ρ < 0 and δ = 0, (10) e -δt N (t) P -→ λ ∞ (δId -A) -1 E(I) if [ρ ≤ 0 and δ > 0] or [ρ > 0 and δ > ρ], (11) e -δt N (t) t P -→ λ ∞ uv E(I) if ρ = δ > 0, (12)
where Id is the identity matrix and ν is a distribution on R k + with CF given by

R k + e <is,x> ν(dx) = exp λ ∞ +∞ 0 [ϕ o y (s) -1] dy , s ∈ R k .
Proof. The proof is presented in Section 3.2.

In the following, we are especially interested in the particular critical case ρ = 0. From Theorems 3 and 4: we have that, if λ(t) is integrable then N (t) converges in distribution to lim t→+∞ Y W t in [START_REF] Cha | Characterization of the generalized Pólya process and its applications[END_REF], and if λ(t) ∼ λ ∞ e δt with δ > 0 then one has the convergence in probability of e -δt N (t) in [START_REF] Cha | Justifying the Gompertz curve of mortality via the generalized Polya process of shocks[END_REF]. We note that an intermediary case is worth exploring when λ(t) does not have an explosive behaviour and roughly speaking, does not converge to 0. This is the case if the associated Cesaro limit lim t→∞ Λ(t)/t = λ ∞ exists, hence some additional convergence result may be obtained. Before detailing this convergence result, we introduce the following quantities:

Q := 1 2 k i, ,n=1 ∂ 2 h i ∂z ∂z n (1, ..., 1)u u n v i , β j := k =1 µ -1 u v u j Q , j = 1, ..., k, (13)
β := k j=1 E(I(j))β j , (14)
c := ( k =1 µ -1 u v ) 2 Q , (15)
where we recall that h i (z) = h i (z 1 , ..., z k ) is the generating function associated with (p i (n)) n∈N k given in [START_REF] Allison | Estimation and testing for a Markov model of reinforcement[END_REF]. We note that the assumption that (Y

(i) 1 , ..., Y (i)
k ) is not zero for some i ∈ {1, ..., k} entails that Q is positive. The following result holds. Theorem 5. Let us assume that the moments of all orders of the random vector (Y

(i) 1 , ..., Y (i) k )
exist for all i = 1, ..., k and the intensity admits a Cesaro nite limit λ ∞ = lim t→∞ Λ(t)/t > 0. When ρ = 0 (critical case), one has the convergence in distribution as

N (t) t D -→ Zv ⊗ µ -1 , t → +∞, ( 16
)
where

Z is a rv distributed as Γ(λ ∞ β, c) with v ⊗ µ -1 = (v 1 µ -1 1 , ..., v k µ -1 k )
, β and c given by ( 14) and ( 15) respectively. Here, Γ(α, θ) denotes the gamma distribution with a shape parameter α and a rate parameter θ.

Proof. The proof is presented in Section 3.3. Thus, it turns out that, in the critical case the support of the limits ( 9), ( 11) and ( 16) are respectively included in the positive half line spanned by v, (δId -A) -1 E(I) and v ⊗ µ -1 . It is worthwhile to point out that all three results concern the critical case, but with dierent assumptions on the intensity function. Remark 6. When the intensity λ(t) is a constant value λ, Theorem 5 is the particular case of [START_REF] Weiner | A multi-type critical age-dependent branching process with immigration[END_REF]Theorem 2] which considers general interarrival times, with a slight change of notation (in that reference µ i stands for the mean lifetime of a type i particle, as opposed to µ -1 i in this paper). See also [START_REF] Vatutin | A critical Bellman-Harris branching process with immigration and several types of particles[END_REF]Theorem 1] for a similar result. When λ(t) converges to some limit λ ∞ , it converges towards the same limit in the sense of Cesaro and the limit in distribution [START_REF] Hyrien | Subcritical Sevastyanov branching processes with nonhomogeneous Poisson immigration[END_REF] corresponds to [START_REF] Mitov | Multitype branching processes with inhomogeneous Poisson immigration[END_REF]Theorem 8]. However, in this paper we introduce a dierent approach to the proof of Theorem 5 (given in Section 3.3). In contrary to [START_REF] Weiner | A multi-type critical age-dependent branching process with immigration[END_REF] which proved the result by showing that the joint moments of N (t)/t converge, our approach does not require renewal arguments and relevant results. In particular, we start directly with the CF (8) which is expressed handily in Lemma 2 and study its convergence. A similar approach was adopted in [START_REF] Mitov | Multitype branching processes with inhomogeneous Poisson immigration[END_REF] although the authors in [START_REF] Mitov | Multitype branching processes with inhomogeneous Poisson immigration[END_REF] started the proof from a uniform estimate from [START_REF] Sevast'yanov | Branching Processes[END_REF] for the probability generating function of {N o (t), t ≥ 0}.

Remark 7. From (52) in the proof of Theorem 5, the limiting distribution of ( 16) admits a similar integral form as the right-hand side of (9) which is available by applying Campbell's formula (the details are given in the beginning of Section 3.1). Indeed, one can check the equality in distribution of Zv ⊗ µ -1 and +∞ 0 e -t dY c t where {Y c t , t ≥ 0} is a compound Poisson process with intensity λ ∞ β and jumps distributed as χv ⊗ µ -1 with χ ∼ E(c).

We now proceed to the proofs of Theorems 3, 4 and 5.

Proof of Theorem 3

We start from the CF in [START_REF] Bühlmann | Mathematical methods in risk theory[END_REF], which entails that the CF of e -ρt N (t) is given by

ϕ t (se -ρt ) = exp t 0 [ϕ o x (se -ρt ) -1]λ(t -x)dx , s ∈ R k .
The main diculty in the proof is to show the following convergence:

t 0 [ϕ o x (se -ρt ) -1]λ(t -x)dx -→ +∞ 0 [ϕ W v (se -ρx ) -1]λ(x) dx, t → +∞, s ∈ R k , ( 17
)
where ϕ W v is the CF of W v.

[ϕ W v (se -ρy ) -1]λ(y)dy} is the CF of ∞ 0 e -ρz dY W z ,
where {Y W t , t ≥ 0} is a nonhomogeneous compound Poisson process with intensity λ(y) and jumps distributed as W v. Hence one has the convergence in distribution of e -ρt N (t) towards ∞ 0 e -ρz dY W z if (17) holds. This proves (9).

In order to prove [START_REF] Greenwood | An inquiry into the nature of frequency distributions representative of multiple happenings with particular reference to the occurrence of multiple attacks of disease or of repeated accidents[END_REF] we need to exploit the convergence N o (y)e -ρy -→ W v a.s. as y → +∞ given in Lemma 1. Studying ( 17) is equivalent to analyzing the limit as t → +∞ of

Q t := t 0 [ϕ o t-x (se -ρt ) -1]λ(x) dx = t 0 {E[exp(< is, N o (t -x)e -ρt >)] -1}λ(x) dx. ( 18
)
That is, Q t may be expressed as

Q t := Q 1,t + Q 2,t , (19)
where

Q 1,t := t 0 E exp < is, N o (t -x) e ρ(t-x) e -ρx > -exp(< is, W ve -ρx >) λ(x) dx, (20)
Q 2,t := t 0 {E[exp(< is, W ve -ρx >)] -1}λ(x) dx. (21)
We then separately examine the limits of ( 20) and ( 21). In the end, it will be shown that [START_REF] Kyprianou | Introductory lectures on uctuations of Lévy processes with applications[END_REF] tends to 0 and (21) tends to the right-hand side of [START_REF] Greenwood | An inquiry into the nature of frequency distributions representative of multiple happenings with particular reference to the occurrence of multiple attacks of disease or of repeated accidents[END_REF].

Limit of Q 1,t in (20) as t → +∞. We shall utilize the following basic inequality in the subsequent proof:

|e ia -e ib | ≤ |a -b|, a ∈ R, b ∈ R, (22)
and also we have

|e ia -e ib | ≤ 2. Hence |e ia -e ib | ≤ |a -b| ∧ 2 for a ∈ R and b ∈ R. We then deduce that |Q 1,t | ≤ t 0 E < s, N o (t -x) e ρ(t-x) e -ρx > -< s, W ve -ρx > ∧ 2 λ(x)dx = ∞ 0 1 [0≤x≤t] E < s, N o (t -x) e ρ(t-x) e -ρx > -< s, W ve -ρx > ∧ 2 λ(x)dx. ( 23
)
By the dominated convergence theorem, we will show in the following that [START_REF] Gat | Extending the Yule process to model recurrent pipe failures in water supply networks[END_REF] tends to zero as t → +∞. From the pointwise convergence as t → +∞ in Lemma 1 with the help of the dominated convergence theorem, one nds that the integrand goes to zero i.e.

1 [0≤x≤t] E e -ρx < s, N o (t -x) e ρ(t-x) > -e -ρx < s, W v > ∧ 2 λ(x) -→ 0, t → ∞, ∀x ≥ 0.
We now nd an upper bound function f (x) ≥ 0 of this integrand such that +∞ 0 f (x)dx < +∞. Recall that u is an eigenvector with positive entries u i for i = 1, ..., k such that Au = ρu (where the elements of the matrix A are dened in [START_REF] Athreya | Branching processes[END_REF]). Since u i > 0 for all i, there exists some constant κ > 0 which is large enough satisfying

0 ≤ |s j | ≤ κu j , ∀j = 1, ..., k, (24)
where we recall that the vector s = (s 1 , ..., s k ) is xed. For example, κ can be chosen as max j=1,...,k |s j |/u j . Also, note that E

[(X + Y ) ∧ 2] ≤ (E[X] + E[Y ]) ∧ 2 for
1 [0≤x≤t] e -ρx E < s, N o (t -x) e ρ(t-x) > + e -ρx E[| < s, W v > |] ∧ 2 λ(x) ≤ 1 [0≤x≤t] e -ρx E < |s|, N o (t -x) e ρ(t-x) > + e -ρx E[< |s|, W v >] ∧ 2 λ(x) ≤ 1 [0≤x≤t] e -ρx κ E < u, N o (t -x) e ρ(t-x) > + e -ρx E[< |s|, W v >] ∧ 2 λ(x), (25)
where the rst equality is due to the fact that

N o (t-x)
e ρ(t-x) and W v have nonnegative entries, and the last inequality is due to [START_REF] Mitov | Multitype branching processes with inhomogeneous Poisson immigration[END_REF]. The rst expectation in ( 25) is E[< u, N o (0)/e ρ×0 >] because of the martingale property and in turn, it is equal to < u, E(I) > because of N o (0) = I. The second expectation is some nite constant. Therefore we conclude that, for some constants K > 0 and K * > 0, [START_REF] Pólya | Sur quelques points de la théorie des probabilités[END_REF] is bounded as

1 [0≤x≤t] e -ρx κE < u, N o (t -x) e ρ(t-x) > + e -ρx E[< |s|, W v >] ∧ 2 λ(x) = 1 [0≤x≤t] [(κ < u, E(I) > e -ρx + Ke -ρx ) ∧ 2]λ(x) ≤ K * e -ρx λ(x) := f (x).
It is now shown that the integrand in [START_REF] Gat | Extending the Yule process to model recurrent pipe failures in water supply networks[END_REF] tends to 0 as t → +∞ for a xed x and is dominated by the function f (x) which is integrable by assumption. Therefore, by the dominated convergence theorem we conclude that [START_REF] Gat | Extending the Yule process to model recurrent pipe failures in water supply networks[END_REF]

goes to 0 as t → ∞, which implies that Q 1,t in (20) veries lim t→∞ Q 1,t = 0. Limit of Q 2,t in (21) as t → +∞. Let us show that |{ϕ W v (se -ρx ) -1}λ(x)| is upper
bounded by some integrable function. With the help of ( 22), the following inequality holds:

exp(< is, W ve -ρx >) -1 ≤ < s, W ve -ρx > ≤ e -ρx < |s|, W v > .
We then arrive at the following bound

{ϕ W v (se -ρx ) -1}λ(x) = {E[exp(< is, W ve -ρx >)] -1}λ(x) ≤ e -ρx λ(x)E [< |s|, W v >] ,
which indeed is integrable by the integrability assumption for e -ρx λ(x). Then by the dominated convergence theorem one nds

Q 2,t -→ +∞ 0 {ϕ W v (se -ρx ) -1}λ(x) dx, t → ∞.
Hence, combining the limits of Q 1,t and Q 2,t , ( 17) is proved. Consequently, this completes the proof.

Proof of Theorem 4

Proving the limit of N (t) in [START_REF] Cha | New shock models based on the generalized Polya process[END_REF] as t → +∞. First, we recall that ρ < 0 and δ = 0, and the intensity satises lim t→∞ λ(t) = λ ∞ . We proceed to show that

t 0 [ϕ o x (s) -1]λ(t -x)dx in Lemma 2 converges to λ ∞ ∞ 0 [ϕ o
x (s)-1]dx as t → +∞ by a dominated convergence argument. Let A λ ≥ 0 be large enough such that C λ := sup y≥A λ λ(y) is nite, and t ≥ A λ . We start by

writing t 0 [ϕ o x (s) -1]λ(t -x)dx = ∞ 0 1 [0≤x≤t] [ϕ o x (s) -1]λ(t -x)dx = ∞ 0 1 [0≤x≤t-A λ ] [ϕ o x (s) -1]λ(t -x)dx + ∞ 0 1 [t-A λ <x≤t] [ϕ o x (s) -1]λ(t -x)dx. ( 26
)
Let us study each term on the right-hand side of [START_REF] Qi | Estimating the propagation of interdependent cascading outages with multi-type branching processes[END_REF]. Using the inequality in [START_REF] Mitov | Multitype branching processes with inhomogeneous Poisson immigration[END_REF], one nds

|ϕ o x (s) -1| = |E[e <is,N o (x)> ] -1| ≤ E[|e <is,N o (x)> -1|] ≤ E[| < s, N o (x) > |] ≤ E[< |s|, N o (x) >] ≤ κE[< u, N o (x) >] = κe ρx E[< u, N o (x)e -ρx >] = κe ρx E[< u, N o (0) >] = κe ρx < u, E(I) >, (27)
where we recall that κ = max j=1,...,k |s j |/u j for example. Since λ(.) is locally integrable (so that A λ 0 λ(x)dx is nite), and e ρx ≤ e ρ(t-A λ ) for t -A λ ≤ x and ρ < 0, the second term on the right-hand side of ( 26) thus veries

∞ 0 1 [t-A λ <x≤t] [ϕ o x (s) -1]λ(t -x)dx ≤ ∞ 0 1 [t-A λ <x≤t] |ϕ o x (s) -1|λ(t -x)dx ≤ κ < u, E(I) > ∞ 0 1 [t-A λ <x≤t] e ρx λ(t -x)dx ≤ κ < u, E(I) > e ρ(t-A λ ) ∞ 0 1 [t-A λ <x≤t] λ(t -x)dx = κ < u, E(I) > e ρ(t-A λ ) A λ 0 λ(x)dx -→ 0, t → ∞, t ≥ A λ . ( 28
)
Concerning the rst term on the right-hand side of (26), from [START_REF] Resing | Polling systems and multitype branching processes[END_REF] and the denition of C λ we nd that

1 [0≤x≤t-A λ ] |ϕ o x (s) -1|λ(t -x) ≤ 1 [0≤x≤t-A λ ] κe ρx < u, E(I) > λ(t -x) ≤ C λ κe ρx < u, E(I) > . ( 29
)
Since ρ < 0, the right-hand side of ( 29) is integrable, and

lim t→+∞ 1 [0≤x≤t-A λ ] [ϕ o x (s) -1]λ(t - x) = [ϕ o
x (s) -1]λ ∞ , the dominated convergence theorem entails that the rst term on the right-hand side of [START_REF] Qi | Estimating the propagation of interdependent cascading outages with multi-type branching processes[END_REF]

converges to λ ∞ ∞ 0 [ϕ o
x (s) -1]dx as t → +∞. Hence, combining with the limit obtained in (28), we deduce that the limit of the left-hand side of [START_REF] Qi | Estimating the propagation of interdependent cascading outages with multi-type branching processes[END_REF]

as t → ∞ is λ ∞ ∞ 0 [ϕ o x (s) -1]dx. We now argue that the function s → exp λ ∞ ∞ 0 [ϕ o x (s) -1]dx is continuous in a neighborhood of 0 ∈ R k such as s ∈ (-1, 1) k , which by Lévy's continuity theorem implies that exp λ ∞ ∞ 0 [ϕ o x (s) -1]dx is the CF of some distribution ν on R k + .
To be more precise, we observe that, when s ∈ (-1, 1) k , the upper bound κ in (29) can be chosen independently from s ∈ (-1, 1) k as max j=1,...,k 1/u j , a dominated convergence argument subsequently implies this continuity property.

Proving the limit of e -δt N (t) in [START_REF] Cha | Justifying the Gompertz curve of mortality via the generalized Polya process of shocks[END_REF] as t → +∞. We assume here that ρ ≤ 0 and δ > 0

or that ρ > 0 and δ > ρ. First, the CF of e -δt N (t) is given as e Rt where R t can be written as below:

R t = t 0 [ϕ o t-x (se -δt ) -1]λ(x) dx = t 0 {E[exp(< ise -δt , N o (t -x) >)] -1}λ(x) dx (30) = R 1,t + R 2,t , R 1,t := ∞ 0 1 [0≤x≤t] E exp < ise -δt , N o (t -x) > -1-< ise -δt , N o (t -x) > λ(x) dx, (31)
R 2,t := ∞ 0 1 [0≤x≤t] E < ise -δt , N o (t -x) > λ(x) dx. ( 32
)
Similar to [START_REF] Jeanblanc | Poisson Processes and Ruin Theory[END_REF], we then study the limits of R

|e ix -1| ≤ |x|) |e ix -1 -ix| ≤ |e ix -1| + |ix| ≤ 2|x| ≤ 2|x| 2 if |x| > 1, so that we have the following general inequality |e ix -1 -ix| ≤ 2|x| 2 , ∀x ∈ R.
The above result, combined with the Cauchy-Schwarz inequality, yields the upper bound for

|R 1,t | given by |R 1,t | ≤ 2 ∞ 0 1 [0≤x≤t] E | < se -δt , N o (t -x) > | 2 λ(x) dx ≤ 2||s|| 2 ∞ 0 1 [0≤x≤t] e -2δt E ||N o (t -x)|| 2 λ(x) dx. ( 33
)
We here separate the cases ρ < 0, ρ = 0 and ρ > 0, with last case requiring the additional constraint δ > ρ. If ρ < 0, the growth rate in Lemma 16 (with proof in Appendix A) implies that E ||N o (tx)|| 2 ≤ Ce ρ(t-x) for some constant C > 0. Also, the assumption λ(x) ∼ λ ∞ e δx as x → ∞, in particular implies that λ(x) is bounded by e δx up to a constant, hence one gets from (33) that for some (dierent) constant C > 0,

|R 1,t | ≤ C ∞ 0 1 [0≤x≤t] e -2δt e ρ(t-x) e δx dx = Ce (-2δ+ρ)t t 0 e (-ρ+δ)x dx = C -ρ + δ [e -δt -e (-2δ+ρ
)t ] -→ 0 as t → +∞.

If ρ = 0, the growth rate in Lemma 16 implies that E ||N o (tx)|| 2 is an O((tx) + 1), hence for some constant C > 0 we have

|R 1,t | ≤ C ∞ 0 1 [0≤x≤t] e -2δt [(t -x) + 1]e δx dx ≤ C(t + 1) t 0 e -2δt e δx dx = C(t + 1) δ [e -δt -e -2δt ] -→ 0 as t → +∞.
Finally, if ρ > 0 then the growth rate in Lemma 16 implies that E ||N o (tx)|| 2 is less than e 2ρ(t-x) up to a constant, hence for some constant C > 0 we have

|R 1,t | ≤ C ∞ 0 1 [0≤x≤t] e -2δt e 2ρ(t-x) e δx dx = Ce 2(-δ+ρ)t t 0 e (-2ρ+δ)x dx, (34)
which can be shown approaching 0 as t → +∞, using δ > ρ.

Limit of R 2,t in [START_REF] Wasserman | Distinguishing between stochastic models of heterogeneity and contagion[END_REF] as t → +∞, and conclusion. From [5, p.208], we know that the mean matrix of the multitype process N o (z) is expressed as E[N o (z)] = e Az E(I) where the matrix A is dened in [START_REF] Athreya | Branching processes[END_REF]. Therefore, R 2,t can be expressed, after some manipulation, as

R 2,t = t 0 e -δt < is, e A(t-x) E(I) > λ(x)dx = ∞ 0 1 [0≤x≤t] e -δt < is, e Ax E(I) > λ(t -x)dx. ( 35
)
We now wish to use the dominated convergence theorem in order to nd the limit in [START_REF] Willmot | Distributional analysis of a generalization of the Polya process[END_REF].

Upper bounding λ(x) by Ce δx for some constant C > 0 results in

1 [0≤x≤t] e -δt < is, e Ax E(I) > λ(t -x) ≤ C < |s|, e (A-δId)x E(I) >,
which is integrable because δ > ρ (in either case ρ ≤ 0 or ρ > 0) so that all eigenvalues of A -δId have negative real parts. Also, the assumption that λ(x) ∼ λ ∞ e δx as x → ∞, results in 1 [0≤x≤t] e -δt < is, e Ax E(I) > λ(tx) -→ λ ∞ < is, e (A-δId)x E(I) > as t → +∞, for all

x ≥ 0. Hence we deduce from (35) that

R 2,t -→ λ ∞ ∞ 0 < is, e (A-δId)x E(I) > dx =< is, λ ∞ ∞ 0 e (A-δId)x dx E(I) > =< is, λ ∞ (δId -A) -1 E(I) > .
Since R 1,t -→ 0 as t → +∞, one arrives at the convergence of the CF of e -δt N (t) to exp(< is, λ ∞ (δId -A) -1 E(I) >), so that e -δt N (t) converges in distribution (or, equivalently, in probability) towards λ ∞ (δId -A) -1 E(I). Hence (ii) in ( 11) is proved.

Proving the limit of e -δt N (t) t in [START_REF] Durham | A problem concerning generalized age-dependent branching processes with immigration[END_REF] as t → +∞. We assume the supercritical case ρ > 0 and ρ = δ. In this case, instead of (30), we consider the quantity R t :

= t 0 [ϕ o t-x (se -δt /t) - 1]λ(x) dx = t 0 {E[exp(< ise -δt /t, N o (t -x) >)] -1}λ(x) dx such that the CF of e -δt N (t)
t is equal to e Rt , which is similarly decomposed as in [START_REF] Vatutin | A critical Bellman-Harris branching process with immigration and several types of particles[END_REF] and [START_REF] Wasserman | Distinguishing between stochastic models of heterogeneity and contagion[END_REF]

as R t = R 1,t + R 2,t , where R 1,t := ∞ 0 1 [0≤x≤t] E exp < ise -δt /t, N o (t -x) > -1-< ise -δt /t, N o (t -x) > λ(x) dx, (36) R 2,t := ∞ 0 1 [0≤x≤t] E < ise -δt /t, N o (t -x) > λ(x) dx. (37)
Utilizing the inequality in [START_REF] Weiner | A multi-type critical age-dependent branching process with immigration[END_REF], one obtains, for some constant C > 0,

|R 1,t | ≤ 1 t 2 C ∞ 0 1 [0≤x≤t] e -2δt e 2ρ(t-x) e δx dx = C 1 t 2 t 0
e -δx dx -→ 0, thus (36) tends to 0 as t → +∞. Using the expression E[N o (z)] = e Az E(I) as in [START_REF] Willmot | Distributional analysis of a generalization of the Polya process[END_REF] and the change of variable z := 1 -x/t together with the assumption δ = ρ, we can rewrite (37) as

R 2,t = 1 t t 0 e -δt < is, e A(t-x) E(I) > λ(x)dx = < is, 1 0 e At(1-z) e -ρt λ(tz)dz E(I) > = < is, 1 0 e (A-ρId)t(1-z) e -ρtz λ(tz)dz E(I) > . (38)
We now wish to investigate (38) when t → +∞. Since A is regular, the Perron Frobenius theory entails that e (A-ρId)x converges to uv as x → +∞, see e.g. [START_REF] Athreya | Branching processes[END_REF]Limit (17) p.203]. Hence, for all z ∈ (0, 1) one has lim t→∞ e (A-ρId)t(1-z) = uv . Also, the assumption λ(x) ∼ λ ∞ e δx as x → ∞ with δ = ρ implies that lim t→∞ e -ρtz λ(tz) = λ ∞ for all z ∈ (0, 1), so that by the dominated convergence theorem we may let t → +∞ in (38) and obtain

R 2,t -→< is, λ ∞ uv E(I) >, t → +∞.
Therefore, since R t = R 1,t + R 2,t with lim t→∞ R 1,t = 0, one concludes the convergence (12).

Proof of Theorem 5 in the critical case ρ = 0

We start with Lemma 2, from which we deduce that the CF of N (t)/t admits the expression

E exp < is, N (t) t > = E[exp(< t -1 is, N (t) >)] = exp t 0 [ϕ o t-y (t -1 s) -1]λ(y)dy . ( 39
)
We thus study

t 0 [ϕ t-y (t -1 s) -1]λ(y)dy = t 0 E exp < is, N o (t -y) t > -1 λ(y)dy = Λ(t) 0 E exp < is, N o (t -Λ -1 (y)) t > -1 dy = 1 0 Λ(t) E exp < is, N o (t -Λ -1 (Λ(t)x)) t > -1 dx := - 1 0 γ t (x)dx, (40)
where Λ -1 (.) is the inverse of the function Λ(.) (invertible as it is assumed that λ(t) > 0 for all t ≥ 0), the second last equality is due to a change of variable with x := y/Λ(t) and γ t (x)

is given by

γ t (x) := Λ(t) E 1 -exp < is, N o (t -Λ -1 (Λ(t)x)) t > . (41)
We note that the assumption lim t→+∞ Λ(t

)/t = λ ∞ implies that lim t→+∞ Λ -1 (t)/t = λ -1 ∞ ,
which is in turn equivalent to

Λ -1 (t) ∼ λ -1 ∞ t, i.e. Λ -1 (t) = λ -1 ∞ t + η(t)t, (42)
where lim t→∞ η(t) = 0.

In the following, we shall prove by the dominated convergence theorem that the right-hand side of (40) has the following limit

- 1 0 γ t (x)dx -→ λ ∞ β ∞ 0 E[exp(< is, e -y X >) -1]dy, t → +∞, ( 43
)
where β is given by [START_REF] Hyrien | Stochastic modeling of stress erythropoiesis using a two-type age-dependent branching process with immigration[END_REF]. In here we have

X = χv ⊗ µ -1 ∈ [0, +∞) k , (44)
where χ ∼ E(c) for c > 0 given by ( 15) and the survival function of X given by

P (X > z) = exp -c max i=1,...,k z i v i µ -1 i , z = (z 1 , ..., z k ) ∈ [0, +∞) k . ( 45
)
The proof is decomposed in the following steps.

Step 1: Dominating the integrand in (40). Using the basic inequality |e ix -1| ≤ x, we have for all t ≥ 0 and x ∈ (0, 1),

s = (s 1 , ..., s k ) ∈ R k that 1 -exp < is, N o (t -Λ -1 (Λ(t)x)) t > ≤ < |s|, N o (t -Λ -1 (Λ(t)x)) t > .
Using [START_REF] Mitov | Multitype branching processes with inhomogeneous Poisson immigration[END_REF], one nds < |s|,

N o (t-Λ -1 (Λ(t)x)) t >≤ κ < u, N o (t-Λ -1 (Λ(t)x))
t > . Hence, taking the expectation and multiplying by Λ(t) on both sides result in

|γ t (x)| ≤ Λ(t)κ E < u, N o (t -Λ -1 (Λ(t)x)) t > ≤ C λ κ E[< u, N o (t -Λ -1 (Λ(t)x)) >] =C λ κ E[< u, N o (0) >] = C λ κ < u, E(I) >, t ≥ t 0 , (46)
where the rst equality is obtained by the martingale argument and C λ := sup t≥t 0 Λ(t)/t which is nite for a large enough t 0 . Since C λ and κ are constants (independent of t and x), the integrand in (40) is dominated by some constant independent from t ≥ t 0 and x ∈ (0, 1).

Step 2: Almost sure limit of the integrand in (40). Let us now prove the following convergence for (41):

γ t (x) -→ γ(x) := λ ∞ β 1 -x E[1 -exp(< is, (1 -x) X >)], (47)
for a xed x ∈ (0, 1) and s = (s 1 , ..., s k ) ∈ R k as t → +∞, where β is dened by [START_REF] Hyrien | Stochastic modeling of stress erythropoiesis using a two-type age-dependent branching process with immigration[END_REF]. To show this, we rst express γ t (x) in (41) using the representation (3) as

γ t (x) = n 0 =(n 0 (1),...,n 0 (k))∈N k γ t (x, n 0 )p n 0 , (48)
where

γ t (x, n 0 ) := Λ(t) 1 -E exp < is, N o (t -Λ -1 (Λ(t)x)) t > N o (0) = n 0 = Λ(t) 1 - k j=1 ψ j (t, x) n 0 (j) , (49)
ψ j (t, x) := E exp < is, N j (t -Λ -1 (Λ(t)x)) t > , j = 1, ..., k. (50)
In the following, we dene for all J ⊂ {1, ..., k} and z = (z 1 , ..., z k ) ∈ R * k + the vector z J of which the jth entry z J j is z j if j ∈ J, and some arbitrary negative value (e.g. -1) otherwise. The following lemma, proved in Appendix B, gives the asymptotic behaviour of ψ j (t, x), j = 1, ..., k, as t → ∞, which helps us prove (47). Lemma 8. The following limit holds for all j = 1, ..., k and x ∈ (0, 1):

t[ψ j (t, x) -1] -→ 1 1 -x β j [1 -E (exp[< is, (1 -x)X >])] , t → ∞. ( 51
)
Note that (51) in particular implies that ψ j (t, x) -

1 = O(1/t) -→ 0 as t → ∞, which, plugged into (49) together with Λ(t) ∼ λ ∞ t as t → +∞, yields γ t (x, n 0 ) ∼ λ ∞ k j=1 n 0 (j) t[ψ j (t, x) -1] , t → +∞.
This implies, along with (51), that (48) converges as

γ t (x) -→ λ ∞ 1 -x k j=1 β j n 0 ∈N k n 0 (j)p n 0 [1 -E (exp[< is, (1 -x)X >])] , t → +∞,
and thus (47) holds.

Step 3: Proof of (43). Thanks to (46) and (47), by the dominated convergence theorem, one thus deduces that (40) converges as t → +∞ to -

1 0 γ(x)dx = - 1 0 λ ∞ β 1 -x E[1 -exp(< is, (1 -x) X >)]dx,
which results in (43) after a change of variable y := -ln(1 -x).

Step 4: End of proof. From (39) with the convergence results of (40) towards (43), one nds that

E exp < is, N (t) t > -→ exp λ ∞ β ∞ 0 E[exp(< is, e -y X >) -1]dy , t → +∞, ( 52
) for s ∈ R k . Since X = χv ⊗ µ -1 with χ ∼ E(c), one computes that E[exp(< is, e -t X > ) -1] = E[exp(< is, v ⊗ µ -1 > χe -t ) -1] = e -t <is,v⊗µ -1 >
c-e -t <is,v⊗µ -1 >

. In turn, changing of variable z := e -t < is, v ⊗ µ -1 > yields that the right-hand side of the above convergence is the CF equivalent to c c-<is,v⊗µ -1 > λ∞ β , which indeed is the CF of Zv ⊗ µ -1 in ( 16). This completes the proof.

Immigration modelled by a contagious Poisson process (CPP)

As discussed in Section 1, CPP is a special case of GPP which has become a well-known contagion model when the transition intensity in the non-homogeneous birth process is a linear function of the current state multiplied by a function of the current time. In this section, we now assume that the arrival process {S(t), t ≥ 0} is a CPP and note that S(0 -) = 0. We rst consider GPP and obtain the CF of N (t). This will in turn imply that the one for a CPP is available. It is known that GPP (or a positive contagion model in [START_REF] Bühlmann | Mathematical methods in risk theory[END_REF][START_REF] Willmot | Distributional analysis of a generalization of the Polya process[END_REF]), is a particular case of self exciting counting process with intensity rate λ(t) satisfying

λ(t) = [aS(t -) + b]λ t , a > 0, b > 0, (53)
for some underlying function λ t > 0 for t ∈ (0, +∞) which is assumed to be continuous and integrable over nite ranges. Let us denote Λ t = t 0 λ y dy for t ≥ 0. When b = 1, this arrival process was referred to as the Linear Extension of the Yule Process by [START_REF] Gat | Extending the Yule process to model recurrent pipe failures in water supply networks[END_REF]. Hence, the intensity increases linearly with the number of arrivals at time t. This explains why these models could be appropriate for the situations where the arriving particles representing cells infected by a rapidly expanding disease contaminate other cells in an organism modelled by a certain network mechanism, or where the occurrence of shocks causes outages of interconnected lines in a power system as studied in [START_REF] Qi | Estimating the propagation of interdependent cascading outages with multi-type branching processes[END_REF]. In particular, we shall focus on the case when λ t = λ > 0 is constant. In this case, {S(t), t ≥ 0} is called a CPP in [START_REF] Allison | Estimation and testing for a Markov model of reinforcement[END_REF] and [START_REF] Wasserman | Distinguishing between stochastic models of heterogeneity and contagion[END_REF].

Let us start by establishing the CF of N (t) as obtained in Lemma 2 for the NHPP immigration.

Lemma 9. When the new particle arrives according to a GPP with the intensity rate given in (53), the CF of N (t) in ( 7) admits the following expression

ϕ t (s) = 1 - t 0 [ϕ o t-y (s) -1]aλ y e aΛy dy -b/a . ( 54
)
In particular, when λ t = λ in (53) (i.e. CPP case), ( 54) is simplied as

ϕ t (s) = 1 - t 0 [ϕ o t-y (s) -1]aλe aλy dy -b/a . ( 55
)
Proof. It is known that the marginal distribution of S(t) is expressed as a negative binomial distribution with Λ t = t 0 λ y dy (e.g. [9, Theorem 1(i)]) given by

p t (n) := P(S(t) = n) = Γ(b/a + n) Γ(b/a)n! (1 -e -aΛt ) n (e -aΛt ) b/a ,
where Γ(z) = ∞ 0 x z-1 e -x dx, for all complex number z with positive real part, is the gamma function. The above is a negative binomial distribution (r, p) where r = b/a and p = 1 -e -aΛt . Its probability generating function is

P t (z) = ∞ n=0 z n p t (n) = ( 1-p 1-pz ) r for all complex number z verifying |z| < p -1 .
Then, from [21, Section 3.2], the CF of N (t) can be expressed as a compound negative binomial distribution

ϕ t (s) = P t ( f t (s)), (56)
where the CF of the secondary distribution is given by

f t (s) = t 0 q t (y)ϕ o t-y (s)dy, (57)
with q t (y) = aλ y e aΛy e aΛt -1 , 0 ≤ y ≤ t.

(58) 56) is obtained as That is,

Since P t (z) = ( 1-p 1-pz ) r , (
ϕ t (s) = e -aΛt 1 -(1 -e -aΛt ) f t (s) b/a = e
ϕ t (s) = 1 + t 0 [1 -ϕ o t-y (s)]aλ y e aΛy dy -b/a , (59)
or equivalently (54).

In the case of a constant baseline intensity λ t = λ, taking the expectation on both sides of (53 We note that there is some resemblance between this exponential expression in (60) in the GPP case and the exponential asymptotic form λ(t) ∼ λ ∞ e δt of the (deterministic) intensity appeared in Theorems 3 and 4 in the NHPP case. However, due to the randomness in time feature of the intensity in this case, we can observe dierent limiting behaviors of the branching process N (t) with the CPP immigration. More precisely, in the following it is shown that the distributional behaviour changes depending on whether the largest eigenvalue ρ of A is less than, larger than, or equal to aλ. The main result of this section is given in the following theorem.

) yields E[λ(t)] = aλE[S(t -)]+λb. Since E[S(t -)] = E[S(t)] and {S(t)- t 0 λ(s)ds, t ≥ 0} is a martingale (see e.
Theorem 10. Suppose that λ t = λ , a is dened by (53) and ρ is the largest eigenvalue of the matrix A with elements a ij dened by [START_REF] Athreya | Branching processes[END_REF].

(1) If ρ > aλ, then e -ρt N (t)

D -→ Z T v, t → +∞, (61)
where v is the left eigenvector of the matrix A, T ∼ Γ(b/a, 1) and {Z t , t ≥ 0} is an independent Lévy process with characteristic exponent ψ(x)

:= R (1 -exp [ixz]) Π(dz), x ≥ 0. Here, Π(.)
is dened by

Π(dz) := E W aλ/ρ 1 [W ≥z] aλ ρ z -aλ/ρ-1 1 [0<z<+∞] dz, (62)
and we recall that W is characterized by its CF given by ( 6).

(2) If ρ < aλ, then e -aλt N (t)

D -→ Z γ, t → +∞, (63)
where Z is a rv distributed as Γ(b/a, 1) and γ is the vector dened by

γ := aλ(aλId -A) -1 E(I) (64)
(3

) If ρ = aλ, then N (t) t e -aλt D -→ Zv, t → +∞, (65)
where Z is a rv distributed as Γ(b/a, E[W ]aλ).

Remark 11. In the case ρ > aλ we may note that, since Π(.) dened by (62) has support on (0, +∞) and veries (0,+∞) min(1, z)Π(dz) < +∞ (precisely because of the condition ρ > aλ), the underlying Lévy process {Z t , t ≥ 0} in (61) belongs to the class of subordinators according to [20, Lemma 2.14, p.55].

Remark 12. As shown in (60), the expected intensity in the present CPP case has an exponential form, hence it is natural to compare the limiting convergence results in Theorem 10 to those in Theorems 3 and 4 in the NHPP case with asymptotic intensity λ(t) ∼ λ ∞ e δt for δ > 0. With an analogy between aλ and δ, Table 1 summarizes the dierent directions of the supports of the limiting distributions obtained in ( 11), ( 12) and ( 9) for ρ < δ, ρ = δ and ρ > δ in the NHPP case, to those in (63), ( 65) and (61) for ρ < aλ, ρ = aλ and ρ > aλ in the CPP case. Each value in Table 1, that all belong to R k + , roughly shows the position in which the renormalized process N (t) is located asymptotically in the corresponding case. Interestingly,

ρ < ξ ρ = ξ ρ > ξ NHPP: ξ = δ (ξId -A) -1 E(I) u v CPP: ξ = aλ (ξId -A) -1 E(I) v v
Table 1 Direction of limiting distribution the directions are the same except for ρ = δ in the NHPP case and ρ = aλ in the CPP case, which are respectively given by the vectors u and v.

The proofs of each case in Theorem 10 are provided in the following Section 4.1, Section 4.2, and Section 4.3 respectively.

Proof of Theorem 10 in the case ρ > aλ

In (55), with the normalizing function g(t) = e ρt we get

ϕ t (s/g(t)) = ϕ t (se -ρt ) = 1 + t 0 [1 -ϕ o t-y (se -ρt )]aλe aλy dy -b/a (66)
for all s ∈ R k . The proof is divided into two steps as follows.

Step 1: Studying the convergence of ϕ t (se -ρt ) as t → +∞. It is convenient to introduce the function

Ξ t,s := t 0 [1 -ϕ o t-y (se -ρt )]aλe aλy dy = ∞ 0 1 [0<y<t] E 1 -exp < is, N o (t -y
)/e ρ(t-y) > e -ρy aλe aλy dy, (67) so that ϕ t (se -ρt ) = {1 + Ξ t,s } -b/a , t ≥ 0. Thus, studying the limit of ϕ t (se -ρt ) as t → +∞ essentially requires nding lim t→+∞ Ξ t,s , which will be completed by the dominated convergence theorem. First note that for all y ∈ (0, +∞) one has that N o (ty)/e ρ(t-y) -→ W v, t → ∞, a.s. from Lemma 1. One has by the dominated convergence theorem for a xed y ∈ (0, +∞) and that

E 1 -exp < is, N o (t -y)/e ρ(t-y) > e -ρy -→ E 1 -exp < is, v > W e -ρy , t → +∞. (68)
Using [START_REF] Laskin | Fractional Poisson process[END_REF], the integrand in (67) is upper bounded in modulus as -y) > e (aλ-ρ)y .

1 [0<y<t] E 1 -exp < is, N o (t -y)/e ρ(t-y) > e -ρy aλe aλy ≤ 1 [0<y<t] E < |s|, N o (t -y)/e ρ(t-y) > e -ρy aλe aλy = aλ1 [0<y<t] E < |s|, N o (t -y)/e ρ(t
By the similar martingale argument applied to the one leading to (29) for example, and since Assumption (A) holds, one can show that 1 [0<y<t] E < |s|, N o (ty)/e ρ(t-y) > is upper bounded by some constant K which is independent of t and y. That is,

0 ≤ 1 [0<y<t] E 1 -exp < is, N o (t -y
)/e ρ(t-y) > e -ρy aλe aλy ≤ aλKe (aλ-ρ)y , (69) which is integrable over y ∈ (0, +∞) when ρ > aλ. Hence, using (68), (69) and the dominated convergence theorem, we arrive at

Ξ t,s -→ Ξ ∞,s := ∞ 0 E 1 -exp < is, v > W e -ρy aλe aλy dy, t → +∞, (70)
so that the renormalized CF in (66) converges as

ϕ t (se -ρt ) -→ φ(s) := {1 + Ξ ∞,s } -b/a , t → +∞. (71)
Step 2: Identifying the CF φ(s). In order to interpret (71) as the convergence towards some known distribution, we use the following elementary Lemma (its proof is given in Appendix C):

Lemma 13. Let {Z t , t ≥ 0} be a Lévy process with characteristic exponent ψ(x) such that E[e ixZt ] = e -tψ(x) for x ∈ R, and let T be a rv distributed as Γ(ζ, 1), independent from {Z t , t ≥ 0}. Then the CF of Z T is given by

E[e ixZ T ] = {1 + ψ(x)} -ζ , x ≥ 0. ( 72
)
The aim is now to write φ(s) in (71) in the form of (72). We rst write Ξ ∞,s in (70) as

Ξ ∞,s = ∞ 0 ∞ 0 1 -exp < is, v > we -ρy aλe aλy dy P(W ∈ dw).
Performing a change of variable z := we -ρy (i.e. y = -1 ρ ln z w ) within the integral in y, it may be expressed as

Ξ ∞,s = ∞ 0 w 0 (1 -exp [< is, v > z]) aλ ρ z w -aλ/ρ dz z P(W ∈ dw) = ∞ 0 (1 -exp [< is, v > z]) ∞ 0 w aλ/ρ 1 [w≥z] P(W ∈ dw) aλ ρ z -aλ/ρ-1 dz = ∞ 0 (1 -exp [< is, v > z]) E W aλ/ρ 1 [W ≥z] aλ ρ z -aλ/ρ-1 dz = R (1 -exp [< is, v > z]) Π(dz),
where the measure Π(dz) on (0, +∞) is dened as (62). Finally, we get the following expression for (71):

φ(s) = 1 + ψ(< s, v >) -b/a , s ∈ R k ,
so that one deduces from Lemma 13 the convergence result in (61).

Proof of Theorem 10 in the case ρ < aλ

After a change of variable y := ty, (55) is rewritten as

ϕ t (s) = 1 + t 0 [1 -ϕ o y (s)]aλe aλ(t-y) dy -b/a , t ≥ 0, s ∈ R k .
Let us consider the normalizing function g(t) = e aλt , so that ϕ t (s/g(t)) = ϕ t (se -aλt ) = 1 +

In the following, the limit of the integral on the right-hand side of (73) is studied in the subcritical case. First, similar to (67), let

Ξ t,s := t 0 [1 -ϕ o y (se -aλt )]aλe aλ(t-y) dy. (74)
To apply the dominated convergence theorem, let us dene

Ξ t,s,y := 1 [0<y<t] [1 -ϕ o y (se -aλt )]aλe aλ(t-y) = 1 [0<y<t] E[1 -e <is,N o (y)>e -aλt ]aλe aλ(t-y) . (75) Since |1 -e <is,N o (y)>e -aλt | ≤ < |s|, N o (y) > e -aλt , (76)
We recall from [5, p.202] that the mean matrix of the multitype process N o (t) is expressed as E[N o (y)] = e Ay E(I) where the matrix A is dened in [START_REF] Athreya | Branching processes[END_REF]. For the case ρ < aλ, the integral +∞ 0 e (A-aλId)y dy is convergent because all eigenvalues of the matrix A -aλId have negative real part in the case of ρ < aλ.

ϕ t (se -aλt ) -→ 1-< is, aλ(aλId -A) -1 E(I) > -b/a , t → +∞,
for all s ∈ R k , which entails (63) with the vector γ dened as (64).

Proof of Theorem 10 in the case ρ = aλ

We consider here the renormalizing function g(t) := te ρt = te aλt . As in (66) and (67), after a change of variable y := y/t we have for all s ∈ R k ,

ϕ t (s/g(t)) = ϕ t (se -aλt /t) = 1 + t 0 [1 -ϕ o t-y (se -aλt /t)]aλe aλy dy -b/a = 1 + t 0 1 -E exp(< is, N o (t -y) > e -aλt /t) aλe aλy dy -b/a = 1 + 1 0 t 1 -E exp(< is, N o (t(1 -y)) > e -aλt /t) aλe aλty dy -b/a = {1 + Ξ t,s } -b/a , (79)
where Ξ t,s is now dened by Ξ t,s :=

In the following we shall determine the limits of Ξ (81)

Step 1: Studying the convergence of Ξ 1 t,s as t → +∞. It is readily obtainable that using the inequality [START_REF] Laskin | Fractional Poisson process[END_REF], one has for all t ≥ 0 and y ∈ (0, 1) that t 1 -exp(< is, W v > e -aλty /t) aλe aλty ≤ t < |s|, W v > (e -aλty /t)aλe aλty =< |s|, W v > aλ, which is integrable, so that for a xed y ∈ (0, 1) one has by the dominated convergence theorem that t 1 -E exp(< is, W v > e -aλty /t) aλe aλty -→ -E [< is, W v >] aλ as t → +∞. Likewise:

t 1 -E exp(< is, W v > e -aλty /t) aλe aλty ≤ E [< |s|, W v >] aλ,
is a constant, so that by the dominated convergence theorem one deduces that

lim t→+∞ Ξ 1 t,s = -E [< is, W v >] aλ = -< is, E[W ]aλv > . (82)
Step 2: Dominating Υ 2 s (t, y). In order to study lim t→+∞ Ξ 2 t,s , we again use the dominated convergence theorem. First, it can be shown that |Υ

≤ aλE[| < s, W v > |] + aλE[| < s, N o (t(1 -y)) > e -aλt(1-y) |] ≤ aλE[< |s|, W v >] + aλE[< |s|, N o (t(1 -y)) > e -aλt(1-y) ],
where the rst inequality is obtained from [START_REF] Laskin | Fractional Poisson process[END_REF] and the last equality holds because W and N o (t(1-y)) are non negative or have non negative entries. Using again the constant κ satisfying [START_REF] Mitov | Multitype branching processes with inhomogeneous Poisson immigration[END_REF] and the martingale argument one thus obtains, together with the above result, that |Υ 2 s (t, y)| is upper bounded by some constant as

|Υ 2 s (t, y)| ≤ aλE[< |s|, W v >] + aλκ < u, E (
I) >, ∀t ≥ 0, ∀y ∈ (0, 1).

Step 3: Pointwise convergence of Υ 2 s (t, y) towards 0 as t → +∞. Let y ∈ (0, 1) be xed. Since R k can be decomposed as the direct sum of Ru and (Rv) ⊥ (the orthogonal vector space of Rv for the euclidean inner product), there exists some (unique) α ∈ R and s 0 ∈ (Rv) ⊥ such that s = αu + s 0 . Since < s 0 , v >= 0, it follows that (83) is expressed as

|Υ 2 s (t, y)| ≤ aλE[| < s, W v > -< s, N o (t(1 -y)) > e -aλt(1-y) |) = aλE(|α < u, W v > -α < u, N o (t(1 -y)) > e -aλt(1-y) -< s 0 , N o (t(1 -y)) > e -aλt(1-y) |) ≤ aλE(|αW -α < u, N o (t(1 -y)) > e -aλt(1-y) |) +E(| < s 0 , N o (t(1 -y)) > e -aλt(1-y) |), (84)
where the last line follows from the triangle inequality and the fact that < u, W v >= W < u, v >= W.1 = W . We next show that both terms in the above inequality tend to zero as t → +∞. From this decomposition of s along Ru and (Rv) ⊥ , it holds that the rst term is linked to the martingale {< u, N o (t)e -ρt >, t ≥ 0} with ρ = aλ, whereas in the second term the behaviour of {< s 0 , N o (t)e -ρt >, t ≥ 0} is determined precisely thanks to the estimates given in [START_REF] Athreya | Limit theorems for multitype continuous time Markov branching processes, II. The case of an arbitrary linear functional[END_REF] and the fact that s 0 ∈ (Rv)

E(| < s 0 , N o (t(1 -y)) > e -aλt(1-y) |) = n 0 =(n 0 (1),...,n 0 (k))∈N k E( | < s 0 , N o (t(1 -y)) > e -aλt(1-y) | N o (0) = n 0 )p n 0 ≤ n 0 =(n 0 (1),...,n 0 (k))∈N k k j=1 n 0 (j)E(| < s 0 , N j (t(1 -y)) > e -aλt(1-y) |) p n 0 ≤ n 0 =(n 0 (1),...,n 0 (k))∈N k k j=1 n 0 (j) p n 0 max j=1,...,k E(| < s 0 , N j (t(1 -y)) > e -aλt(1-y) |) = k j=1 E(I(j)) max j=1,...,k E(| < s 0 , N j (t(1 -y)) > e -aλt(1-y) |) ≤ k j=1 E(I(j)) max j=1,...,k E(| < s 0 , N j (t(1 -y)) > | 2 e -2aλt(1-y) |) 1/2 (85)
where the last line is obtained thanks to the Cauchy-Schwarz inequality. From [4, Proposition 3] together with < s 0 , v >= 0, there exists some real number a(s 0 ) < ρ = aλ as well as an integer γ(s 0 ) (both depending on s 0 , see their precise denitions in [4, (9a) and (9b)]) such that one of the three following situations occurs for all j = 1, ..., k:

E < s 0 , N j (t) > 2 =    O(e 2a(s 0 )t t 2γ(s 0 ) ) if 2a(s 0 ) > ρ = aλ, O(e 2a(s 0 )t t 2γ(s 0 )+1 ) if 2a(s 0 ) = ρ = aλ, O(e ρt ) = O(e aλt ) if 2a(s 0 ) < ρ = aλ.
Here the above three cases are corresponding to [4, a), b) and c) of Proposition 3] respectively.

In all cases, since a(s 0 ) veries a(s 0 ) < ρ = aλ, one can check easily that E

[| < s 0 , N j (t) > | 2 ]e -2ρt = E[| < s 0 , N j (t) > | 2
]e -2aλt tends to 0 as t → +∞. Since Assumption (A) holds, (85) thus tends to 0 as t → +∞ (for a xed y ∈ (0, 1)). Combining all the above results, we thus prove that both terms on the right-hand side of (84) converge to 0. Therefore, it is concluded that (81) goes to zero as t → +∞ for all y ∈ (0, 1).

Step 4: Convergence of Ξ

ϕ t (se -aλt /t) -→ 1-< is, E[W ]aλv > -b/a , t → +∞.
Hence we have proved (65).

Transient expectation when k = 2

We shall hereafter consider two-type branching processes (i.e. k = 2) to study transient expectation of the number of particles at time t. Assume that the lifetime of type j particles for j = 1, 2 is exponentially distributed as E(µ i ). The branching mechanism is given by the following generating functions in (1)

h 1 (z 1 , z 2 ) = p 1 (0, 0) + p 1 (0, 1)z 2 , h 2 (z 1 , z 2 ) = p 2 (0, 0) + p 2 (1, 0)z 1 , (z 1 , z 2 ) ∈ [0, 1] 2 ,
where probabilities p 12 := p 1 (0, 1) and p 21 := p 2 (1, 0) in (0, 1] satisfy p 12 p 21 < 1, which means that type 1 particle (resp. 2) produces a type 2 (resp. 1) particle with probability p 12 (resp. p 21 ), or else dies with probability p 1 (0, 0) = 1 -p 12 (resp. p 2 (0, 0) = 1 -p 21 ). We also suppose that there is only one immigrant of type

E[N 1 (t)] = t 0 m(t) -m(t -z) + t 0 [m(t -v) -m(t -v -z)] µ 1 µ 2 p 12 p 21 1 ζ 1 (ζ 2 -ζ 1 ) e ζ 1 v + 1 ζ 2 (ζ 1 -ζ 2 ) e ζ 2 v dv µ 1 e -µ 1 z dz, t ≥ 0, (86)
where ζ 1 and ζ 2 are given by

ζ 1 := 1 2 -(µ 1 + µ 2 ) + (µ 1 -µ 2 ) 2 + 4µ 1 µ 2 p 12 p 21 , (87)
ζ 2 := 1 2 -(µ 1 + µ 2 ) -(µ 1 -µ 2 ) 2 + 4µ 1 µ 2 p 12 p 21 . (88)
Note that the expression (86) depends on the renewal function m(t), which is explicitly available in many processes. For example, m(t) = t 0 λ(s) ds when the immigration process is a NHPP with intensity λ(•) whereas m(t) = ( b a ) 1-e -aΛ t e -aΛ t when the immigration process is a GPP with parameters (a, b, λ t ). In addition to these two processes considered in this paper, we remark that (86) for the transient rst moment is also available for other non Poisson arrival processes where their renewal functions are known. Typical examples include the case when {S(t), t ≥ 0} is a fractional Poisson process with parameter β ∈ (0, 1) (where m(t) = Ct β for some constant C > 0, see [START_REF] Laskin | Fractional Poisson process[END_REF]Expression (26)]), or when the interarrival times T i -T i-1 , i ≥ 1, follow matrix exponential distributions (in which case m(t) is explicit and given by [3, Theorem 3.1]).

Proof. The key idea is to consider the successive passage times from type 2 to type 1 of the ith particle arriving at T i , i ∈ N * which is type 1. The type of particles is changing between 1 and 2 while it remains in the same type during an exponentially distributed lifetime as long as it is alive (i.e. it has not left the system). Let G i be the number of sojourn times as type 1 before dying, so that the (G i ) i∈N * are iid with distribution

G i ∼ G(1 -p 12 p 21 )
where G denotes the geometric distribution. In other words, G i = r if the ith particle survives r times as a type 1 particle then dies, or survives one last time as a type 2 particle then dies. Let us introduce the sequence (V (r) i ) r∈N representing the successive time instants of this particle (arriving at time T i ) changing back to type 1 after being type 2 prior to its death. In other words, this ith particle becomes type 1 again at times T i + V represents the rth sojourn time of type j particle for j = 1, 2 such that the rv Y (r) j,i , j = 1, 2, i ∈ N * , r ∈ N * , are independent, with distributions given by D(Y (r) j,i ) = E(µ j ). See Figure 1 for an illustration. Then, N 1 (t) has the following expression ] corresponds to the event that type 1 particle which arrived at time T i is again type 1 at time t after its rth return time. Taking the expectation in (89), interchanging the order of summation and using the independence of G i from the arrival process {S(t), t ≥ 0} and (Y

T i T i + V (1) i T i + V (2) i T i + V (r) i T i + V (r+1) i Y (
P(T i + V (r) i ≤ t < T i + V (r) i + Y (r+1) 1,i
).

Let us denote by G (r) (•) the cumulative distribution function of V (r) i , in other words, the rth convolution of the sum of two exponential variables with mean µ 1 and mean µ 2 , and also denote G (0) (ds) = δ 0 (ds). By the independence of the arrival process {S(t), t ≥ 0}, V = P N j (t -Λ -1 (Λ(t)x)) t -Λ -1 (Λ(t)x) > t t -Λ -1 (Λ(t)x) z N j (t -Λ -1 (Λ(t)x)) > 0 -→ P X > 1 1 -x z = P((1 -x)X > z), t → +∞, for a xed x ∈ (0, 1). The latter convergence along with (98) and (99) entail that the components of the integrand in (97) satisfy t P N j (t -Λ -1 (Λ(t)x)) t > z J = t t -Λ -1 (Λ(t)x) P N j (t -Λ -1 (Λ(t)x)) t > z J N j (t -Λ -1 (Λ(t)x)) > 0 × [t -Λ -1 (Λ(t)x)] P(N j (t -Λ -1 (Λ(t)x)) > 0)

-→ 1 1 -x β j P((1 -x) X > z J ), t → +∞, ∀J ⊂ {1, ..., k}.

It is important to note that from the convergence result in (99), [t -Λ -1 (Λ(t)x)] P(N o (t -Λ -1 (Λ(t)x)) > 0) is bounded uniformly in t ≥ 0 and x ∈ (0, 1) by some constant. Furthermore, (100) says that A t converges in distribution towards (1 -x)X as t → ∞, where A t is a random vector such that D(A t ) = D N j (t-Λ -1 (Λ(t)x)) t N j (t -Λ -1 (Λ(t)x)) > 0 for all t ≥ 0. We recall that X admits the expression (44), and we thus deduce that (1 -x)X is light tailed, i.e. that there exists some vector w 0 ∈ R k close to 0 with positive entries such that E(exp(< (1-x)X , w 0 >)) is nite. The fact that w 0 has positive entries together with Cherno 's bound thus yields that P N j (t -Λ -1 (Λ(t)x)) t > z J N j (t -Λ -1 (Λ(t)x)) > 0 = P(A t > z J ) ≤ P(< A t , w 0 > > < z J , w 0 >) ≤ E(exp(< A t , w 0 >))e -<z J ,w 0 > -→ E(exp(< (1 -x)X , w 0 >))e -<z J ,w 0 > , t → ∞.

Also, t t-Λ -1 (Λ(t)x) is upper bounded in t ≥ 0 by some constant that depends on x as it is convergent towards 1 1-x as t → +∞. Therefore, from the above ndings we conclude that the following function is bounded by t P N j (t -Λ -1 (Λ(t)x)) t > z J ≤ K j x e -<z J ,w 0 > , ∀J ⊂ {1, ..., k}, ∀t ≥ 0, where K j x is some constant independent from t ≥ 0 and z ∈ R * k + . Since (R * + ) Card(J) e -<z J ,w 0 > dz J is nite for all J ⊂ {1, ..., k}, one nds by the dominated convergence theorem that the integrand in (97) satises (R * + ) Card(J) l∈J [is l exp(is l z l )] t P N j (t -Λ -1 (Λ(t)x))/t > z J dz J -→ (R * + ) Card(J) l∈J [is l exp(is l z l )] 1 1 -x β j P((1-x) X > z J ) dz J , t → +∞, ∀J ⊂ {1, ..., k}

  g. [18, Proposition 8.3.2.1 p.467]), we arrive at E[λ(t)] = aλ t 0 E[λ(s)]ds+λb for all t ≥ 0, from which the expected intensity has the closed form E[λ(t)] = bλe aλt , t ≥ 0.

t 0 [ 1

 01 ϕ o y (se -aλt )]aλe aλ(t-y) dy -b/a .

( 75 )

 75 is bounded in modulus by|Ξ t,s,y | ≤ 1 [0<y<t] E[< |s|, N o (y) >]aλe -aλy ≤ E[< |s|, N o (y) >]aλe -aλy =< |s|, E[N o (y)] > aλe -aλy := Ξ * s,y .

  from r ∈ N * with V

  ) r∈N * , j = 1, 2, yield E[N 1 (t)] = ∞ r=0 B r (p 12 p 21 ) r ,

  tv) -m(tvz)]G (r) (dv)µ 1 e -µ 1 z dz.

  It then follows from (90) that E[N 1 (t)] is given by E[N 1 (t)] = t 0 t 0 [m(tv) -m(tvz)]Ψ(dv)µ 1 e -µ 1 z dz, t ≥ 0,(91)

  nonnegative random variables X and Y . Using these results and the martingale property of {< u, N o (t)e -ρt >, t ≥ 0}, also denoting |s| := (|s 1 |, ..., |s k |), we conclude that the integrand is bounded as

  1,t and R 2,t separately as t → +∞. Limit of R 1,t in (31) as t → +∞. It will be shown that R 1,t tends to 0 as t → +∞. Let us rst note that for all x ∈ R we have |e ix -1 -ix| ≤ |x| 2 if |x| ≤ 1, and (using

  aλ y e aΛy dy = e aΛt -1, one nds from (57) and (58) that e aΛt -(e aΛt -1) f t (s) =

			-b/a	.
	Using	t 0 t	t
		aλ y e aΛy dy + 1 -	aλ y e aΛy ϕ o t-y (s)dy
		0	0
		= 1 +

aΛt -(e aΛt -1) f t (s) t 0 [1 -ϕ o

t-y (s)]aλ y e aΛy dy.

  In turn, one concludes that |1 -e <is,N o (y)>e -aλt |e aλt is upper bounded by < |s|, N o (y) > which has a nite expectation. Finally, because of the bound for the integrand Ξ t,s,y obtained in (77) and the pointwise limit in (78), one deduces that (74) Consequently, it follows that (73) converges to

	for a xed y ∈ (0, +∞) one nds	+∞ 0	Ξ * s,y dy converges. Besides,
	E[1 -e <is,N o (y)>e -aλt ]e aλt -→ -E[< is, N o (y) >],	t → +∞	(78)
	by the dominated convergence theorem. Indeed, from (76) converges to
	+∞	+∞
	Ξ t,s -→ -	E(< is, N o (y) >)aλe -aλy dy = -	< is, e Ay E(I) > aλe -aλy dy
	0	0
		+∞
	= -< is,	aλe (A-aλId)y dy E(I) >=< is, aλ(aλId -A) -1 E(I) >,
		0

as t → +∞.

  1 t,s and Ξ 2 t,s separately as t → +∞. For notational convenience, let Ξ 2 t,s := := t E exp(< is, W v > e -aλty /t) -E exp(< is, N o (t(1 -y)) > e -aλt /t) aλe aλty .

	1 0 Υ 2 s (t, y)dy where
	Υ 2 s (t, y)

  ⊥ . Indeed, one has from (94) that D(t) = O(e 2ρt ),resulting in E[||N o (t)|| 2 e -2ρt ] t≥0 is uniformly upper bounded with ρ = aλ here. Since E[| < u, N o (t) > e -aλt | 2 ] is upper bounded by E[||N o (t)|| 2 e -2ρt] up to a constant for all t ≥ 0, one deduces that the martingale {< u, N o (t)e -aλt >, t ≥ 0} is uniformly square integrable, hence converges in mean square towards W as t → +∞; and in turn, the rst term on the right-hand side of (84) converges to 0 as t → +∞. Concerning the second term, we have, using the representation (3) and triangular inequality, that

  1 entering the system at each time T i , i ∈ N, i.e. that I ∼ δ (1,0) . We shall see in Remark 15 that this assumption can be relaxed to a general incoming immigration vector I in[START_REF] Altman | On stochastic recursive equations and innite server queues[END_REF]. Also note that in this model, each particle produces one ospring only, and only one ospring of the other type. Finally, we denote by m(t) = E[S(t)] for t ≥ 0 the renewal function associated to the immigration process {S(t), t ≥ 0} with the convention m(t) = 0 when t < 0.

Theorem 14. At time t, the transient expectation E[N 1 (t)] for type 1 particle is given by

  Fig 1. Evolution of ith particle.

		1) 1,i	Y	(1) 2,i			Y	(2) 1,i	Y	(2) 2,i	Y	(r+1) 1,i	Y	(r+1) 2,i
	Arrival of ith particle
					Type 1				Type 2
								∞	G i -1
				N 1 (t) =	i=1	r=0	1 [T i +V (r) i	≤t<T i +V i (r)	+Y 1,i (r+1)	] ,	(89)
	as [T i + V i (r)	≤ t < T i + V	(r) i	(r+1) 1,i + Y

where Ψ(ds) is a measure dened by Ψ(ds) = ∞ r=0 (p 12 p 21 ) r G (r) (ds), (92) which remains to be determined. Since G (r) (ds) is the distribution of the sum of two independent Erlang distributions with respective parameters (r, µ 1 ) and (r, µ 2 ), its Laplace transform is given by

, r ≥ 0, x ≥ 0, so that, taking the Laplace transform on both sides of (92), one obtains

where ζ 1 and ζ 2 are dened by ( 87) and (88). Inverting (93) then yields

e ζ 2 s ds, s ≥ 0 which, by substituting in (91), yields (86).

Remark 15. Similar analysis is available to obtain a transient expression for E[N 2 (t)] for type 2 particles, as well as the expected number of particles, say E[M 1 (t)] of type 1 when immigration has the distribution I ∼ δ (0,1) , i.e. when one particle of type 2 arrives at each instant T i . By the superposition principle (3)), one deduces that, for a general immigration vector I = (I(1), I(2)) in (2), the total number of expected particles of type 1 is then given by

Appendix A: A growth rate for the second order moment of N o (t)

) be the second order matrix at time t ≥ 0 of the the baseline multitype branching process {N o (t), t ≥ 0}, with branching mechanism in [START_REF] Allison | Estimation and testing for a Markov model of reinforcement[END_REF] and immigration vector at time 0 in (2) described in Section 2. The following lemma is a consequence of the growth rates for the second moments provided in [START_REF] Athreya | Branching processes[END_REF]Section 7.4 p.203] which are the results obtained when N o (0) is deterministic. Lemma 16. The following growth rate holds entrywise under Assumption (A) as t → ∞:

We study the terms on the right-hand side of (95) in the following cases for the eigenvalue ρ, yielding the growth rates (94).

If ρ < 0, we have from [5, (19) p.204] that E(N j (t)N j (t) ) = O(e ρt ) entrywise and uniformly in all j = 1, ..., k, and from [5, (17)

meaning that each entry of D (n 0 ) (t) grows at most as k j=1 n 0 (j) 2 Ce ρt for some constant C independent from n 0 . This nally implies, thanks to Assumption (A), that

proving (94) when ρ < 0. If ρ = 0, we have from [5, (19) p.204] that E(N j (t)N j (t) ) = O(t) entrywise and uniformly in all j = 1, ..., k, and from [5, (17) p.203] that E(N j (t)) = O(1) entrywise and uniformly in all j = 1, ..., k. Thus similar analysis using (95) and Assumption (A) implies the growth rate (94) when ρ = 0. Finally, if ρ = 0, we have from [5, (19) p.204] that E(N j (t)N j (t) ) = O(e 2ρt ) entrywise and uniformly in all j = 1, ..., k, and from [5, (17) p.203] that E(N j (t)) = O(e ρt ) entrywise and uniformly in all j = 1, ..., k, implying in turn that E(N j (t))E(N j (t)) and E(N j (t))E(N p (t)) are O(e 2ρt ). Similar analysis implies the growth rate (94) when ρ > 0.

Appendix B: Proof of Lemma 8

Before proceeding with the proof, we recall the multidimensional version of Polya's theorem, which will be used later on. Lemma 17. Let {X t , t ≥ 0} be a sequence of random variables with values in R k converging in distribution towards X ∈ R k , such that x ∈ R k → P(X ≤ x) is continuous. Then one has for all x ∈ R k that lim t→+∞ P(X t > x t ) = P(X > x),

where lim t→∞ x t = x, x t lying in R k , and '≤' and '>' are understood componentwise.

We now turn to the proof of Lemma 8. First, for the process N j (t) = (N j 1 (t), ..., N j k (t))

described in Section 2 one has that

for l = 1, ..., k. Together with an expansion formula, we get that

where J⊂{1,...,k} is the sum over nonempty sets J ⊂ {1, ..., k}. The results in [START_REF] Weiner | On a multi-type critical age-dependent branching process[END_REF] will be repeatedly used in the following leading to the convergence (47). Taking expectations on both sides of (96), it follows that ψ j (t, x) in (50) may be expressed as

which, using the notation z J introduced just before Lemma 8 and multiplied by t, drives a more compact form given by

where, for two vectors v 1 and v 2 , v 1 > v 2 means that each entry of v 1 is larger than the corresponding one in v 2 , and dz J := l∈J dz l . Next, let us observe that, for a xed x ∈ (0, 1), from (42) it follows that t -

Since the above result entails that t -Λ -1 (Λ(t)x) -→ +∞ as t → +∞, from [33, Theorems 1 and 5], we nd for x ∈ (0, 1) that

as t → +∞ and for all z = (z 1 , ..., z k ) ∈ R * k + , where we recall that β j is given by [START_REF] Feller | On a general class of contagious distributions[END_REF] and that X has a distribution given by (45). Here again, the relation '>' is understood entrywise.

It is noted that (100) simply states that the distribution of

given that N j (t -Λ -1 (Λ(t)x)) > 0 converges to the distribution of X . Also, since z ∈ R k → P(X > z) is continuous (extending the denition in (45) from z ∈ R * k + to z ∈ R k by putting P(X > z) = 1 if max i=1,...,k z i ≤ 0), and lim t→+∞ t t-Λ -1 (Λ(t)x) = 1 1-x from (98), one has from Lemma 17 for all z that

for a xed x ∈ (0, 1). Putting this into (97) yields the following limit as t → ∞ t[ψ j (t, x) -1] -→ J⊂{1,...,k} (R * + ) Card(J) l∈J [is l exp(is l z l )] 1 1 -x β j P((1 -x) X > z J ) dz J .

By an expansion formula similar to (96) (with (1 -x)X here instead of N j (t -Λ -1 (Λ(t)x))/t), one can check that the right-hand side of the above limit is equal to


Vivien Robert

Hugues Talbot

DOES SUPER-RESOLUTION IMPROVE OCR PERFORMANCE IN THE REAL WORLD? A CASE STUDY ON IMAGES OF RECEIPTS

Keywords: Single image super-resolution, deep learning, task-based performance evaluation, image restoration

Recently, many deep learning methods have been used to handle single image super-resolution (SISR) tasks and often achieve state-of-the-art performance. From a visual point of view, the results look convincing. Yet, does it mean that those techniques are reliable and robust enough to be implemented in real business cases to enhance the performance of other computer vision tasks? In this article, we investigate the use of SISR to construct higher-resolution images of real receipt photos sent by a company's customers and evaluate its impact on the performance of an OCR task (receipt information retrieval). Using built-in task-based performance evaluation methods, we show that the use of SISR can significantly improve OCR performance in the case where recognition was poor in low-resolution, but can also deteriorate the performance for receipts that were already successfully recognized. As a conclusion, we provide recommendations on how to best use SISR in a production environment.

INTRODUCTION

When dealing with images in real applications, acquisition procedures may sometime result in low-resolution (LR) images, which are an imperfect rendering of reality. Single-Image Super-Resolution (SISR) aims at constructing higherresolution (HR) images from observed LR images which should be as similar as possible to reality. Over the last few years, AI researchers have started to use powerful deep learning (DL) algorithms for SISR. Many different methods have been used, ranging from Convolutional Neural Networks (CNN) [START_REF] He | Learning a deep convolutional network for image super-resolution[END_REF][START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF][START_REF] Kim | Accurate image super-resolution using very deep convolutional networks[END_REF][START_REF] Shi | Real-time single image and video super-resolution using an efficient sub-pixel convolutional neural network[END_REF] to Generative Adversarial Networks (GAN) [START_REF] Ledig | Photo-realistic single image super-resolution using a generative adversarial network[END_REF][START_REF] Wang | Esrgan: Enhanced super-resolution generative adversarial networks[END_REF][START_REF] Zhang | Self-attention generative adversarial networks[END_REF]. There are some issues with the current approaches to SISR. First, these algorithms mainly focus on improving the visual and perceptual image quality. From a human visual point of view, the generated images look cleaner and do feature greater resolution. Another problem is that SISR is usually trained and evaluated on synthetic datasets, where the LR images are deduced from the HR ones by down-sampling, and sometime adding various amounts of artificial noise. As pointed out recently in [START_REF] Cai | Toward real-world single image superresolution: A new benchmark and a new model[END_REF], this is not realistic, and current SISR algorithmic performances do not translate to real-world applications. The authors of [START_REF] Cai | Toward real-world single image superresolution: A new benchmark and a new model[END_REF] have proposed to build LR/HR pairs by acquiring images with different cameras and various levels of focal lengths.

Here, we evaluate the effectiveness of SISR by a different mean, specifically, can we improve the performance of other computer vision tasks by increasing the resolution of the input images with a SISR model? It is not at all obvious as the super-resolved images are not always similar to the true HR images. Indeed, they can suffer from the creation of unwanted patterns and hallucinating artifacts [START_REF] Culley | Quantitative mapping and minimization of super-resolution optical imaging artifacts[END_REF], which can lead to misinterpretation and errors. While the benefits of adopting a SISR pre-processing step have been demonstrated on four popular vision tasks for LR images, it has up to now not been shown to provide benefits on a real vision application [START_REF] Dai | Is image super-resolution helpful for other vision tasks?[END_REF]. In this article, we provide a real production use case in which SISR is applied to improve a computer vision process and evaluate the usefulness of this approach to meet a company business need.

The use case is based on a loyalty program called Scanobar and launched by the company Heineken. The main idea of this program is to reward customers when they buy Heineken products in supermarkets or bars via the Scanobar application. The latter is a chatbot available on Facebook Messenger. Customers send receipt pictures containing Heineken products to the Messenger chatbot. A computer vision algorithm detects the products and related prices on the receipts and credits a certain amount of loyalty points proportional to the price of the Heineken purchases. Google Vision's OCR is then used to extract the information (products and prices) on the received receipts. This extraction is a critical step to award the right number of loyalty points to the customers. Thus, the readability of the receipt images should be good enough so that the OCR software can read the information properly. Yet, the readability of the receipt images can be poor for 3 main reasons:

1. There is no uniformity in the pictures dataset as image sizes range from less than 5 kiB up to over 3 MiB. This can lead to information loss, especially when one needs to resize all the images.

2. The photographs are not professionally taken or framed.

They can be blurred due to motion, lack of illumination or malfunctioning autofocus, with numerous artifacts such as shadows or flash reflections. As it is typically difficult to take the whole receipt in the same picture because it is often too big to fit. As a result, the photographs is taken at a large distance from the receipts, decreasing readability significantly.

3. Finally, the photos are often heavily compressed. The pictures may be downgraded automatically when downloaded from the Scanobar bot during peak hours, resulting in much lower picture resolutions.

While there are many different causes behind the poor quality of the received pictures, in many cases, the pictures remain readable from a human point of view. However, in about 30 percent of the cases, the receipt images are readable for a human, but the OCR fails to extract the relevant information.

In this work, we develop a task-based SISR solution for the Scanobar program. It aims at improving the OCR product and price detection by enhancing the resolution of the input images which have been degraded by the 3 causes described above. We assume that the images which are not readable for humans are out of scope as they are simply too degraded. Thus, the real target is to increase the resolution of the images that are readable for humans but not for the OCR. The question is then: does the SISR model help OCR recognize information on those receipt images?

MODEL

Even though many state-of-the art super-resolution results on natural images have been produced using GANs, in our work we do not use such a solution due to a lack of data. Moreover, a GAN model may lack stability and can be highly sensitive to minor hyper-parameter changes as well as degradation models. From a production perspective, one needs a fast and robust model capable of generating simple patterns (figures and characters). This is why we use a simpler model based on two learning strategies, content and texture losses, also called perceptual loss in the original paper [START_REF] Johnson | Perceptual losses for real-time style transfer and superresolution[END_REF]. Its main objective is to enhance the images by "decrappifying" [START_REF] Antic | Decrappification, deoldification, and super resolution[END_REF] them, meaning that the algorithm removes noises and increases resolution. It corresponds well to the Scanobar use case where multiple image degradation types exist. As described in [START_REF] Antic | Decrappification, deoldification, and super resolution[END_REF], there are two main steps in the algorithm:

1. The generative model (gen) is a U-NET with an encoder based on a Resnet 34 backbone architecture which has been pre-trained on ImageNet. Blur is used in this model to avoid checkerboard artifacts at each layer. Cross-connection with the direct input of the model are applied (skip connection). This model takes degraded images as input and cleans them. The degraded images are constructed from HR images to which a custom degradation function D is applied.

2. For the learning strategy, one uses a VGG-16 classification model pre-trained on ImageNet. This loss network is used as inference to compare the generated images ŷ and its corresponding target y. The idea is to be sure that the created images have the same "concept" and "style" (in our case, mainly figures and characters) than the real images. The loss function L model is divided into 3 parts, each with its specific goal:

• To measure the overall proximity between the 2 images, an 1 pixel loss pixel is used. • To capture feature closeness, a content loss content is expressed. We hook the activation layers of the VGG 16 intermediate blocks noted J (blocks 2 to 4) for the created image ŷ and its corresponding target image y. Then, 1 differences are computed for each activation. We only use blocks 2 to 4 because these blocks are not specific while later blocks are too focused on the ImageNet classification task. • Finally, for the same feature maps J, texture loss texture is used to capture style proximity. Gram Matrices are built on the feature maps J and their 1 distances are computed.

To conclude, the final loss is a combination of the 1 loss from the generative model gen and the 1 losses from the VGG activation layers J computed each time in two different ways (content and textures losses).

L model = gen pixel (y, ŷ) + vgg16,J content (y, ŷ) + vgg16,J texture (y, ŷ)

Dataset

The training set should be as close as possible to the use case. Thus, and only for training, We used HR images of receipts which we then degraded to create the input images to our model. Fortunately, we have access to a sufficient number of real receipt images sent by customers since the launch of the application in early 2019. We combine complementary approaches:

1. We selected manually 1000 receipt images that were readable for a human and then degraded them through a custom degradation function designed to bring as close as possible to use cases (see next section). Even with data augmentation, this dataset proved too small to yield satisfying results.

2. We then added to this first dataset a set of 1000 receipt images from the ICDAR 2019 competition [START_REF] Huang | Icdar 2019 robust reading challenge on scanned receipts ocr and information extraction[END_REF]. These receipts were much too clean to correspond exactly to our business case, but they helped with training once degraded.

3. We selected a second subset of 6000 receipt images from the Scanobar database but without any visual check on the readability of the images. In spite of this, training was still improved.

Thus, in total, the final training set is composed of 8,000 (LR, HR) pairs images to which we further apply data augmentation (random rotation, random zoom and random symmetric warp).

Degradation function

The custom degradation function applied to the images should be as similar as possible to the real degradation described in introduction: 1) we shrunk the size of the input image with bilinear interpolation, then 2) we compressed the image using the JPEG lossy standard coding; 3) finally, we applied a Gaussian Blur.

The combination of the three degradation functions corresponds to the small sensor size, default compression and poor optics of many smartphone cameras used in challenging lighting and shooting distance circumstances. In our test, a resizing to 500 × 500, a quality of 15 and a blur radius of 1 were the most effective.

Normalization

Normalization is important as the receipt images are heterogeneous in terms of size, brightness and contrast. In our tests, we scaled the input with a normalization based on ImageNet statistics.

RESULTS

Task-based evaluation method

SISR models can be evaluated using simulations, subjective testing or task-based testing.

Simulations are the most common way SISR methods are evaluated. A HR image database is degraded in some way, then enhanced via super-resolution, and the results are compared with the original HR images in terms of PSNR or SSIM.

Subjective testing is also often used. Done correctly, this can be a reliable way of evaluating image quality. Its main principle is to ask to a group of human testers their opinion on the quality of various images. Several classical methods exist to perform these studies [START_REF]Methodology for the subjective assessment of the quality of television pictures[END_REF][START_REF] Mohammadi | Subjective and objective quality assessment of image: A survey[END_REF]. A challenge with subjective testing is to enroll a sufficient number of fair testers, as it can be time consuming, tiring and expensive.

Finally, task-based testing consists of using the SISR model to improve another computer vision task, which can in turn be evaluated objectively. The prediction of the computer vision task from original images and from SISR improved images are compared to measure the impact of SISR on the overall prediction performance. Task-based quality assessment exploit image attributes that are important for the computer vision task, not for humans. For instance, an object recognition model may focus on high-level semantics while ignoring the image contrast and noise which are two important features for humans [START_REF] Wang | Deep learning for image super-resolution: A survey[END_REF]. Therefore, task-based evaluation may be useful for computer vision domain-specific applications, but less so if the end usage is to improve human perceptual quality. In our production usage case, task-based evaluation is preferable since the goal of the SISR model is to enhance OCR performance.

Best model results and performance evaluation

When we use our SISR model on test images (HR images which are degraded by our custom degradation function), we can see a clear visual improvement. More importantly, when applying our model to real, non-degraded, case images, we also observe a significant increase in quality. For example, in Fig. 2, the recognized text in the low-resolution image is entirely incorrect.

In contrast, OCR on the super-resolved text is correct in 5 out of 7 lines. It reads: "6 TOMATE GRAPPE FRAN 6 MELON PIECE ; 7 6X25CL GOUDALE PRO ; 6 OIR.LIME 6000 IMP ; 7 75CL VOP CHARDONNA ; 6 SUCRINE -1EURD ; 6 INNOCENT PUR. JUS" However, sometimes SISR does not help. In Fig. 3, OCR recognizes correctly the first two lines from the real LR image (top). The super-resolved image make things worse, none of the lines in the SISR image is recognized correctly from the bottom image.

To conduct the task-based evaluation of our SISR model, we need to manually label a sufficient number of receipt im- ages from the test set. We collect the following information: food store name, products, date and postal code. Performance is evaluated on 2 different datasets from the Scanobar database and is measured using the percentages of success on the total products, store names, postal codes and ticket dates additionally detected.

The first dataset (Dataset 1 in Table 1) is composed of 100 Scanobar labeled receipt images for which OCR performed poorly and failed to detect the relevant information. Results are illustrated in Table 1 and show a significant improvement for all types of information. For instance, with the SISR pre-processing step, 10% of the total Heineken products of Dataset 1 are detected additionally by the OCR. Thus, the use of the SISR on previously non-readable images (for Google Vision OCR) facilitates OCR, as illustrated on Fig. 2.

The second dataset (Dataset 2 in Table 1) contains 100 LR labeled receipts images at least partially successfully rec- ognized. Results in Table 1 show a generally smaller decrease in the text detected by the OCR for all types of information. Thus, we observe that SISR can sometimes worsen results. To conclude, the SISR model is generally useful in our use case as it significantly increases the OCR performance, especially with noisy images. Even though OCR can in certain cases create pattern distortions, the gain in performance remains positive.

CONCLUSION

In this paper, we sought to evaluate the influence of Single-Image Super-Resolution (SISR) on OCR performance in challenging circumstances, i.e. with handheld, telephonequality photographs of poorly printed text. We showed that text recognition increased by up to 15%. On the other hand, on some images where the recognition is already good, using SISR can degrade results by up to 9%. Thus, for production implementations, we advise to run OCR both with and without the SISR model. By keeping only the most plausible result (i.e. with correct date, real product names, etc) from the two runs, one can take advantage of the SISR ability to restore noisy receipt images while not degrading performance on clean images.

Fig. 1 .

 1 Fig. 1. System overview

Fig. 2 .

 2 Fig. 2. (a) Zoom on a real raw image from the Scanobar database for which OCR fails to recognize the 7 item references. (b) Corresponding prediction made by the SISR model for which Google Vision OCR is far more accurate and manages to recognize perfectly 5 out of the 7 item references.

Fig. 3 .

 3 Fig. 3. (a) Zoom on a real raw image from the Scanobar database. Google Vision OCR detects perfectly the first two lines. (b) Corresponding prediction made by the SISR model. Google Vision OCR is less accurate and hardly detects any item reference.

Table 1 .

 1 Impact of the SISR pre-processing step on the OCR detection performance. With SISR,

	Information type	Dataset 1	Dataset 2
	Heineken Products	+10%	-5%
	Receipt Date	+5%	-1%
	Postal Code	+15%	-9%
	Food Store Name	+7%	-6%