Jean-Yves Fortin
email: jean-yves.fortin@univ-lorraine.fr

Origin of the spontaneous oscillations in a simplified coagulation-fragmentation system driven by a source

Keywords: numbers: 05.20.Dd, 05.45.-a, 36.40.Qv, 36.40.Sx

We consider a system of aggregated clusters of particles, subjected to coagulation and fragmentation processes with mass dependent rates. Each monomer particle can aggregate with larger clusters, and each cluster can fragment into individual monomers with a rate directly proportional to the aggregation rate. The dynamics of the cluster densities is governed by a set of Smoluchowski equations, and we consider the addition of a source of monomers at constant rate. The whole dynamics can be reduced to solving a unique non-linear differential equation which displays self-oscillations in a specific range of parameters, and for a number of distinct clusters in the system large enough. This collective phenomenon is due to the presence of a fluctuating damping coefficient and is closely related to the Liénard self-oscillation mechanism observed in a more general class of physical systems such as the van der Pol oscillator.

Introduction

The processes of aggregation and fragmentation are found in a large variety of systems at different levels of scale. We can mention various physical phenomena that fit into this class of problems, such as foam dynamics with coalescence of soap bubbles, the stability of atmospheric particles of different sizes which aggregate due to the van der Waals forces, or even internet large scale networks where nodes can acquire additional links according to their attachment preference weight [START_REF] Krapivsky | A Kinetic View of Statistical Physics[END_REF]. On larger scales, it is known that aggregation plays an important role in the formation of planetary rings [START_REF] Cuzzi | [END_REF] or dust clouds induced by gravitational forces, and in the formation of galaxy clusters [3,4]. In contrast, fragmentation acts as a counterweight to the aggregation phenomenon, leading in many cases to a stationary state for the mass or size distribution of the particle clusters. The mechanisms of fragmentation have diverse origins. For example, the shattering of heavy particles that collide leads to smaller particles, or tidal forces between large masses can break them apart because of the inhomogeneity of the gravitation between different parts of the same body. Thermal fluctuations are also responsible for the desegregation of large polymers in solutions.

A general kinetic theory that incorporates both the aggregation and fragmentation is usually constructed from a set of Smoluchowski type equations that govern the dynamics of the mass densities of each individual object. In the long time limit, the steady state distribution exhibits in many models, according to the type of kernel chosen for the aggregation or fragmentation processes, a power law behavior associated with an exponential decay [5,6,[START_REF] Da Costa | Mathematical Aspects of Coagulation-Fragmentation Equations Math[END_REF][START_REF] Krapivsky | [END_REF]9]. A summary of several analytical solutions for the density distribution is presented in Table 1, for different kinds of kernels and eventually in the presence of a particle source. Apart from these time independent steady states, leading to an equilibrium distribution of masses, one can also observe collective oscillations in a narrow window of parameter range [10,11,12], see Table 2. However it is not yet clear how these oscillations arise from the complex dynamics involving different clusters of particles, since it is a highly non-linear process and the Smoluchowski equations are usually quadratic in the densities and exact solutions can not be found in general, except for small sets of equations [START_REF] Buchstaber | Algebraically integrable quadratic dynamical systems[END_REF][START_REF] Calogero | [END_REF], and this is an ongoing challenge in the field of mathematics.

In this paper we would like to investigate analytically the nature and origin of these oscillations in a simplified model of aggregation and fragmentation process. Here we consider a set of N ≫ 1 possible distinct masses M k = kM 1 , k ≤ N, with their density n k , and we restrict in this manuscript the dynamics to the aggregation process between the monomer particles and bigger masses M 1 + M k → M k+1 , and the fragmentation of large masses into the smallest units M 1 + M k → (k + 1)M 1 , see figure 1. We also add a source of monomers M 1 at constant rate and eventually a sink for the largest clusters M N , which conserves the total mass of the system.

The manuscript, which tries to detail the dynamics of this model, is organized as follow. In section 2, we describe the general model and write the set of Smoluchowski master equations for the dynamics of each cluster density. These equations can be reduced to a set of self-consistent integral equations, depending solely on the solution for the monomer density which satisfies a non-linear integral equation of Volterra type. The structure of these integrals is dependent on the eigenstates of the main coefficient matrix with eigenvalues in the complex plane. As a first simple example, we demonstrate in section 3 that the monomer density can only be positive and develops a complex time signal. In section 4, we consider the approximation case of a single harmonic oscillator with negative damping for representing the eigenstate whose eigenvalue has the largest positive real part and we rewrite the Volterra equation for the monomer density as a second order non-linear ordinary differential equation (ODE). The monomer density reaches either a constant value in the long time limit or displays steady oscillations, depending on the input parameters of the model. The frequency of these oscillations is independent on the current and can be computed perturbatively around the oscillations onset. In section 5, we extend this method to the main model in order to obtain a more precise value of the oscillation frequency, and in the conclusion we explain why this oscillating phenomenon is related to the Liénard self-oscillations observed in non-linear mechanical systems with a negative damping.

k . θ = ν -µ ≥ 0 and β = µ + ν.
The mass density is normalized to unity. Kernels K ij , F ij , F k , and current I k Asymptotic distribution References

K ij = 1, constant kernel, F ij = 0 n k (t) ≃ 4t -2 e -2k/t
[15]

K ij = (i + j)/2, additive kernel, F ij = 0 n k (t) ≃ e -2t g(ke -t ), g(z) = e -z/2 (2πz 3 ) -1/2
[15]

K ij = ij, multiplicative kernel, F ij = 0 n k (t) ≃ (2π) -1/2 k -5/2 t -1
[15]

K ij = ij, F ij = 0, F k = νk, ν ≥ 2 (gelation for ν ≤ 2) n k ≃ (ν-2)k+(1+α) 2 √ π e -k log(ν/2) k -7/2
[16]

K ij = i µ j ν + i ν j µ , F ij = 0, I k = Iδ k1 , θ < 1 and β < 1 n k ≃ I(1-θ 2 ) cos(πθ/2) 4π 1/2 k -(3+β)/2
[17, 18]

K ij = i µ j ν + i ν j µ , F ij = λK ij , µ = ν = β/2 n k ≃ n 1 (1+λ) 2 √ π(1+2λ) e -k k -(3+β)/2
[11]

K ij = 1, F ij = λ, F k = ν n k ≃ c(4π) -1/2 e -νk/2 k -3/2 , c = (λ -ν + η)/(1 + 2λ) η = (λ -ν) 2 + 2ν(1 + 2λ) [9] K ij = (ij) µ , F ij = λ, F k = 0 n k ≃ λ 1-2µ Γ(1/2-µ) e -λ 2 k k -3/2-µ [9] K 1i = K i1 = 1, F 1i = F i1 = λ, λ > λ c n k ≃ λ 1+λ e -k log(1+λ)
[8]

K 1i = K i1 = i, F 1i = F i1 = λi, λ > λ c n k ≃ λe -k log(1+λ) k -1 [8]
K ij = i µ j ν + i ν j µ , I k = Iδ k1 , mass cutoff M ≫ 1 |θ| = |ν -µ| > 1 [10, 11] K ij = i -a j a + i a j -a , F ij = λK ij , I k = 0 Brownian kernel a = ν = -µ, a ≃ 1 -, λ ≪ 1 [12] K ij = i µ j ν + i ν j µ , F ij = λK ij , I k = 0 1 < θ = ν -µ ≃ 2 -(1/2 < a ≃ 1 -) and λ < λ c (θ) ≪ 1 [12]

Main model

We consider a set of clusters containing k particles, 1 ≤ k ≤ N, each with a density n k . Only one-particle clusters react with other clusters of size k ≥ 1, see figure 1. A particle can aggregate with a cluster of size k to form a cluster of size k + 1 with probability K k . Or it can shatter this cluster entirely into k single particles with probability λK k , as considered in previous works [12]. Usually it is assumed that λ ≪ 1. The rates K k grow according to the cluster size k, as a power law with exponent a, such that K k ∝ k a . A source of monomer particles is finally introduced with rate µ, and a sink term with rate -µ/N for the N-particle clusters in order to conserve the total mass. The Smoluchowski equations therefore read

∂n 1 ∂t = -2n 2 1 - N -1 k=2 (1 -kλ)K k n 1 n k + NλK N n 1 n N + µ, ∂n k ∂t = K k-1 n 1 n k-1 -(1 + λ)K k n 1 n k , 2 ≤ k ≤ N -1 ∂n N ∂t = K N -1 n 1 n N -1 -λK N n 1 n N -µ/N. (1)
This set of equations satisfies the balance relation for the mass current ∂ t k≥1 kn k = 0. It is convenient to redefine the cluster densities by introducing the rescaled variables v k = K k n k and a new definition for the time τ = τ (t) such as dτ = n 1 dt. We obtain in this case a series of quasi-linear ODEs

∂v 1 ∂τ = -2v 1 - N -1 k=2 (1 -kλ)v k + Nλv N + µ/v 1 , K -1 k ∂v k ∂τ = v k-1 -(1 + λ)v k , 2 ≤ k ≤ N -1 K -1 N ∂v N ∂τ = v N -1 -λv N -µ(Nv 1 ) -1 .
(2)

Algebraic analysis

In order to solve formally the system (2), we rewrite it in the form

∂ τ v = Av + j, with v = (v 1 , • • • , v N ),
A the N ×N coefficient matrix of the system, and

j = µ(1/v 1 , 0, • • • , 0, -K N (Nv 1 ) -1 )
Figure 1: Aggregation and fragmentation process for the main model. A one-particle collides and aggregates with a k-particle cluster of size k with probability K k , or fragment it into k individual particles with probability λK k .

the current vector. The matrix A is given explicitly by

A =             -2 -1 + 2λ -1 + 3λ • • • -1 + (N -1)λ Nλ K 2 -(1 + λ)K 2 0 0 • • • 0 0 K 3 -(1 + λ)K 3 0 • • • 0 . . . 0 K 4 -(1 + λ)K 4 0 0 0 0 0 . . . . . . 0 0 0 • • • 0 K N -1 -(1 + λ)K N -1 0 0 0 • • • 0 K N -λK N             . (3)
The general ODE for v(τ ) can be solved formally using the diagonalization of A. This leads directly to the integral solution

v(τ ) = τ 0 P e D(τ -τ ′ ) P -1 j(τ ′ )dτ ′ + P e Dτ P -1 v(0), (4)
where A = P DP -1 , D is the diagonal matrix containing the complex eigenvalues D = diag(z 1 , • • • , z N ), and P the matrix of the eigenvectors. As initial condition, we take a system consisting of only monomer particles, v(0) = (v 1 (0), 0, • • • , 0), with v 1 (0) = n 1 (0) = M 0 and M 0 being the total initial mass that we can take to be unity for simplification. The eigenvalues z i are roots of the polynomial equation

(z + 2) N -1 k=2 (1 + λ + zK -1 k ) + N -1 k=2 (1 -kλ) N -1 l=k+1 (1 + λ + zK -1 l ) = λN λ + zK -1 N , (5)
with the convention that when k = N -1 in the sum, the product inside the sum is equal to one. We can check that the eigenvalue z = 0 is always solution of this polynomial, and corresponds physically to a steady state. In the limit N ≫ 1, with K k = k a , λ > 0 and a > 0, the product

N -1 k=2 (1 + λ + zK -1 k
) is generally diverging and the secular equation ( 5) tends to the implicit series

1 = ∞ k=1 (λk -1) k l=1 1 + λ + zK -1 l -1 . (6)
This equation has real solutions near the poles z = -(1 + λ)K k < 0, as well as the steady state solution z = 0. However it is difficult to extract complex solutions in this limit, even numerically. For the simplest case when N is finite and a = 0, or K k = 1, the equation ( 5) has N -1 solutions located on the unit circle of center -1 -λ, with the value z = -λ excluded, in addition to the trivial solution z = 0

z k = -1 -λ + e 2iπk/N , k = 1, • • • , N -1, and z 0 = 0. ( 7
)
All these solutions except the zero solution have strictly negative real part. To obtain complex solutions with a strictly positive part, a has to be positive, and N sufficiently large. In figure 2 are represented the complex eigenvalues for different sizes N. The set of solutions deviates from the circle [START_REF] Da Costa | Mathematical Aspects of Coagulation-Fragmentation Equations Math[END_REF], and one observes that two complex conjugated eigenvalues converge to a finite value with a positive part in the asymptotic limit, see the inset of figure 2. In particular, this occurs when N > 10. When N = 10, all the non-zero eigenvalues have negative real parts, and in this case we do not observe any oscillatory behavior, as we will see in the following. A necessary condition, but a priori not sufficient, to observe oscillations is to have at least one pair of eigenvalues with a strictly positive real part, which occurs already for N = 20 in figure 2. A more precise analysis will be made in section 5. Each component of the vector v solution of the system (4), and satisfying the previous initial condition, can be expressed as Figure 2: Set of eigenvalues of the matrix A, with N = 10, 20, 50, 100, 200, exponent a = 0.9, and fragmentation factor λ = 0.04. Each eigenvalue comes in conjugate pair in the complex plane, and the zero eigenvalue is always solution of the polynomial (5). In inset are represented the eigenvalues close to the origin, which converge as N becomes large. Only two of them have a strictly positive real part.

v k (τ ) = µ τ 0 Q k (τ -τ ′ ) v 1 (τ ′ ) dτ ′ + Qk (τ ), (8)
with the kernel functions Q k and Qk defined by the eigenvectors of the matrix

A Q k (τ ) = N l=1 P kl (P -1 ) l1 - K N N (P -1 ) lN e z l τ , Qk (τ ) = N l=1 P kl (P -1 ) l1 e z l τ , (9)
and initial conditions

Q k (0) = δ k,1 -(K N /N)δ k,N , Qk (0) = δ k,1 . The behavior of Q k (τ ) or Qk (τ )
is dominated after some transit time by the set of conjugate eigenvalues with the largest real part.

If this one has strictly positive real and non zero imaginary parts, Q k and Qk will oscillate and diverge exponentially. Even if the Q k s and Qk s are exponentially diverging with time, we will show that the density solutions are finite. The above equations are all determined by v 1 (τ ), which satisfies the following self-consistent functional equation

v 1 (τ ) = µ τ 0 Q 1 (τ -τ ′ ) v 1 (τ ′ ) dτ ′ + Q1 (τ ), ( 10
)
and which is usually called a Volterra integral equation of the second kind in the Hammerstein form [START_REF] Burton | Volterra Integral and Differential Equations[END_REF]. It is in general difficult to obtain an exact solution for any arbitrary function Q 1 as the equation ( 10) is explicitly a non-linear form and generally non-local. The Laplace transform can be applied to (10) since it involves a convolution product, however this does not give any more information as there is no useful relations between the Laplace transforms of the function v 1 (τ ) and its reciprocal v 1 -1 (τ ).

Steady state solution

Before analyzing the time evolution of the equation ∂ τ v = Av + j, we can search for constant solutions in the long time limit. Since A is singular, solving directly (2) with ∂ τ v i = 0 does not lead to any definite answer. We consider instead the solution of the equation Av = 0, which is a one-dimensional vector, and after some computation we find that the kernel ker(A) is generated by the vector

v = v 1 1, (1 + λ) -1 , (1 + λ) -2 , • • • , (1 + λ) 2-N , λ -1 (1 + λ) 2-N ⊤ , (11)
where v 1 is any arbitrary constant. It is easy to show that if we add in the last component of this vector the term -µ(Nλv 1 ) -1 , then the new vector satisfies Av + j = 0. In order to determine v 1 , we use the mass conservation equation N k=1 kv k /K k = 1 and obtain a second order equation for

v 1 N -1 k=1 kv 1 K k (1 + λ) k-1 + Nv 1 K N λ(1 + λ) N -2 - µ K N λv 1 = 1. ( 12
)
This equation always possesses only one positive solution v * 1 , and in particular, for N ≫ 1, we obtain that v * 1 of µ in this limit. The densities n k then satisfy the equilibrium values

n k = k -a (1 + λ) k-1 v * 1 ≃ e -λ(k-1) k a v * 1 . (13)

Numerical analysis

Numerically, we can solve the system of equations ( 2) with a 4 th order Runge-Kutta algorithm for N different densities v k as a function of the effective time τ . In the following we consider only the time τ since the real time t is an increasing and monotonic function of τ . Indeed n 1 is found to be always strictly positive. In figure 3 is represented the first 10 densities v k=1,•••,10 when N = 100 for the set of parameters a = 0.9, λ = 0.04, and µ = 0.002. They all display collective oscillations, with a retarded phase shift between each consecutive curve. For this number N, the largest pair of eigenvalues of the matrix A is found to be approximately equal to z ± max ≃ 0.169 ± 1.531i, and has therefore a strictly positive real part. All the solutions of system ( 2) are oscillating at a frequency Ω estimated after Fourier transform of the v 1 signal to be equal to Ω ≃ 1.226, which is lower than the imaginary part Im(z + max ). We therefore can not relate directly the observed frequency to Im(z + max ), which is a natural frequency of the problem, and a more detailed analysis has to be performed in the next sections.

Simplified model

In order to investigate the previous ODE (10), and show in particular that v 1 is always a positive quantity with time, we consider a very simplified kernel version of Q 1 (τ ) and Q1 (τ ) which is represented by the general sawtooth or triangular function with equidistant time steps τ n≥0 such that T = τ n+1 -τ n with τ 0 = 0, and increasing slope coefficients α n > 0

Q 1 (τ ) = 1 -α 0 τ, 0 ≤ τ ≤ T Q 1 (τ ) = Q 1 (τ n ) + (-1) n+1 α n (τ -τ n ), τ n ≤ τ ≤ τ n+1 , (14)
where the extrema are given by Q

1 (τ n ) = 1 + T n-1 k=0 (-1) k+1 α k .
Replacing Q 1 and Q1 in equation ( 10) by ( 14), and differentiating with respect to τ , we obtain the following first order and non-local integro-differential equations for each time interval

v ′ 1 (τ ) = µ v 1 (τ ) -µα 0 τ τ -τ 1 dτ ′ v 1 (τ ′ ) + µα 1 τ -τ 1 τ -τ 2 dτ ′ v 1 (τ ′ ) -• • • + (-1) n+1 µα n τ -τn 0 dτ ′ v 1 (τ ′ ) + (-1) n+1 α n , τ n ≤ τ ≤ τ n+1 . ( 15
)
Deriving a second time with respect to τ , we obtain a set of piecewise second order ODEs

v ′′ 1 (τ ) = - µv ′ 1 (τ ) v 2 1 (τ ) - µα 0 v 1 (τ ) + µ(α 0 + α 1 ) v 1 (τ -τ 1 ) - µ(α 1 + α 2 ) v 1 (τ -τ 2 ) + • • • + (-1) n+1 µ(α n-1 + α n ) v 1 (τ -τ n ) , τ n ≤ τ ≤ τ n+1 . (16)
Each solution on the intervals τ n ≤ τ ≤ τ n+1 depends on the previous solutions v 1 (τ -τ k ) for k = 0, • • • , n and can be recursively computed from the first solution in the interval 0 ≤ τ ≤ T which satisfies

v 2 1 (τ )v ′′ 1 (τ ) = -µ (v ′ 1 (τ ) + α 0 v 1 (τ )) . ( 17
)
Initial conditions are v 1 (0) = 1 and v ′ 1 (0) = µ -α 0 ‡. Solving this system of ODEs requires the continuity of v 1 at the interval boundaries:

v 1 (τ + 1 ) = v 1 (τ - 1 )
. There also exists a relation between the first derivatives at the boundaries v ′ 1 (τ + 1 ) = v ′ 1 (τ - 1 ) + α 0 + α 1 , which can be derived directly from equation (15). We use the following generalized conditions at τ = τ n :

v 1 (τ + n ) = v 1 (τ - n ), v ′ 1 (τ + n ) = v ′ 1 (τ - n ) + (-1) n+1 (α n-1 + α n ). ( 18
)
‡ We could also rescale the time variable s = α 0 τ with ṽ1 (s) = v 1 (τ ) to express (17) with only one parameter κ = µ/α 0 , such that in the first time interval ṽ2 1 (s)ṽ ′′ 1 (s) = -κ(ṽ ′ 1 (s) + ṽ1 (s)).

Solution in the time interval 0 ≤ τ ≤ T

In the first interval 0 ≤ τ ≤ T , we assume that the parameter κ = µ/α 0 is small. The initial conditions are given by v 1 (0) = 1 and v ′ 1 (0) = µ -α 0 < 0, as seen previously. We expect in this case that v 1 (τ ) decreases initially with Q 1 (τ ). The overall exact solution is given by an implicit relation between τ and v 1 . We first assume that v 1 can be expressed implicitly through the functional integral

τ = 1 v 1 (τ ) v Z(v) dv, (19)
where Z(v) is an unknown function to be determined. By differentiating the previous expression with respect to τ , we have the identity 17) as functions of Z(v 1 (τ )) and v 1 (τ ), we obtain a new ODE of reduced order

Z(v 1 (τ )) = -v 1 (τ )v ′ 1 (τ ) > 0. Expressing v ′ 1 (τ ) and v ′′ 1 (τ ) in the equation (
Z ′ (v) - Z(v) v = µ 1 v - α 0 v Z(v) , (20)
where Z ′ (v) is the derivative of Z with respect to v and not τ . The initial condition at τ = 0 is replaced by Z(1) = α 0 -µ > 0, according to the previous values of v 1 (0) and v ′ 1 (0). We then perform the substitution W (v) = (Z(v) + µ)/v in order to simplify even further the previous ODE

W ′ (v) = µα 0 µ -vW (v) . (21)
This equation belongs to the class of Abel ODEs of the second kind [START_REF] Polyanin | Handbook of exact solutions for ordinary differential equations[END_REF][START_REF] Murphy | Ordinary Differential Equations and Their Solutions[END_REF]. A way to resolve it is to use an inverse transformation, solving v as a function of W instead [START_REF] Cheb-Terrab | [END_REF]. This yields the linear equation

µα 0 v ′ (W ) = µ -W v(W ), (22)
whose solutions v(W ) are simply given by the integrating factor

v exp W 2 2µα 0 -α -1 0 W 0 exp u 2 2µα 0 du = C, (23)
which determines W and Z as a function of v implicitly. C is a constant which can be evaluated by considering the identity W (v 1 (0)) = α 0 at τ = 0. This gives the following value of C and its expansion for κ = µ/α 0 small

C = e 1/2κ - 1 0 exp u 2 2κ du ≃ (1 -κ)e 1/2κ . (24)
A similar expansion for (23) can be performed when κ is small, assuming that µ is small compared to

W v - µα 0 W (v) exp W 2 (v) 2µα 0 ≃ C. (25)
We can then obtain a further approximate value for Z(v), assuming that Z(v) always stays positive

Z(v) ≃ v 2µα 0 log (C/v) -µ. ( 26
)
This implies finally that the density v 1 (τ ) is given implicitly by

τ ≃ 1 v 1 (τ ) vdv v -2µα 0 log (e -1/2κ v) -µ . ( 27
)
In equation [START_REF] Burton | Volterra Integral and Differential Equations[END_REF], we have assumed that v 1 (τ ) ≤ 1 is decreasing and that the function Z(v) has to be strictly positive, otherwise τ may not be defined if Z vanishes and becomes negative. For small values of τ , we can approximate (27) with Z(v) ≃ α 0 v and therefore obtain that v 1 (τ ) ≃ 1 -α 0 τ , which follows the initial behavior of the function Q 1 and decreases as τ increases. However, for τ large, the inspection of equation ( 27) implies that v 1 (τ ) should reach asymptotically from above a limiting value v * > 0 determined by the singularity of Z, and solution of Z(v * ) = 0:

v * -2µα 0 log (e -1/2κ v * ) = µ. (28)
In the limit of κ small, we obtain that v * ≃ κ + κ 2 log κ. As τ increases, the integrand becomes indeed larger when v approaches v * from above. This is equivalent to implying that v 1 (τ ) → v * > 0 when τ ≫ 1, with T or α 0 large.

Numerical solution in the other intervals

τ n ≤ τ ≤ τ n+1
We study here the global solutions ( 16) for the successive time intervals, and for different sawtooth geometries such as the slope of each linear parts, T = 1 being fixed. In regard to the time behavior of (9), we consider a diverging amplitude for the oscillating function Q 1 (τ ), by taking α n = e βn α 0 , with β > 0. Figure 4 displays the solutions for different values of α 0 with β = 0.1. We have checked that using numerical algorithms for solving directly the non-local integral equation (10), or solving the ODE (16) with the Runge-Kutta-Fehlberg RKF45 method via Maple, leads to the same results. For small α 0 , v 1 follows the sawtooth behavior of Q 1 but stays positive. When α 0 increases, v 1 develops into a more complex periodic signal. It also takes very small values at times τ k = (2k + 1)T , and its maximum amplitude is also smaller as additional frequency modes appear in the signal. It is more difficult to obtain solutions when α 0 > 2, as v 1 takes very small values in time intervals that cannot be resolved precisely with the algorithms. In summary, v 1 is periodic and stays positive when we apply the approximate sawtooth kernel function Q 1 which is a diverging oscillating function with large negative and positive values. However additional modes appear in the time Fourier spectrum as α 0 increases, and which seem to be harmonics of the fundamental frequency π/T .

Model with a driven harmonic oscillator and negative damping parameter

In this section, we replace the kernels Q 1 (τ ) and Q1 (τ ) by a function representing an oscillator characterized by a negative damping parameter. The idea is that to approximate Q 1 (τ ) and Q1 (τ ) defined in equation ( 9) by a simple cosine function of frequency ω multiplied by a diverging exponential factor at rate γ. The quantities ω and γ can be viewed as the imaginary and real parts of the largest eigenvalue z + max . Although this approximation is correct asymptotically, it is not the case at small times since all the other eigenstates will give a finite contribution. However we expect to extract enough information about the physics of the system. We also introduce a phase shift ϕ ∈] -π/2, π/2[ such that

Q 1 (τ ) = e γτ cos(ωτ -ϕ) cos ϕ , cos ϕ = 0, (29)
with initial condition Q 1 (0) = 1, independent of the phase. The phase diagram is actually symmetrical by the reflection operation ϕ → -ϕ, so that we will restrict afterwards the analysis to the case when ϕ ≥ 0. The derivative of this function at the origin is equal to Q ′ 1 (0) = γ + ω tan ϕ > 0. Substituting (29) into equation (10), we obtain

v 1 (τ ) = µ cos ϕ τ 0 e γ(τ -τ ′ ) cos(ω(τ -τ ′ ) -ϕ) v 1 (τ ′ ) dτ ′ + e γτ cos(ωτ -ϕ) cos ϕ . ( 30
)
We can use the fact that Q ′′ 1 -2γQ ′ 1 + (γ 2 + ω 2 )Q 1 = 0 to eliminate the dependence in Q 1 by considering the first and second derivative of the previous expression. After combining the different terms we finally obtain a second order ODE independent of the non-local integral

v ′′ 1 (τ ) -2γv ′ 1 (τ ) + (γ 2 + ω 2 )v 1 (τ ) = -µ v ′ 1 (τ ) v 2 1 (τ ) -(γ -ω tan ϕ) µ v 1 (τ ) . ( 31
)
The initial conditions are v 1 (0) = 1 and v ′ 1 (0) = γ + ω tan ϕ + µ. Notice that we can remove the dependence in µ by substituting v 1 → √ µv 1 . The left hand side is the linear ODE of an harmonic oscillator that Q 1 (τ ) satisfies, with frequency ω and negative damping parameter -γ, and the right hand side are non-linear terms. This equation has a time-independent solution in the long time limit provided that tan ϕ > γ/ω and which is equal to

v * 1 = µ ω tan ϕ -γ γ 2 + ω 2 1/2 . ( 32
)
In the opposite case, tan ϕ ≤ γ/ω, the solution v 1 (τ ) decays to zero exponentially with τ , since for small v 1 the two terms on the right hand side of equation ( 31) are dominant, and we have instead

v ′ 1 (τ ) ≃ -(γ -ω tan ϕ)v 1 (τ ). (33)
This yields v 1 (τ ) ∼ e -(γ-ω tan ϕ)τ → 0. To study the stability of the solution (32), we introduce a correction w(τ ) to the steady state solution v 1 (τ ) = v * 1 + w(τ ), and linearize the equation (31)

w ′′ (τ ) -2γ - µ v * 1 2 w ′ (τ ) + 2(γ 2 + ω 2 )w(τ ) = 0. ( 34
)
The solution w will grow if the damping factor is positive, or 2γ ≥ µ/v * 1 2 . This will occur if the phase ϕ is larger than a critical value ϕ c which satisfies

tan ϕ ≥ 3γ 2 + ω 2 2γω = tan ϕ c > γ ω . (35)
In this case, we observe steady oscillations with a frequency Ω determined approximately at this order by

Ω 2 ≃ 2γ 2 + 2ω 2 , ( 36
)
and which is independent of the current µ as expected. In table 3 we have summarized the phase diagram as function of ϕ, according to the properties of (34).

Table 3 0 ≤ tan ϕ ≤ γ/ω v 1 (τ ) goes to zero exponentially (33) γ/ω < tan ϕ ≤ (3γ 2 + ω 2 )/2γω v 1 (τ ) goes to the equilibrium value v * 1 (32) (3γ 2 + ω 2 )/2γω < tan ϕ v 1 (τ ) oscillates with a finite amplitude and frequency Ω (36)

In figure 5(a) is represented numerically the time evolution of v 1 as function of different values of the phase ϕ. There are two domains according to whether ϕ is lower or greater than ϕ c , for which v 1 reaches asymptotically the equilibrium value v * 1 or self-oscillates respectively. The oscillation frequency Ω can be extracted by Fourier transform of the time signals and represented as a curve function of ϕ for different values of ω, see figure 5

(b).

There is no simple relation between Ω and the other parameters, except at the critical angle ϕ c for which Ω = 2γ 2 + 2ω 2 , see equation (36).

Near ϕ c , we can assume that the oscillation amplitude is small compare to the constant solution v * 1 . The previous linear expansion (34) shows that the damping factor 2γ -µ/v * 1 2 must be positive in order to generate oscillations. However, this expansion leads to an unstable solution as w diverges with time. We can study the next terms of the expansion and verify whether these terms stabilize the oscillation amplitude. If we expand (31) up to the third order in w, we obtain a non-linear equation with a time dependent damping factor

w ′′ (τ ) -2γ - µ v * 1 2 + 2 µ v * 1 3 w(τ ) -3 µ v * 1 4 w 2 (τ ) w ′ (τ ) +(2γ 2 + 2ω 2 )w(τ ) + (γ -ω tan ϕ) µ v * 1 3 w 2 (τ ) -(γ -ω tan ϕ) µ v * 1 4 w 3 (τ ) = 0. ( 37
)
We have checked that this ODE displays stable oscillations close to the exact solution for ϕ ϕ c , see figure 6(a). These oscillations are not symmetric, and their amplitude depends on ω, see figure 6(b). For example, as ω decreases, the amplitude increases. Restricting the expansion to the second order would still lead to non stable oscillations, due in particular to the presence of the linear term w(t) in the damping coefficient, in factor of w ′ (t). The stability comes instead from the quadratic term w 2 (t) in the same damping coefficient and from the last term in w 3 (t). To analyze the stability of (37), we consider an expansion in the small parameter 0 < tan ϕ-tan ϕ c ≃ δϕ ≪ 1.

The technique we use in the following is described in the reference book [START_REF] Nayfeh | Nonlinear oscillations Wiley Classics Library[END_REF] and can be applied to different ODEs with negative damping containing a linear or quadratic amplitude in addition to the constant term. This is usually useful to study self-excited oscillators. We first simplify the previous ODE in order to reduce the number of parameters. We chose to rescale the amplitude

w(τ ) = c w(τ ), with c = v * 1 ωγδϕ/(γ 2 + ω 2 ), so that w′′ (τ ) + Ω 2 0 w(τ ) = 4γǫ 1 + w(τ ) - 3ǫ 2 w2 (τ ) w′ (τ ) + ǫ 2 Ω 2 0 w2 (τ ) -ǫ w3 (τ ) , ( 38
)
with Ω 2 0 = 2(γ 2 + ω 2 ), and ǫ = ωγδϕ/(γ 2 + ω 2 ) ≪ 1 the rescaled small parameter. This yields in particular c = ǫv * 1 . When ǫ = 0, we can write the general solution (38) as w(τ ) = u cos(Ω 0 τ + β). When ǫ = 0, the general solution can still be expressed in the same form, provided the coefficients u and β are time dependent. Since we have two functions satisfying one ODE, we can chose a constraint between u and β. It is usually convenient to impose the identity w′ (τ ) = -uΩ 0 sin(Ω 0 τ + β). This implies in particular that

u ′ (τ ) cos φ -uβ ′ (τ ) sin φ = 0, ( 39
)
where we used the notation φ = Ω 0 τ + β(τ ). Equation (38) can then be formally rewritten as w′′ (τ ) + Ω 2 0 w(τ ) = ǫf ( w, w′ ) = ǫf (u cos φ, -uΩ 0 sin φ),

with function f defined by

f ( w, w′ ) = 4γ 1 + w - 3ǫ 2 w2 w′ + 1 2 Ω 2 0 w2 -ǫ w3 . ( 41
)
After computing the second derivative w ′′ (τ ) from the expression w ′ (τ ) = -uΩ 0 sin φ, we obtain

w′′ (τ ) + Ω 0 w(τ ) = -u ′ (τ )Ω 0 sin φ -uΩ 0 β ′ (τ ) cos φ. ( 42
)
The right hand side can be identified with (40), and using the relation between u and β, see (39), we obtain the following set of first order ODEs

u ′ = - ǫ Ω 0 sin φf (u cos φ, -uΩ 0 sin φ), β ′ = - ǫ Ω 0 u cos φf (u cos φ, -uΩ 0 sin φ). ( 43
)
These equations are identical to (38), since no approximation have been made yet. We then assume that φ varies rapidly compared to the variations of u and β. We can replace in this approximation the previous equations by their time average over a period 2π/Ω 0 , after considering the quantities u ′ and β ′ quasi constant. This yields

u ′ ≃ - ǫ 2πΩ 0 2π 0 dφ sin φf (u cos φ, -uΩ 0 sin φ), β ′ ≃ - ǫ 2πΩ 0 u 2π 0 dφ cos φf (u cos φ, -uΩ 0 sin φ). (44)
After performing the integrations over φ, we obtain a set of coupled ODEs for the amplitude and phase shift

u ′ = γǫ(2u - 3ǫ 4 u 3 ), β ′ = 3 16 Ω 0 ǫ 2 u 2 . (45)
We should notice that the evolution of the amplitude u does not depend on the terms in factor of Ω 2 0 in equation (41), after averaging over φ, but is determined by the damping terms in factor of w′ , especially the constant and quadratic term w2 . However, the change in frequency depends instead on the cubic term in factor of Ω 2 0 . All the terms which are a odd function of φ do not contribute to (45). The presence of the quadratic term in the damping factor leads to the existence of a finite equilibrium limit for u which converges exponentially to a positive constant u 0 determined by the non-zero solution of u ′ = 0 in (45), or explicitly u 0 = 8/(3ǫ). This leads, after rescaling to the original variable, to the equilibrium amplitude w 0 = 8ǫ/3v * 1 , and therefore we can conclude that the effect of the additional quadratic term in the damping factor is to stabilize the oscillations that are generated self-consistently. The frequency Ω is modified accordingly after solving (45) for β ′ in the asymptotic limit, and is shifted by a small amount such that the corrected frequency is given by Ω ≃ Ω 0 + β ′ or Ω ≃ Ω 0 (1 + ǫ/2). It is, as expected, independent of µ, contrary to the amplitude which is proportional to v * 1 ∝ √ µ. As an example, if we take the parameters from figure 5(b), ω = 3, γ = 1 and µ = 0.04, we obtain for ϕ = 1.17 ϕ c ≃ 1.107 the overall estimated amplitude of the signal 2w 0 ≃ 0.185, which is close to the numerical value extracted from the analysis of the ODE (38) or 0.217, with ǫ = 0.128. We also obtain Ω ≃ 4.759, to be compare with Ω 0 ≃ 4.472 and the numerical estimation from the Fourier transform of the signal Ω ≃ 4.634. The present model in this section does not consider all the eigenstates with eigenvalues that have negative real part, since we have assumed that these states have a negligible contribution as they correspond to the terms in (9) that vanish exponentially with time. However their contribution at small times may be relevant to the general solution as the convolution integral equations ( 8) and ( 10) are non local. In the following section, we will generalize the previous method to the main model which incorporates N oscillators.

Generalization to a set of harmonic oscillators

The generalization of the previous results to the main model of section 2 is important in order to analyze the onset of the self-oscillations, and to extract the value of the frequency observed for example in figure 3 and which can not be approximated with a simple oscillator. We first write a general ODE for v 1 by removing all dependence in the Q 1 and Q1 functions, as well as the convolution integrals. We expect that for N oscillators, v 1 satisfies a N th -order ODE. As in the previous section, we can write an ODE satisfied by Q 1 and Q1 , since they are a sum of N harmonic oscillators (or N -1 precisely because of the presence of the zero eigenvalue eigenstate). From equation (9), it is easy to establish a linear ODE satisfied by Q 1 and Q1 using the characteristic polynomial of the coefficient matrix A whose coefficients are given by det (zI

N -A) = N k=0 a k z k , (46)
with a N = 1 and a 0 = 0. Since the eigenvalues are roots of this polynomial, we have directly an ODE for Q 1 and Q1 after differentiating ( 9)

N k≥1 a k Q (k) 1 (τ ) = 0. ( 47
)
By differentiating N times the identity (10), and multiplying each successive equations by a k , we can eliminate all the time dependence involving Q 1 and Q1 , and therefore all the convolution integrals disappear. This yields

N k=1 a k v (k) 1 -µ k-1 l=0 Q (k-1-l) 1 (0) ∂ l ∂τ l 1 v 1 = 0. ( 48
)
This equation contains only derivatives in v 1 since a 0 = 0 and also because there is no contribution in v 1 coming from the l = 0 terms in the previous sum as This comes from the fact that one of the factors α l in the sum (9), Q 1 (τ ) = N l=1 α l e z l τ , vanishes when l = l 0 , with l 0 associated to the zero eigenvalue z l 0 = 0. This has been checked numerically. Otherwise, it is easy to see that the previous equation (48) would not be consistent with the existence of the time-independent solution that we found in section 2.2. There is therefore no v 1 term in the ODE, contrary to the case of the single harmonic oscillator that we studied in the previous section. Each derivative term of the ODE (48) can be rearranged by increasing differential order of v -1

N k=1 a k Q (k-1) 1 (0) = 0. (49)
1 N k=1 a k v (k) 1 -µ 1 v 1 ′ a 2 Q 1 (0) + a 3 Q ′ 1 (0) + • • • + a N Q (N -2) 1 (0) -µ 1 v 1 ′′ a 3 Q 1 (0) + a 3 Q ′ 1 (0) + • • • + a N Q (N -3) 1 (0) -• • • -µ 1 v 1 (N -1) a N Q 1 (0) = 0. ( 50
)
We can define (N -1) additional coefficients b

k = N -k-1 l=0 a k+1+l Q (l) 1 (0) such that the previous ODE becomes N k=1 a k v (k) 1 -µ N -1 k=1 b k 1 v 1 (k) = 0. ( 51
)
We can, as before, eventually remove the dependence on µ by rescaling v 1 → √ µv 1 . The complex coefficients b k can be computed recursively as the inverse powers of the eigenvalues, using the factors α l . For the first two terms, we have indeed

b 1 = -a 1 l =l 0 α l z -1 l , b 2 = -a 1 l =l 0 α l z -2 l -a 2 l =l 0 α l z -1 l , (52)
and for the general case b k = -k i=1 a i l =l 0 α l z -k-1+i l . From these results, we can analyze the stability of the steady state solution v * 1 of the equation ( 12) if we linearize (51) by setting

v 1 = v * 1 + w, |w| ≪ v * 1 , which yields N k=1 a k w (k) + µ v * 1 2 N -1 k=1 b k w (k) = 0. ( 53
)
The stability of the solution can be studied by finding the roots of the following polynomial

R(z) = N -1 k=1 a k + µ v * 1 2 b k z k + a N z N . ( 54
)
This polynomial is different from the characteristic polynomial of A due to the additional terms b k . The perturbation w(τ ) will grow if at least one complex root of R(z) has a strictly positive real part. In figure 7(a) we have represented the roots of R(z) for different N, with a fixed set of parameters. As N increases, one pair of conjugate roots moves into the positive part of the real axis, and self-oscillations occur, see figure 7(b). The imaginary part of the root z + max of R(z) with the largest real part gives an estimate of the oscillation frequency. We find for N = 50 z + max ≃ 0.092 + 1.454i, and the Fourier transform of the signal of figure 7(b) gives approximately Ω ≃ 1.332 for this value of N. Im(z + max ) is closer to Ω than the value given by the characteristic polynomial of A which is equal to Im(z + max ) ≃ 1.561. A more precise analysis of (51) requires to estimate the contributions from all the non-linear terms. The successive derivatives of v -1 1 can be ordered as function of the powers of v * 1 -1 . This can be done recursively by using a combinatorial method

N k=1 a k w (k) + µ N -1 k=1 b k 1 v * 1 2 w (k) (1 + w/v * 1 ) 2 - 1 v * 1 3 k-1 l 1 =1 k l 1 w (l 1 ) w (k-l 1 ) (1 + w/v * 1 ) 3 (55) + 1 v * 1 4 k-1 l 1 =1 k-l 1 -1 l 2 =1 k l 1 k -l 1 l 2 w (l 1 ) w (l 2 ) w (k-l 1 -l 2 ) (1 + w/v * 1 ) 4 -• • • + (-1) k+1 v * 1 k+1 k!w (1) k (1 + w/v * 1 ) k+1 = 0.
This expression is exact. This gives a formal expansion in v * 1 -1 . When v * 1 is large enough, we can approximate this ODE by retaining only the first few terms. As in section 4, we can analyze the stability of the solution associated to the pair of conjugate roots z ± max = Λ 0 ±iΩ 0 of the polynomial (54) with Λ 0 assumed in the following to be strictly positive but smaller than Ω 0 . We consider the perturbation amplitude w(τ ) = e Λ cos(Ω 0 τ + β), after rescaling w = √ µ w and set ǫ = √ µ/v * 1 which is the small parameter of the expansion. The quantities Λ = Λ 0 τ + a(τ ) and β depend on time and a perturbative solution can be given as an expansion at the lowest order in ǫ. We also choose a constraint between the two functions Λ and β and impose the following identity for the derivative: w′ (τ ) = e Λ (Λ 0 cos φ -Ω 0 sin φ), where we used the notation φ = Ω 0 τ + β(τ ). This implies the following condition

a ′ (τ ) cos φ -β ′ (τ ) sin φ = 0. ( 56
)
After some calculations and keeping only the linear terms in a ′ and β ′ , we establish the following successive derivatives of w

w(k) (τ ) ≃ e Λ [Z k (φ) + c k (a ′ (τ ) sin φ + β ′ (τ ) cos φ)] , Z k (φ) = Re (z + max ) k e iφ , c k = -Im z + max (z + max ) k-1 -Λ k-1 0 z + max -Λ 0 . ( 57
)
We then replace the derivatives of w in the ODE (55) by their expression from (57), and organize the different terms as a series expansion in ǫ. We obtain a single equation for the combination a ′ (τ ) sin φ + β ′ (τ ) cos φ:

(a ′ (τ ) sin φ + β ′ (τ ) cos φ) N -1 k=1 c k (a k + ǫ 2 b k + O(ǫ 3 )) ≃ ǫ 3 e Λ N -1 k=1 b k k-1 l=1 k l Z l Z k-l + 2Z k cos φ -ǫ 4 e 2Λ N -1 k=1 b k k-1 l 1 =1 k-l 1 -1 l 2 =1 k l 1 k -l 1 l 2 Z l 1 Z l 2 Z k-l 1 -l 2 +3 k-1 l=1 k l Z l Z k-l cos φ + 3Z k cos 2 φ + O(ǫ 5 ). (58)
As before, we can multiply this expression by sin φ and cos φ, and use the relation (56) to obtain independently an equation for a ′ and β ′ respectively. A time average over the period 2π/Ω 0 is then performed, after we assumed that a ′ and β ′ are quasi constant quantities. After some algebra, we obtain the set of ODEs

a ′ (τ ) N -1 k=1 c k (a k + ǫ 2 b k ) ≃ - 3ǫ 4 8 e 2Λ N/2 k=1 (-1) k b 2k-1 Ω 2k-1 0 β ′ (τ ) N -1 k=1 c k (a k + ǫ 2 b k ) ≃ - 3ǫ 4 8 e 2Λ (N -1)/2 k=1 (-1) k b 2k Ω 2k 0 , (59)
where we have dropped the dependence in Λ 0 ≪ Ω 0 in the coefficients in order to simplify the expressions. We can rewrite the first equation using the substitution u(τ ) = e Λ , and obtain the relation u ′ ≃ Λ 0 u -3ǫ 3 a 0 u 3 /8, with the constant a 0 given by

a 0 = N/2 k=1 (-1) k b 2k-1 Ω 2k-1 0 N -1 k=1 c k (a k + ǫ 2 b k ) . ( 60
)
u reaches an asymptotic non-zero limit if a 0 > 0, given by u 2 0 = 8Λ 0 /(3a 0 ǫ 4 ). If Λ 0 < 0, u 0 = 0 and no oscillation is observed. If Λ 0 = 0 and a 0 < 0, additional corrections are necessary in order to verify the existence of oscillations. The equation for the frequency shift can be written as β ′ ≃ -3ǫ 4 u 2 β 0 /8, with constant β 0 given by

β 0 = (N -1)/2 k=1 (-1) k b 2k Ω 2k 0 N -1 k=1 c k (a k + ǫ 2 b k ) . ( 61
)
The shift in frequency will therefore be given asymptotically by β ′ = -Λ 0 β 0 /a 0 , independent of ǫ. The sign of the shift depends on the sign of β 0 , in the case where a 0 > 0. We can extend the expansion (58) up to the ǫ 6 order, and find that β ′ ≃ -ǫ 4 u 2 β 0 (3/8 + 5ǫ 2 u 2 /16) §, which should give a negligible correction. For example, for the parameters given in figure 7 and N = 80, we obtain numerically ǫ ≃ 0.319, which is small enough to apply the previous expansion. In particular, after solving for the roots of the polynom R, we find z + max = Λ 0 + iΩ 0 ≃ 0.026 + 1.669i, and therefore Λ 0 /Ω 0 ≪ 1. The frequency Ω from the numerical resolution of the ODEs (2) was evaluated by Fourier analysis and found to be equal to Ω ≃ 1.562, while the correction Ω 0 + β ′ ≃ 1.517 gives a closer value to the numerical estimate than Ω 0 . In this case we find that a 0 > 0 and β 0 > 0. The amplitude of the oscillations, see figure 7(b), corresponding to √ µu 0 , was found after a non-linear curve fitting with a simple cosine function to be equal approximately to 0.071 while the amplitude given by the expansion (59) gives 0.047, for v * 1 ≃ 0.044, which is of the same order of magnitude. The expansion (59) gives a first estimate of the oscillation characteristics, and shows that these oscillations can be stabilized by considering the first non-linear terms of the ODE (51).

Conclusion

In this paper, we have considered a simple model of coagulation-fragmentation of particle clusters which displays collective oscillations for a distribution of masses broad enough. The massdependent kernel leads to a coefficient matrix that has at least two conjugate complex eigenvalues with strictly positive real part, provided the size of the matrix is sufficiently large. This is a necessary condition for the emergence of a self-oscillating mechanism. This model allows us to evaluate all the densities from the monomer density solution v 1 that satisfies an implicit non-local Volterra integral equation. This equation depends on two kernels or excitation functions (which are identical in the large N limit) that have a strongly exponential divergence with positive and negative values. However v 1 appears to oscillate with a finite amplitude and with a minimum positive value proportional to the current. The oscillation amplitude depends on the parameter values such as current, the exponent of the coagulation kernel and the fragmentation rate. When § The numerical coefficients 3/8 and 5/16 are the angle average of cos 4 φ and cos 6 φ respectively.

we approximate the kernels with a linear or sawtooth function, we have shown that v 1 can never be negative, and is larger than µ when µ is small in the first approximation.

When we consider a harmonic function as an approximation of the kernels in the long time limit, i.e. with a diverging oscillatory behavior, the model shows self-excited oscillations that can be explained by the presence of a fluctuating negative damping in the expansion of the corresponding ODE around the threshold beyond which oscillations appear. This phenomenon is very similar to the behavior of a van der Pol oscillator in non-conservative systems which has a quadratic contribution in the damping factor. When the amplitude increases, the damping factor tends to be positive because of the minus sign in front of the quadratic term, and therefore the amplitude is reduced. And inversely, when the amplitude is small, the damping factor is negative and tends to increase the overall amplitude. The non-conservative oscillator can sustain self-oscillations only in the presence of a current, which drives the amplitude. The period of these oscillations does not depend on the current, however no exact solution can be found for its dependence on the other parameters, except around the oscillation threshold using a perturbative calculus, where it is found to be equal in this approximate model to √ 2 times the modulus of the largest eigenvalue of the coefficient matrix. The generalization to N oscillators leads to a more complex ODE that can be solved using the same perturbative method for some range of small parameters, but it is also the presence of negative damping factors that gives rise to self-oscillations.

Several questions remain concerning the presence of self-oscillations observed in more general aggregation kernels and without source [12]. In the present manuscript, the oscillations are driven by a source of monomers, but a more comprehensive study is needed for closed and conservative systems yielding oscillations without any external source, including extended fragmentation processes between clusters of different masses.

Figure 3 :

 3 Figure 3: Numerical solution of the system (2) for the parameters N = 100, a = 0.9, λ = 0.04, and µ = 2.10 -3 . The successive oscillating densities from v 1 (τ ) (bottom) to v 10 (τ ) (top) are shown and shifted upwards for more clarity. The vertical dashed line shows the progressive phase shift between the densities.

Figure 4 :

 4 Figure 4: Solutions of the ODEs (16) for a regular sawtooth function Q 1 with α n = α 0 e βn , β = 0.1, µ = 0.04, T = 1, and for different slope parameters α 0 .

Figure 5 :

 5 Figure 5: (a) Solution of the ODE (31) for different angles ϕ. The parameters are ω = 5, γ = 1, and µ = 0.04. The critical angle is ϕ c = arctan(14/5) ≃ 1.228 > arctan(γ/ω) = arctan(1/5). (b) Oscillation frequency Ω as function of ϕ for several values of ω, with parameters γ = 1, and µ = 0.04. The dashed line is the threshold ϕ = ϕ c for which oscillations appear and which delimitates the domain of existence of Ω. This corresponds to Ω 2 = 2γ 2 + 2ω 2 . In the limit ϕ = ϕ c = π/2, ω → 0 and Ω = √ 2γ.

Figure 6 :

 6 Figure 6: (a) Approximation corresponding to (37) (red dashed line) and comparison with the exact solution (31) (black line) for ω = 5, γ = 1, and µ = 0.04. ϕ = π/2 -0.3 is chosen to be close to ϕ c : δϕ = ϕ -ϕ c ≃ 0.043. (b) Oscillation cycles for three different parameters ω, with γ = 1, ϕ = π/2 -0.1, and µ = 0.04. The dot points correspond to the equilibrium solution v * 1 .

Figure 7 :

 7 Figure 7: (a) Plot of the complex roots of R(z), equation (54), for different values of N. The parameters are µ = 2.10 -4 , λ = 0.06 and a = 0.9. In inset is represented the roots near the origin. Only the cases N = 30, N = 50, and N = 80 have a pair of roots with strictly positive real parts. (b) Time dependence of v 1 (τ ) for the previous cases. Oscillations are observed for N = 30, 50, 80. For the other cases N = 10, 20, v 1 reaches the asymptotic limits v * 1 ≃ 0.047 and v * 1 ≃ 0.045 respectively, see equation (12).

Table 1 :

 1 Asymptotic distributions for several models of aggregation-fragmentation in the steady state. The notations are the following: Aggregation kernel K ij , fragmentation kernel F ij , spontaneous fragmentation rate F k , and source of k-mass clusters I

Table 2 :

 2

Models with stable oscillations Kernels K ij , F ij , and current I k Parameter range for stable oscillations References

Acknowledgments

I would like to thank the Department of Physics and Astronomy at the Seoul National University where part of this work was carried out for its hospitality.


Gabriel Marchesan Almeida

Sameer Varyani
email: varyani@lirmm.fr

Rémi Busseuil
email: busseuil@lirmm.fr

Nicolas Hebert
email: hebert@lirmm.fr

Gilles Sassatelli
email: sassatelli@lirmm.fr

Pascal Benoit
email: benoit@lirmm.fr

Lionel Torres
email: torres@lirmm.fr

Michel Robert
email: robert@lirmm.fr

Providing Better Multi-Processor Systems-on-Chip Resources Utilization by Means of Using a Control-Loop Feedback Mechanism

Keywords: MPSoC, homogeneous, adaptive, NoC, distributed memory, RTOS

In this paper we propose a strategy for better exploiting Multi-Processor Systems-on-Chip resources utilization by means of using a control-loop feedback mechanism. We apply the proposed techniques in a purely distributed memory MPSoC architecture that is composed of a frequency scaling module responsible for tuning the frequency of processors at run-time. Results show very promising in terms of adaptation capabilities for system with dynamic workload. Performance results demonstrate the effectiveness of the proposed approach when workload requirements for applications may vary, affecting the overall performance of the system. For validating the proposed approach we have implemented a multi-thread MJPEG decoder application and created an architecture model with/without perturbations in the system.

I. INTRODUCTION

Over the years scientists in all domains have put huge efforts for better exploiting and optimizing available resources. In massively parallel systems, especially in Multi-Processor Systems-on-Chip (MPSoC) resource management is mostly modeled targeting power optimization, performance and fault-tolerance strategies in such systems.

Since this work targets massively parallel on-chip multiprocessor systems, adaptability is a major concern in the approach. In the direction of scalable systems we have developed an homogeneous MPSoC architecture with distributed memory making use of a message passing programming model. Given the large variety of possible use cases that these platforms must support and the resulting workload variability, offline approaches are no longer sufficient as application mapping paradigms, because they do not allow coping with time changing workloads. Furthermore, the large number of parameters that play a role on the platform performance makes it difficult to estimate the best system response at design time.

In [START_REF] Almeida | An adaptive message passing mpsoc framework[END_REF] we have proposed an adaptive strategy that is responsible for making decisions at run-time. Decisions are taken by processors in a distributed fashion and relate mostly to application performance. This paper goes into the direction of adaptive strategies for homogeneous MPSoC architectures. We propose a strategy for better exploiting architecture resources utilization by means of using a control-loop feedback mechanism and mostly addresses both the benefits brought by this technique as well as the associated optimization. The new architecture relies on a Real-Time Operating System (RTOS) with support for semaphores, mutexes and task priority based scheduling algorithm. Moreover, a tiny implementation of a communication stack comprising UDP TCP/IP protocols is used aimed at providing important concepts such as ports and enabling advanced features such as reordering, rerouting, error detection, etc. [START_REF] Almeida | Evaluating the impact of task migration in multi-processor systems-on-chip[END_REF].

The rest of the paper is organized as follows. The next section presents related works. Section III describes the proposed architecture. Section IV discusses the proposed controller while experimental evaluations are presented in Section V. Finally some conclusions are drawn in Section VI.

II. RELATED WORKS A. Adaptation in MPSoCs Architectures

Recently, researchers have put focus on adaptation techniques in order to handle with dynamic and unpredictable behaviors that can appear in nowadays embedded systems. This section presents some works that have been conducted in this direction using techniques as such as task migration mechanisms and dynamic voltage and frequency scaling.

1) Task Migration Mechanisms: A number of work in the literature based on distributed memory systems has been using shared memory as a mean for enabling task migration [START_REF] Bertozzi | Supporting task migration in multi-processor systems-onchip: A feasibility study[END_REF][4] [START_REF] Pittau | Impact of task migration on streaming multimedia for embedded multiprocessors: A quantitative evaluation[END_REF]. In [START_REF] Acquaviva | Assessing task migration impact on embedded soft real-time streaming multimedia applications[END_REF] each core runs a single operating system instance in its logical private memory. Processor cores execute tasks from their private memory and explicitly communicate with each other by means of shared memory. The target platform uses a shared bus as interconnect. In [START_REF] Greiner | Tsar: a scalable, shared memory, many-cores architecture with global cache coherence[END_REF] authors propose to implement a scalable shared memory many-cores architecture with global cache coherence. The architecture is built around 4096 cores which makes use of a logically shared memory but physically distributed, with cache coherence enforced by hardware, using a directorybased protocol.

In the case of distributed memory MPSoCs, in [START_REF] Almeida | An adaptive message passing mpsoc framework[END_REF] we have proposed an adaptive strategy that is responsible for making decisions at run-time. Decisions are taken by processors in a distributed fashion and relate mostly to application performance. In [START_REF] Almeida | Evaluating the impact of task migration in multi-processor systems-on-chip[END_REF] authors address both the benefits brought by task migration mechanisms as well as the associated performance penalty in a purely distributed homogeneous MPSoC architecture.

Taking into account the future homogeneous MPSoC systems, scalable architectures with purely distributed memory system are suitable. To the best of our knowledge, our architecture is the only one with a purely distributed memory system in which does not rely on using shared memory for enabling task migration [START_REF] Almeida | Evaluating the impact of task migration in multi-processor systems-on-chip[END_REF]. Instead, it uses the NoC as a communication link where the tasks are transmitted during the migration process.

2) Dynamic Voltage and Frequency Scaling (DVFS): DVFS has been a widely applied technique for reducing power consumption in many domains, especially those concerning micro-architectures such as embedded systems. In [START_REF] Ghasemazar | Minimizing the power consumption of a chip multiprocessor under an average throughput constraint[END_REF] authors show a technique for minimizing the total power consumption of a Chip Multiprocessor (CMP) while maintaining a target average throughput. The proposed solution relies on a hierarchical framework, which employs core consolidation, coarse-grain dynamic voltage and frequency scaling (DVFS). The problem of optimally assigning dependent tasks in multiprocessor system has not been addressed in their paper.

In this paper we propose a novel solution targeting both intra/extra-communication between tasks. Intracommunications are handled by means of using local FIFOs, while extra-communication makes use of a communication stack. The model is based on Message Passing Interface (MPI) and tasks communicate with each other by exchanging messages over a Network-on-Chip (NoC) [START_REF] Moraes | Hermes: an infrastructure for low area overhead packetswitching networks on chip[END_REF].

Adaptability has been explored under a number of aspects in embedded systems, ranging from adaptive modulation used in the future 3GPP-LTE standard (SDR for software defined radio) [START_REF] Clermidy | An open and reconfigurable platform for 4g telecommunication: Concepts and application[END_REF] to adaptive cores instantiation in dynamically reconfigurable FPGAs [START_REF] Puschini | Dynamic and Distributed Frequency Assignment for Energy and Latency Constrained MP-SoC[END_REF]. Two of the foremost popular techniques are here discussed. This paper presents two major contributions: 1) Propose a novel and purely distributed memory architecture with adaptation capabilities by means of using a Proportional Integral Derivative (PID) controller; 2) Analyze and discuss the benefits of using such strategy for better exploiting architecture resources; We demonstrate how soft-real time application requirements in homogeneous MPSoC architectures can be handled by better exploiting architecture resources using a PID controller and how we can obtain a model for the system that is essential to determine the PID controller parameters. The proposed strategy is validated by scenarios where processors frequency can be tuned (reducing/increasing) without violating soft-real time constraints.

III. SHOP ARCHITECTURE

The key motivations of our approach are scalability and adaptability; the system presented in the rest of this paper is built around a distributed memory/message passing system that provides efficient support for task migration. For these reasons the architecture is named SHoP (Selfadaptive Homogeneous Platform). This system aims at achieving continuous, transparent, and decentralized runtime task placement on an array of processors for optimizing application mapping according to various potentially timechanging criteria. Figure 1 presents a structural view of the SHoP architecture. The architecture is made of a homogeneous array of Processing Elements (PEs) communicating through a packet-switching network. For this reason, the PE is called Network Processing Unit (NPU). The NPU is built of two main layers, the network layer and the processing layer. The network layer is essentially a compact routing engine based on the Hamiltonian Routing Algorithm [START_REF] Lin | Deadlock-free multicast wormhole routing in 2-d mesh multicomputers[END_REF]. The Network-on-Chip (NoC) used in this work was proposed in [START_REF] Moraes | Hermes: an infrastructure for low area overhead packetswitching networks on chip[END_REF].

The processing layer is based on a simple and compact RISC microprocessor, its static memory, and a few peripherals (one timer, one interrupt controller, one UART). The processor used has a compact instruction set comparable to a MIPS-1 [START_REF]Mips technologies[END_REF]. It has 3 pipelines stages, no cache, no Memory Management Unit (MMU), and no memory protection support in order to keep it as small as possible. A multitasking tiny Real-Time Operating System (RTOS) implements the support for time-multiplexed execution of multiple tasks. The original version of the RTOS as well the RTL description of the processor used as part of this work are available at [START_REF] Rhoads | Plasma -most mips i(tm)[END_REF].

IV. PID CONTROLLER A Proportional Integral Derivative (PID) controller is a generic control loop feedback mechanism widely used in industrial control systems. A PID controller calculates an error value as the difference between a measured process variable and a desired setpoint. The controller attempts to minimize the error by adjusting the process control inputs. In the absence of knowledge of the underlying process, a PID controller is the optimal controller [START_REF] Bennett | A History of Control Engineering[END_REF]. However, for best performance, the PID parameters used in the calculation must be tuned according to the nature of the system while the design is generic, the parameters depend on the specific system.

In this paper we propose the usage of a PID controller for adjusting the appropriated frequency of the processors at the same time as deadline miss ratio is reduced. Figure 2 summarizes a traditional PID controller. The PID controller algorithm involves three-term control: the Proportional, the Integral and Derivative values. The proportional value determines the reaction to the current error, the integral value resolves the reaction based on the sum of recent errors, and the derivative value sets the reaction based on the rate at which the error has been changing [START_REF] Bennett | A History of Control Engineering[END_REF]. Heuristically, these values can be interpreted in terms of time: P depends on the present error, I on the accumulation of past errors, and D is a prediction of future errors, based on current rate of change.

The proportional, integral, and derivative terms are summed to calculate the output of the PID controller. Defining u(t) as the controller output, the final form of the PID algorithm is:

u(t) = M V (t) = K p e(t) + K i t 0 e(τ )dτ + K d d dt e(t) (1)
where: Proportional gain, K p : larger values typically mean faster response since the larger the error, the larger the proportional term compensation. An excessively large proportional gain will lead to process instability and oscillation.

Integral gain, K i : larger values imply steady state errors are eliminated more quickly. The trade-off is larger overshoot: any negative error integrated during transient response must be integrated away by positive error before reaching steady state.

Derivative gain, K d : larger values decrease overshoot, but slow down transient response and may lead to instability due to signal noise amplification in the differentiation of the error.

V. CASE STUDY AND RESULTS

A. Introduction

In order to demonstrate the efficiency of our approach we have performed our experiments based on an implemented MJPEG decoder algorithm mapped onto the SHoP platform. The application is basically composed of four main tasks: the sender (T 1 ), responsible for feeding the application with the input stream and three additional tasks which belong to the MJPEG application (IVLC (T 2 ), IQ (T 3 ) and IDCT (T 4 )). The initial application mapping is the following: T 1 is mapped onto NPU 1 while T 2 , T 3 and T 4 are hosted in NPU 2 . Figure 3 shows the MJPEG application task graph and its initial mapping.

B. System Characterization

In order to ensure optimized PID control loop, the process characterization has been performed. It has to be either linear or pseudo-linear and the dynamic response of the process has to be evaluated.

The proposed PID controller is modeled as a service provided by the RTOS and it is available in each NPU of the architecture in a distributed fashion. We will consider the process as a black box where input is the current processor frequency and output the application throughput. Aiming to validate its linearity, a scenario composed of steps in frequency, changing from 55 to 355MHz has been used. The observing window is consisted of the time interval between two frequencies changes. By applying this method we can consider the process at a steady state before performing a state changing.

Figure 4 shows a linear interpolation of the process according to the MJPEG performance vs processor frequency. Marks in blue represent obtained throughput during the time considered as steady state. The standard deviation is approximately 150KB/s, which is less than 12% of the average throughput. The linearity consideration of the proposed process can be assumed. The dynamic response of the process has been evaluated by changing the frequency from 55MHz to 355MHz. Figure 5 shows the process reactiveness when changing processor frequency. The rise time of the process is calculated in approximately 400ms (from 0.49s to 0.53s). This value can be explained by the average necessary to reject noises due to burst of incoming frames.

Indeed, as three tasks of the MJPEG are mapped onto the same NPU, the computation power utilized by T 2 , T 3 and T 4 is shared. This creates a considerable delay with no processed frames in the output.

C. Task Migration

One of the purpose of the feedback controller is to be able to adapt the processor frequency when a perturbation occurs in the system, and in particular a task migration. To prove the efficiency of the controller, a scenario using a task migration has been created. In this scenario, two tasks communicating with each other have been mapped onto NPU 3 . Later, one of those two tasks is migrated to the NPU 2 , generating an overhead in the application performance of the NPU 2 . Figure 6 shows the throughput variation due to this migration. The NPU frequency has been maintained at 55MHz. The task is migrated at 1.86s, represented by the dashed line.

1.7 The throughput goes from 3.15MB/s down to 0.45MB/s, reducing the performance of the system. This behavior can be considered as instantaneous and it is model as perturbation. For such scenarios we are capable of measuring the cost in terms of reaction time for having a dynamic response generated by the PID controller. Hence, tuning the PID controller to reject the perturbation will be the same as tuning it for the dynamic response of the process.

D. PID Controller Tuning

Using the Takahashi method [START_REF] Takahashi | Method of designing control system based on a partial knowledge of control object[END_REF], we can optimize the controller only using the dynamic response as shown in Figure 5. The principle of the Takahashi method is to use the step response of a process to compute the PID coefficient, K p , K i and K d . It takes into account the latency of the process and the rise time to compute the coefficient which minimize the overall error:

k≥0 | k |
where k is the error at sampled time k.

We demonstrate in the next section how the controller regulates the throughput of the MJPEG application in two different scenarios: in the first one the process is constant and no perturbation is modeled. In the second scenario a task migration is performed in order to generate a perturbation in the system and the behavior of the controller is analyzed.

E. Regulation Without Perturbation

The PID controller has been implemented as a task in the architecture. Thanks to an absolute timer not affected by the frequency scaling, the controller function is executed at fixed time slots. The information about packets arrival is sent to the controller by the last task of the application.

In the first scenario we have modeled a MJPEG decoder which requires a minimum throughput of 1MB/s. This specification is used as a setpoint for the PID controller. A set of 400 frames is computed by the MJPEG application and the controller is responsible for adjusting the frequency of the processor at run-time. The system has to ensure a minimum throughput of the application according to the setpoint. Figure 7 shows the frequency regulation and the obtained throughput. We can observe a big overshoot at the beginning of the regulation. This phenomenon can be explained by two factors. First, there is a burst of incoming packets induced by the internal organization of the system that reports to the controller instantaneous information of the system. Second, traditional PID controllers produce this kind of error which is minimized at the time. These two phenomena lead to a snowball effect: as no packet arrives, the controller increases the processor frequency, which results to a faster computation of frames by the IVLC and IQUANT tasks, generating a burst of packets to the input of the last task (IDCT) of the application. This overshoot should occurs only at the beginning of the computation, or in scenarios with huge variations of setpoints values, which can be considered as scarce phenomena.

F. Regulation With Perturbation

In order to validate the adaptability of our proposed strategy, a different scenario has been modeled with perturbations that are triggered at run-time. The scenario is similar to the one presented in the previous section. The difference is that a task migration to NPU 2 is performed in the middle of the regulation. Figure 8 shows the behavior of the proposed controller in a scenario with system perturbation. In this scenario, task migration is triggered at 1.9s. We can observe the efficiency of the regulation: the cost in terms of performance due to the migration process is rapidly amortized by the process: after approximately 250ms the application setpoint is again guaranteed by applying changes in the frequency. The maximum overshoot is around 330KB/s, which is coherent with the value obtained for the scenario without regulation. The different between obtained frequency and the expected frequency is due the fact that values are rounded to the closest available frequency in the architecture (55, 105, 155, 205, 255, 305 and 355MHz).

VI. CONCLUSIONS

In this paper we have proposed a strategy for better exploiting architecture resources utilization in a purely distributed memory homogeneous MPSoC architecture. Results show very promising regarding to the adaptability of the system in scenarios with dynamic workload changing.

We have demonstrated the efficiency of the proposed PID controller by presenting two different scenarios with/without perturbations. For validating our approach we have implemented a multi-thread version of the MJPEG decoder which utilizes a model based on message passing interface.

Results shown that by using the proposed strategy obtained throughput is very close to the expected throughput in scenarios where perturbations are not considered. In scenarios where perturbations are taken into account, we have shown that the system is capable of reacting to the throughput dropping by means of using a frequency scaling module proposed in the architecture.

As perspectives we plan to investigate the behavior of the system for different types of regulators i.e non-linear regulators. We also aim to scale the system by adding more NPUs and analyze the efficiency of the proposed approach for applications with different processing requirements.

Figure 1 .

 1 Figure 1. Structural View of the SHoP Architecture

Figure 2 .

 2 Figure 2. PID Controller

Figure 3 .

 3 Figure 3. MJPEG Task Graph and Initial Mapping

Figure 4 .

 4 Figure 4. MJPEG Application Performance

Figure 5 .

 5 Figure 5. System Reactiveness

THROUGHPUTFigure 6 .

 6 Figure 6. Migration Cost Without PID Controller Regulation

Figure 7 .

 7 Figure 7. Frequency Regulation