Christophe Chalons
email: christophe.chalons@uvsq.fr

Samuel Kokh
email: samuel.kokh@cea.fr

Maxime Stauffert
email: maxime.stauffert@uvsq.fr

An all-regime and well-balanced Lagrange-projection type scheme for the shallow water equations on unstructured meshes

come

An all-regime and well-balanced Lagrange-projection type scheme for the shallow water equations on unstructured meshes

Introduction

We are interested in the numerical approximation of the shallow water equations (SWE)

   ∂ t h + ∇ • (hu) = 0, ∂ t (hu) + ∇ • (hu ⊗ u) + ∇ gh 2 2 = -gh∇z, (1a) (1b)
where x ∈ R 2 → z(x) denotes a given smooth topography and g > 0 is the gravity constant. Both the water depth h and the velocity u = (u 1 , u 2 ) ∈ R 2 depend on the space and time variables, namely x ∈ R 2 and t ∈ [0, ∞). We assume that the initial water depth h(x, t = 0) = h 0 (x) and velocity u(x, t = 0) = u 0 (x) are given.

Let us briefly properties of system (1) in the case ∇z = 0 : the system is strictly hyperbolic over the phase space Ω = {(h, hu) ∈ R 3 | h > 0}. Moreover, if n ∈ R 2 is an arbitrary unit vector, the eigenstructure of (1) is composed by two genuinely nonlinear characteristic fields associated with the eigenvalues {u T n -c, u T n + c}, where c := √ gh is the sound speed, and a linearly degenerated field associated with the eigenvalue u T n. We recall also that the regions where (u T n) 2 < c 2 (resp. (u T n) 2 > c 2 ) are called subcritical or subsonic (resp. supercritical or supersonic).

We are interested in this work in developing a numerical scheme that satisfies the well-balanced property. More specifically we want our scheme to strictly preserve the "lake at rest" steady solutions, that are the states satisfying

h + z = constant, u = 0.
For a review on numerical schemes that satisfy the so-called well-balanced property we refer the reader to the pioneering work [START_REF] Bermudez | Upwind methods for hyperbolic conservation laws with source terms[END_REF], books [START_REF] Bouchut | Nonlinear stability of finite Volume Methods for hyperbolic conservation laws: And Well-Balanced schemes for sources[END_REF] and [START_REF] Gosse | Computing qualitatively correct approximations of balance laws[END_REF]. We also refer to [START_REF] Chalons | A large time-step and well-balanced Lagrange-projection type scheme for the shallow water equations[END_REF] where the authors focus on the 1D case and already propose a well-balanced Lagrange-projection strategy. At last, in [START_REF] Díaz | A fully wellbalanced Lagrange-projection-type scheme for the shallow-water equations[END_REF] the proposed Lagrange-projection scheme is exact for a full set of equilibrium solutions (and not only the lake at rest). The Lagrange-projection methodology is especially well suited for subsonic or near low-Froude number flows. We use an implicit-explicit strategy that allows to keep a stable scheme under a CFL time step limitation which is driven only by (slow) material waves and not by (fast) acoustic waves. The implicit-explicit Lagrange-projection [START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF] scheme is designed following the pioneering work [START_REF] Coquel | Entropy-satisfying relaxation method with large time-steps for Euler IBVPs[END_REF]. More recent works are concernet with the case of Euler systems in the large friction or low-Mach regimes [START_REF] Chalons | Large time step and asymptotic preserving numerical schemes for the gas dynamics equations with source terms[END_REF][START_REF] Chalons | Operator-splitting based AP schemes for the 1D and 2D gas dynamics equations with stiff sources[END_REF][START_REF] Chalons | An all-regime Lagrangeprojection like scheme for the gas dynamics equations on unstructured meshes[END_REF] for single or two-phase flow models. The treatment of the low Froude number is considered through the so called all-regime (or asymptotic-preserving) property and follows the anti-diffusive technique on the pressure numerical flux introduced in [START_REF] Dellacherie | Analysis of Godunov type schemes applied to the compressible Euler system at low Mach number[END_REF] and also used in [START_REF] Chalons | An all-regime Lagrangeprojection like scheme for the gas dynamics equations on unstructured meshes[END_REF].

The SWE has been largely studied and one can find nice overviews and references in the books [START_REF] Bouchut | Nonlinear stability of finite Volume Methods for hyperbolic conservation laws: And Well-Balanced schemes for sources[END_REF] and [START_REF] Gosse | Computing qualitatively correct approximations of balance laws[END_REF]. The scheme proposed in [START_REF] Chalons | A large time-step and well-balanced Lagrange-projection type scheme for the shallow water equations[END_REF] in one dimension has been studied in the framework of SWE and more specifically its behaviour for low-Froude number flows in [START_REF] Zakerzadeh | On the Mach-uniformity of the Lagrange-projection scheme[END_REF]. A different implicit-explicit methodology in two dimension context has been proposed by [START_REF] Bispen | Imex large time step finite volume methods for low Froude number shallow water flows[END_REF].

In section 2, we study the dimensionless system associated to the SWE (1) and its asymptotic limit in low Froude regimes. In section 3, we present the Lagrange-projection like acoustic / transport decomposition associated to system (1). In section 4, we present the schemes, the finite volume scheme in 1D, the study of its truncation error in low Froude regimes and the proposed correction, and finally the extension towards 2D schemes on unstructured meshes. At last, we show some numerical results in 2D to verify the well-balanced property and illustrate the behaviour of the scheme in different regimes, especially in the low Froude one.

Low Froude limit for continuous equations

Dimensionless shallow water equations

In this section, we briefly introduce the dimensionless SWE. These equations will be useful to study the low-Froude asymptotic behaviour of the solutions of (1). With this in mind, we define the following dimensionless quantities

t = t T , x = x L , h = h h 0 , ũ = u u 0 , z = z z 0 ,
where T , L, h 0 , u 0 and z 0 are respectively reference time, length, water height, velocity and topography such that

u 0 = L T and z 0 = h 0 .
Defining the Froude number Fr by Fr = u 0 c 0 , where c 0 = √ gh 0 is the reference sound speed, easy calculations then give the dimensionless SWE

       ∂ th + ∇ x • hũ = 0, ∂ t hũ + ∇ x • hũ ⊗ ũ + 1 Fr 2 ∇ x h2 2 = - 1 Fr 2 h∇ x z. ( 2a
) (2b)

Asymptotic equations in low Froude limit

In this section, we give the asymptotic behaviour of the solutions of the SWE equations in the low Froude limit. If we omit the tilde notation for the sake of readability in system (2) and, if we introduce the dimensionless pressure function p(h) = h 2 2 , we get

   ∂ t h + ∇ • (hu) = 0, ∂ t (hu) + ∇ • (hu ⊗ u) + 1 Fr 2 ∇p = - 1 Fr 2 h∇z. (3a) (3b)
Let us assume that h and z admit the following expansions in powers of the Froude number:

h = h (0) + h (1) Fr + h (2) Fr 2 + O(Fr 3 ) and u = u (0) + u (1) Fr + u (2) Fr 2 + O(Fr 3 ),
which gives in particular

p = p (0) + p (1) Fr + p (2) Fr 2 + O(Fr 3 ) = p(h (0) ) + h (1) p (h (0) )Fr + O(Fr 2 ).
The governing equations give at order -2 and -1 with respect to the Froude number that

∇p (0) + h (0) ∇z = 0 ⇔ ∇h (0) = -∇z ⇔ h (0) + z = H(t), ∇p (1) + h (1) ∇z = 0 ⇔ h (0) ∇h (1) = 0 ⇔ ∇h (1) = 0 ⇔ h (1) = h (1) (t).
The asymptotic behavior is then given by

∂ t h (0) + ∇ • (h (0) u (0) ) = 0, ∂ t (h (0) u (0) ) + ∇ • (h (0) u (0) ⊗ u (0) ) + ∇p (2) = -h (2) ∇z. (4a) (4b)
Now if we impose one of the following velocity boundary conditions

Ω ∇ • u dΩ = 0 and Ω ∇ • (zu) dΩ = 0 or Ω ∇ • (hu) dΩ = 0 , integrating (4a)
with respect to the space variable gives

0 = Ω ∂ t h (0) + ∇ • (h (0) u (0) ) dΩ = Ω ∂ t (H -z) dΩ + Ω ∇ • ((H -z) u (0) ) dΩ = Ω ∂ t H dΩ + H Ω ∇ • u (0) dΩ - Ω ∇ • (zu (0) ) dΩ = ∂ t Ω H dΩ = |Ω|∂ t (h (0) + z) = |Ω|∂ t h (0)
thus ∂ t h (0) = 0 and h (0) +z = H is constant both in space and time. This leads to ∇•(h (0) u (0) ) = 0 and therefore

∇ • u (0) = ∇ • ( z H u (0) ),
while the evolution of u is given by

1 - z H ∂ t u (0) + ∇ • (u (0) ⊗ u (0) ) + 1 H ∇p (2) = ∇ • ( z H u (0) ⊗ u (0) ) -h (2) ∇ z H .
Notice that when the topography is flat, i.e. z = 0, the three equations

         h (0) + z = H = cst ∇ • u (0) = ∇ • ( z H u (0) ) 1 - z H ∂ t u (0) + ∇ • (u (0) ⊗ u (0) ) + 1 H ∇p (2) = ∇ • ( z H u (0) ⊗ u (0) ) -h (2) ∇ z H
degenerate towards the incompressible Euler equations

         h (0) = cste ∇ • u (0) = 0 ∂ t u (0) + ∇ • (u (0) ⊗ u (0) ) + 1 h (0) ∇p (2) = 0.

An acoustic/transport operator decomposition

Let us first introduce notations related to our discretization. We suppose that the computational domain Ω ⊂ R 2 is covered by N polygonal cells (Ω j ) 1≤j≤N . We consider Γ, a face of the cell j, and we suppose the following admissibility assumptions are satisfied:

• either there exists a single 1 ≤ k ≤ N such that Γ = Ω j ∩ Ω k = ∅. In this case we note Γ = Γ jk and Γ jk can either be a vertex or a single face of the mesh,

• either Γ ⊂ ∂Ω and we suppose that there exists a single k > N that will help to index ghost values for boundary conditions and we shall note Γ = Γ jk .

For 1 ≤ j ≤ N , we note N (j) the set of indices k such that Γ jk is a face of Ω j and if k ∈ N (j) we set n jk to be the unit normal vector to Γ jk pointing out of Ω j .

We can now turn to the acoustic / transport decomposition of the system (1). If we develop the spatial derivatives and isolate the transport terms (u • ∇)ϕ, where ϕ = h, hu, we can use a splitting operator with respect to time to obtain on one hand the acoustic step

∂ t h + h∇ • (u) = 0, ∂ t (hu) + hu(∇ • u) + ∇p = -gh∇z, (5)
and on the other hand the transport step

∂ t h + (u • ∇)h = 0, ∂ t (hu) + (u • ∇)(hu) = 0. (6)
With these notations, the Lagrange-projection algorithm is defined as follows: for a given discrete state (h, hu) n j , j ∈ Z, defining (h, hu) n+1 j is a two-step process defined as follows 1. Update (h, hu) n j to (h, hu) n+1- j by approximating the solution of system (5),

2. Update (h, hu) n+1- j to (h, hu) n+1 j
by approximating the solution of system (6).

Relaxation approximation of the acoustic system. Before entering the details of these two steps in the following section, let us note that if we denote τ = 1/h, by simple manipulations system (5) can be recast into:

∂ t τ -τ (x, t)∇ • u = 0, ∂ t u + τ (x, t)∇p = -τ (x, t) g τ ∇z.
Following [START_REF] Coquel | Entropy-satisfying relaxation method with large time-steps for Euler IBVPs[END_REF], we will choose to approximate the solution of system (5) thanks to a Suliciurelaxation process. More precisely we will solve

       ∂ t τ -τ (x, t)∇ • u = 0, ∂ t u + τ (x, t)∇Π = -τ (x, t) g τ ∇z, ∂ t Π + τ (x, t)a 2 ∇ • u = λ(p EOS (τ ) -Π), (7)
with p = p EOS (τ ) = g/(2τ 2 ), in the regime λ → +∞. The parameter a is a constant that is chosen in agreement with the subcharacteristic stability conditions that will be given later. Over the time interval [t n , t n + ∆t), we can account for the limit λ → +∞ by setting Π(x,

t n ) = p EOS (τ (x, t n )),
and then solving the relaxed system with λ = 0. We add another approximation by supposing that over

[t n , t n + ∆t) it is reasonable to replace τ (x, t)∂ xr by τ (x, t n )∂ xr , r = 1, 2.
Finally, we will define our approximation of the acoustic system (5) by solving

       ∂ t τ -τ (x, t n )∇ • u = 0, ∂ t u + τ (x, t n )∇Π = -τ (x, t n ) g τ ∇z, ∂ t Π + τ (x, t n )a 2 ∇ • u = 0, (8)
over [t n , t n + ∆t), with Π(x, t n ) = p EOS (τ (x, t n )).
Note that the system (8) is rotational invariant. This will allow us in the following to define every flux for the two dimensional problem in the reference frame associated to each face. In this last referential, the problem will be reduced to a quasi-one dimensional problem that we study in the beginning of next section. One can also notice that the eigenstructure of (8) in the phase space (h, hu T , Π, z) ∈ R 5 , h > 0, z > 0 is very simple since it has three eigenvalues {-a, 0, a} all associated to linearly degenerated characteristic fields.

Finite volume approximation

In this paragraph, we present in details the first-order finite volume scheme associated with the acoustic / transport decomposition of section 3.

A well-balanced Lagrange-projection finite volume scheme in 1D

We start by considering one-dimensional problems and briefly recall the method proposed in [START_REF] Chalons | A large time-step and well-balanced Lagrange-projection type scheme for the shallow water equations[END_REF]. In this case the Saint-Venant equations read

         ∂ t h + ∂ x (hu 1 ) = 0, ∂ t (hu 1 ) + ∂ x hu 2 1 + g h 2 2 = -gh∂ x z, ∂ t (hu 2 ) + ∂ x (hu 1 u 2 ) = 0.
The system associated with the acoustic step reads

             ∂ t τ -τ (x, t n )∂ x u 1 = 0, ∂ t u 1 + τ (x, t n )∂ x Π = -τ (x, t n ) g τ ∂ x z, ∂ t u 2 = 0, ∂ t Π + τ (x, t n )a 2 ∂ x u 1 = 0,
and the system that accounts for transport boils down to

∂ t ϕ + u 1 ∂ x ϕ = 0, ϕ ∈ {h, u 1 , u 2 }.
We suppose given a strictly increasing sequence x j+1/2 ∈ R, for j ∈ Z and we consider the set of cells Ω j = [x j-1/2 , x j+1/2 ). The local space step is defined by ∆x j = x j+1/2 -x j-1/2 . We note ∆t > 0 the time step and we set t n = n∆t for n ∈ N.

The following discretization strategy was presented in [START_REF] Chalons | A large time-step and well-balanced Lagrange-projection type scheme for the shallow water equations[END_REF]: the acoustic step (8) is approximated by

                       τ n+1- j = τ n j -τ n j ∆t ∆x j u j+1/2 -u j-1/2 (u 1 ) n+1- j = (u 1 ) n j -τ n j ∆t ∆x j Π L, j+1/2 -Π R, j-1/2 (u 2 ) n+1- j = (u 2 ) n j Π n+1- j = Π n j -τ n j ∆t ∆x j a 2 u j+1/2 -u j-1/2 (9a) (9b) (9c) (9d)
where for all j, Π n j = g

h n j 2
2 and the numerical fluxes u j+1/2 , Π L, j+1/2 and Π R, j-1/2 are defined by

u j+1/2 = u ∆ (U j , U n j , U j+1 , U n j+1 ), Π R, j+1/2 = Π R ∆ (U j , U n j , U j+1 , U n j+1 ), Π L, j+1/2 = Π L ∆ (U j , U n j , U j+1 , U n j+1 ),
where U is the state

    h hu Π z     and with {gh∆z} ∆ (U n L , U n R ) = g h n L + h n R 2 (z R -z L ) Π ∆ (U L , U R ) = Π L + Π R 2 -a (u 1 ) R -(u 1 ) L 2 u ∆ (U L , U n L , U R , U n R ) = (u 1 ) L + (u 1 ) R 2 - Π R -Π L 2a - 1 2a {gh∆z} ∆ (U n L , U n R ) Π L ∆ (U L , U n L , U R , U n R ) = Π ∆ (U L , U R ) + 1 2 {gh∆z} ∆ (U n L , U n R ) Π R ∆ (U L , U n L , U R , U n R ) = Π ∆ (U L , U R ) - 1 2 {gh∆z} ∆ (U n L , U n R )
If one chooses = n (resp. = n + 1 -) the system (9) provides a time-explicit (resp. time-implicit) discretization of the acoustic system (8). The approximation of the transport step is performed thanks to a standard upwind scheme for ϕ ∈ {h, hu 1 , hu 2 }

ϕ n+1 j = ϕ n j - ∆t ∆x j u j+1/2 ϕ n+1- j+1/2 -u j-1/2 ϕ n+1- j-1/2 - ∆t ∆x j ϕ n+1- j u j+1/2 -u j-1/2 , (10)
where

ϕ n+1- j+1/2 = ϕ n+1- j , if u j+1/2 ≥ 0, ϕ n+1- j+1 , if u j+1/2 < 0.
Note that {gh∆z} ∆ that accounts for the gravity source term is always evaluated at time t n , even for the time implicit scheme.

In the above formulas, the parameter a is an approximation of the Lagrangian sound speed hc = h √ gh and must satisfy the sub-characteristic condition a > hc which ensures that the relaxed system (7) is a dissipative approximation of the acoustic step (5) (see [Bou04, CC08, CC05, Des10] and the references therein). In order to limit the numerical diffusion we take a local approximation of the Lagrangian sound speed at every interface, given by a j+1/2 = κ max(h j gh j , h j+1 gh j+1 ), where κ > 1.

For detailed properties of the numerical scheme ( 9)-( 10) we refer the reader to [START_REF] Chalons | A large time-step and well-balanced Lagrange-projection type scheme for the shallow water equations[END_REF], nevertheless let us recall that: the overall discretization is conservative in the usual sense of finite volumes methods with respect to (h, hu 1 , hu 2 ). Moreover, the scheme is also well-balanced for lake at rest conditions: if u n j = 0 and h n j + z n j = h n j+1 + z n j+1 for all j ∈ Z, then h n+1 j = h n j and u n+1 j = u n j , j ∈ Z. At last, the time-implicit scheme is stable under a condition which does not depend either on the acoustic system or the sound speed c, but which only depends on the transport step and its material velocity u which is of particular interest in the low-Froude regime.

Truncation error in the low-Froude regime

In this paragraph, we consider the dimensionless shallow-water equation and we motivate a correction of the above scheme in order to make it efficient in low-Froude regimes. The correction is similar to the one in [START_REF] Chalons | An all-regime Lagrangeprojection like scheme for the gas dynamics equations on unstructured meshes[END_REF] for low-Mach regimes and we focus on the explicit case = n.

In the following we will say that the flow is in the low Froude regime if Fr 1 and ∂ x p+h∂ x z = O(Fr 2 ). Regarding the dimensionless equations (2), we can observe that, in this regime, the variations of the discharge hu remain of order 1 as expected.

We can express the fluxes given in the previous section, using the dimensionless quantities, which leads to

u n j+1/2 = 1 2 (u n j + u n j+1 ) - 1 2aFr Π n j+1 -Π n j + h n j + h n j+1 2 (z j+1 -z j ) , Π L,n j+1/2 = Π n j Fr 2 + 1 2Fr 2 Π n j+1 -Π n j + h n j + h n j+1 2 (z j+1 -z j ) - a 2Fr (u n j+1 -u n j ), Π R,n j+1/2 = Π n j+1 Fr 2 - 1 2Fr 2 Π n j+1 -Π n j + h n j + h n j+1 2 (z j+1 -z j ) - a 2Fr (u n j+1 -u n j ),
if one focuses on the time-explicit scheme for the sake of simplicity.

If we compute the truncation errors in the fluxes above, using the fact that

Π n j+1 -Π n j + h n j + h n j+1 2 (z j+1 -z j ) = O(Fr 2 ∆x),
we obtain:

u n j+1/2 = 1 2 (u n j + u n j+1 ) + O(Fr∆x), Π L,n j+1/2 = Π n j Fr 2 + 1 2Fr 2 Π n j+1 -Π n j + h n j + h n j+1 2 (z j+1 -z j ) + O( ∆x Fr ), Π R,n j+1/2 = Π n j+1 Fr 2 - 1 2Fr 2 Π n j+1 -Π n j + h n j + h n j+1 2 (z j+1 -z j ) + O( ∆x Fr ).
At this stage, it is clear that the consistence errors are not uniform with respect to the Froude number in the pressure fluxes. In order to avoid large errors in the numerical diffusion terms when the Froude tends to zero, we propose to correct the flux formula of Π ∆ (U L , U R ) by:

Π θ ∆ (U L , U R ) = Π L + Π R 2 -θa (u 1 ) R -(u 1 ) L 2
which amounts to reduce the numerical diffusion on the pressure gradient in the low Froude regime. Indeed, we now get

Π L,n,θ j+1/2 = Π n j Fr 2 + 1 2Fr 2 Π n j+1 -Π n j + h n j + h n j+1 2 (z j+1 -z j ) + O( θ j+1/2 ∆x Fr ), Π R,n,θ j+1/2 = Π n j+1 Fr 2 - 1 2Fr 2 Π n j+1 -Π n j + h n j + h n j+1 2 (z j+1 -z j ) + O( θ j+1/2 ∆x Fr ),
and as long as we take θ j+1/2 = O(Fr), we recover the uniform consistency of the global scheme with respect to the Froude number. In practice, we will set θ j+1/2 = min

|u n j+1/2 | max(cj ,cj+1) , 1 .

The Lagrange-projection scheme on 2D unstructured meshes

We now extend the Lagrange-projection scheme in two dimensions. Let n ∈ R 2 be a unit vector and U T = (h, hu T ), we define

R n = n 1 n 2 -n 2 n 1 , T n U =     h h(R n u) Π z     .
Following standard lines, we take advantage of the rotational invariance of the acoustic system (8) to define the two-dimensional fluxes numerical fluxes (see for example [START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF][START_REF] Bouchut | Nonlinear stability of finite Volume Methods for hyperbolic conservation laws: And Well-Balanced schemes for sources[END_REF]). This leads to

                   τ n+1- j = τ n j + τ n j ∆t k∈N (j) σ jk u jk , u n+1- j = u n j -τ n j ∆t k∈N (j) σ jk Π ,θ jk n jk , Π n+1- j = Π n j -τ n j ∆t k∈N (j) σ jk (a jk ) 2 u jk , (11a) (11b) (11c)
where

u jk = u ∆ (T n jk U j , T n jk U n j , T n jk U k , T n jk U n k ), Π ,θ jk = Π L,θ ∆ (T n jk U j , T n jk U n j , T n jk U k , T n jk U n k ),
that is to say

u jk = 1 2 n T jk (u j + u k ) - 1 2a jk (Π k -Π j ) - 1 2a jk {gh∆z} n jk , Π ,θ jk = 1 2 (Π j + Π k ) - a jk θ jk 2 n T jk (u k -u j ) + 1 2 {gh∆z} n jk , with
a jk ≥ max[(hc) n j , (hc) n k ], {gh∆z} n jk = g h n j + h n k 2 (z k -z j ).
The source term is accounted for by the terms Π ,θ jk since the fluxes Π ,θ jk n jk in Equation (11b) are not symmetric, indeed Π ,θ jk n jk = -Π ,θ kj n kj (even if a jk = a kj and θ jk = θ kj , which is the case in practice).

As far as the transport step is concerned and in order to discretize the system (6), we use an explicit scheme between times t n+1-and t n+1-+ ∆t, where the fluxes are chosen upwind with respect to the sign of u jk . If ϕ ∈ {h, hu}, the scheme for the transport step reads

ϕ n+1 j = ϕ n+1- j -∆t k∈N (j) σ jk ϕ n+1- jk u jk + ∆tϕ n+1- j k∈N (j) σ jk u jk , (12)
where

ϕ n+1- jk = ϕ n+1- j , if u jk ≥ 0, ϕ n+1- k , if u jk < 0.
Note that one can rewrite the transport step (12) as follows

ϕ n+1 j = L j ϕ n+1- j -∆t k∈N (j) σ jk u jk ϕ n+1- jk ,
where L j = 1 + ∆t k∈N (j) σ jk u jk . Therefore, replacing the quantities ϕ n+1- j using system (11) gives the following update formulas which take into account the acoustic and transport steps together:

           h n+1 j = h n j -∆t k∈N (j) σ jk h n+1- jk u jk , (hu) n+1 j = (hu) n j -∆t k∈N (j) σ jk (hu) n+1- jk u jk + Π ,θ jk n jk , (13a) (13b)
where the quantities u jk and Π ,θ jk are computed with u and Π , and the quantities h n+1- jk and (hu) n+1- jk with τ n+1-and u n+1-from system (11).

Stability and well-balanced properties

Let us first notice that the scheme is conservative with respect to the water height h, and with respect to hu if the topography is flat (z = cste). In particular, it degenerates towards the scheme proposed by [START_REF] Chalons | An all-regime Lagrangeprojection like scheme for the gas dynamics equations on unstructured meshes[END_REF] adapted to the framework of barotropic Euler system when the bottom is flat. Next, recall that from section 4.2, if θ is chosen to be like O(Fr) when Fr goes to 0, the truncation erro of the numerical scheme is uniform with respect to Fr. At last, assuming that the time step ∆t is such that the CFL conditions associated to the acoustic step

∆t max 1≤j≤N τ n j max k∈N (j) σ jk a jk ≤ 1 2 ,
and to the transport step

∆t max 1≤j≤N   k∈N (j),u n jk <0 σ jk |u n jk |   ≤ 1,
hold true, the water height h n j is positive for all j and n > 0 provided that h 0 j is positive for all j, for the time-explicit scheme corresponding to = n. Indeed, notice that L n j turns out to be positive while the transport step correspond to a convex combination of states at time t n+1-.

As far as the mixed implicit-explicit scheme corresponding to = n + 1-is concerned, the same properties hold true under the transport CFL condition

∆t max 1≤j≤N    k∈N (j),u n+1- jk <0 σ jk |u n+1- jk |    ≤ 1.
Notice that the acoustic step is implicit and therefore is free of CFL condition. Now we show that the schemes are well-balanced. We begin with the explicit scheme = n.

Proposition 1. The full explicit scheme ( = n) is well-balanced on 2D unstructured mesh in the sense that (h 0 j + z j = H = cste and u 0 j = 0) =⇒ (h n j + z j = cste and u n j = 0) ∀j. Proof. Assume that h 0 j + z j = H = cste and u 0 j = 0 ∀j. We have

{gh∆z} 0 jk = g h 0 j + h 0 k 2 (H -h 0 k ) -(H -h 0 j ) = Π 0 j -Π 0 k , u 0 jk = - 1 2a jk Π 0 k -Π 0 j - 1 2a jk {gh∆z} 0 jk = 0, Π 0,θ jk = 1 2 Π 0 j + Π 0 k + 1 2 {gh∆z} 0 jk = Π 0 j .
Injecting those values in the acoustic step gives:

           u 1- j = u 0 j -τ 0 j ∆t k∈N (j) σ jk Π 0 j n jk = -τ 0 j Π 0 j ∆t k∈N (j) σ jk n jk = 0, Π 1- j = Π 0 j , τ 1- j = τ 0 j .
Next, since u 0 jk = 0, ∀j, k, the transport step is trivial and the variables h and hu are unchanged at time t 1 .

Proposition 2. The implicit-explict scheme ( = n + 1-) is well-balanced on 2D unstructured mesh in the sense that (h 0 j + z j = H = cste and u 0 j = 0) =⇒ (h n j + z j = cste and u n j = 0) ∀j.

Proof. With the same calculus as in the explicit case, we can verify that the vector (u 1- j , Π 1- j ) = (0, Π 0 j ) is the only solution of the coupled system over (u, Π). Indeed we have in this case

{gh∆z} 0 jk = g h 0 j + h 0 k 2 (H -h 0 k ) -(H -h 0 j ) = Π 0 j -Π 0 k , u 1- jk = - 1 2a jk Π 1- k -Π 1- j - 1 2a jk {gh∆z} 0 jk = 0, Π 1 -,θ jk = 1 2 Π 1- j + Π 1- k + 1 2 {gh∆z} 0 jk = Π 0 j .
and therefore

     u 1- j = u 0 j -τ 0 j ∆t k∈N (j) σ jk Π 1 -,θ jk n jk = -τ 0 j Π 0 j ∆t k∈N (j) σ jk n jk = 0, Π 1- j = Π 0 j .
Then we easily get τ 1- j = τ 0 j so that the Lagrangian step is well-balanced. As before, we can conclude that the transport step is also well-balanced since u 1- jk = 0, ∀j, k.

Numerical experiments

We present several test cases that aim at testing our scheme against classical flow configurations on unstructured meshes and also in the low Froude regime.

In the following, the EXEX scheme will refer to the full explicit scheme ( = n) with the time step ∆t defined by

∆t EXEX = K CFL 2 max j k∈N (j) |Γ jk | |Ωj | max k∈N (j) max(v Acou jk , v Trans jk )
,

where K CFL = 0.9, v Acou jk = τ j a jk , a jk = 1.01 max(h j c j , h k c k ), c j = gh j and v Trans jk = |u n jk |. The IMEX scheme will refer to the implicit-explicit scheme ( = n + 1-) with the time step ∆t defined by Except if otherwise stated, we will always use the corrected numerical fluxes with

∆t IMEX = K CFL 2 max j k∈N (j) |Γ jk | |Ωj | max k∈N (j) (v Trans jk ) ρ -0.5 ∞ u ∞ EXEX 2.
θ jk = min |u n jk | max (c j , c k )
, 1 , so that θ approximates a local Froude number on every edge. On the other hand, for the sake of comparison we take θ = 1 on every edge for the fluxes without correction.

Test of the well-balanced property

In order to test the well-balanced property of the scheme, we first consider the following lake at rest initial condition:

h(x, y, 0) = H -z(x, y), u(x, y, 0) = 0,
where H = 0.5 is constant and the topography z is a smooth bump defined by

z(x, y) = 0.3 ×                   
0.5 exp(2 -0.1 x-0.325 ), if 0.325 < x ≤ 0.375, 1 -0.5 exp(2 -0.1 0.425-x ), if 0.375 < x < 0.425, 1, if 0.425 ≤ x ≤ 0.575, 1 -0.5 exp(2 -0.1

x-0.575 ), if 0.575 < x < 0.625, 0.5 exp(2 -0.1 0.675-x ), if 0.625 ≤ x < 0.675, 0 otherwise.

(

) 14
The physical domain [0, 1] × [0, 1] is discretized over a 20 000-cell triangular mesh. We impose Neumann boundary conditions and we observe the solution at final time T f = 0.1.

For both EXEX and IMEX schemes, the errors between the numerical and the exact solution, which is also the initial stationary condition, are machine epsilon as we can observe in Table 1.

Planar dam break test problem

We are interested in the behaviour of our schemes with regard to the propagation of a rarefaction wave and a shock wave. We use the same triangular mesh, boundary conditions and T f value as in section 5.1. The topography is also kept identical to the one given in ( 14), the velocity initialized to zero and the initial total water height H = h + z is defined as follows:

H(x, y, 0) = 0.5 if x ≤ 0.5, 1 otherwise.

In Figure 1 we present the results for both the EXEX and IMEX schemes. We have performed a cut of the solution along the y = 0.5 axis and compared it with the one computed by a genuine 1D code with a 200-cell uniform grid. We can observe that for both EXEX and IMEX the results of the 2D simulations are in agreement with the 1D results although they were computed with an unstructured mesh. It is worth noting that the 2D simulation manages to fairly preserve the planar structure of the approximate solutions.

Traveling vortex with flat bottom

In order to challenge our schemes with low Froude regimes, we consider a traveling vortex as in [START_REF] Bispen | Imex large time step finite volume methods for low Froude number shallow water flows[END_REF]. The exact solution of this test is detailed in [START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF]. For this test case we consider a flat bottom and we use a regular cartesian mesh of 160 × 160 cells that discretizes the physical domain [0, 1] × [0, 1]. The boundary conditions imposed are periodic along the x-direction and absorbing boundaries along the y-direction. The initial conditions are given by: Due to the periodic boundary conditions, the exact solution is periodic with period T = 5 3 and given at any time t > 0 by:

h(x, y, 0) = 110 + Γ 2 gω 2 (k(ωr c ) -k(π)) if ωr c ≤ π, 0 otherwise, u(x, y, 0) = 0.6 + Γ (1 + cos(ωr c )) (0.5 -y) if ωr c ≤ π, 0 otherwise, v(x, y, 0) = 0 + Γ (1 + cos(ωr c )) (x -0.5) if ωr c ≤ π, 0
h(x, y, t) = h(x -t/T, y, 0), u(x, y, t) = u(x -t/T, y, 0), v(x, y, t) = v(x -t/T, y, 0).
We present the results of both the EXEX and IMEX schemes, with (θ = O(Fr)) and without correction (θ = 1) using = 0.05. The mapping of the velocity magnitude is displayed in Figure 2 and we can observe that the accuracy of the solution is really improved by the low-Froude correction. Furthermore, the accuracy of the solution between the EXEX and the IMEX scheme with low-Froude correction is comparable whereas it took about 100 times less time steps and 10 times less CPU time computation to reach the final time with the IMEX than with the EXEX scheme as we can see in Table 2.

Traveling vortex with non-flat bottom

We extend the physical domain of the traveling vortex test above to the rectangle [0, 2] × [0, 1]. The boundary conditions and initial conditions for h and u are the same as in section 5.3. However we consider here a topography defined by z(x, y) = 10 exp -5(x -1) 2 -50(y -0.5) 2 following the idea of [START_REF] Bispen | Imex large time step finite volume methods for low Froude number shallow water flows[END_REF].

We do not have exact analytical solution because of the non-flat bottom but we still can compare in figure 3 the results between EXEX and IMEX schemes, with or without low Froude correction. Here again, the vortex structure of the flow is completely destroyed by numerical diffusion without low-Froude corrections θ = O(Fr), with both schemes. The mapping of the Froude number is not given here, but is similar to the one of the velocity magnitude, with a range of values from 1.6•10 -3 to 1.1 • 10 -2 . Finally, we can remark that the EXEX scheme took about 20 times more iterations and 15 times more CPU times than the IMEX scheme, both with low-Froude correction, as we can see in Table 3.

Conclusion

We have proposed a large time step and well-balanced scheme for the shallow-water equations in two dimensions for unstructured meshes. We studied the truncation of the scheme with respect to the Froude number Fr and gave a correction in accordance to the source term. By studying the one-dimensional case, we obtained proposed of modification of the scheme that allows to obtain a uniform truncation error with respect to the Froude number for one-dimensional flows.

Moreover, we showed that the semi-implicit scheme yields good numerical results for flows from low to high Froude values since its CFL condition is based on (slow) material waves only.

Further developments shall include extensions to high-order methods in multiple-dimensions, following for example what has already been achieved with Finite Volume or discontinuous Galerkin methods in 1D. We also intend to adapt the method presented in this work to other compressible flows models involving non-conservative terms.

Figure 1 :

 1 Figure 1: Dam break test case at T f = 0.1: mapping of the total water height with the IMEX scheme (top) and profile of H and z along the y = 0.5 axis obtained with both the EXEX and IMEX with 2D and 1D simulations (bottom).

  otherwise, where r c = x -(0.5, 0.5) , Γ = 15.0, ω = 4π, Traveling vortex test case with flat bottom. Numbers of iterations and CPU times with low-Froude correction.

Figure 2 :

 2 Figure 2: Traveling vortex test case with flat bottom. Mapping of the velocity magnitude at T f = 0.1 obtained with the EXEX scheme (top) and the IMEX scheme (bottom). We used the values θ = 1 (left) and θ = O(Fr) (center). The right column displays the exact solution.

  Traveling vortex test case with non-flat bottom. Numbers of iterations and CPU times with low-Froude correction.

Figure 3 :

 3 Figure 3: Traveling vortex test case with non-flat bottom. Mapping of the velocity magnitude at instant T f = 0.1 obtained with the EXEX scheme (top) and IMEX scheme (bot) with θ = 1 (left) and θ = O(Fr) (right).

Table 1 :

 1 Well-balanced property. Errors for EXEX and IMEX schemes.

				6 10 -16	1.3 10 -13
			IMEX	2.6 10 -16	3.9 10 -8
	where v Trans jk	= |u n+1-jk	|. Thus, the time step of the IMEX scheme is not constrained by the acoustic
	waves.


archive ouverte pluridisciplinaire HAL, est destinée au dépôt

Polyelectrolytes Multilayer Films on titanium and titanium alloys for biomedical applications

INTRODUCTION:

The aim of this research was to optimize titanium (Ti) and titanium alloy (NiTi) surfaces coated with Polyelectrolytes Multilayer Films (PMF) [1] . We studied different biomedical parameters related to this surface treatment in order to develop specific biomedical applications in the dental field (dental implants, endodontic instruments, orthodontic arches). Firstly, we determined if PMF had a detectable physisorbtion on Ti and NiTi. Secondly, we studied biocompatibility of poly-ethyleneimin (PEI), the precursor basedlayer of PMF. Finally, sterilization tests were realized to analyze their potential impact on the physico-chemical structure of PMF. METHODS: Two types of films were caracterized : polystyrenesulfonate/polyallylamine hydrochloride films PEI-(PSS/PAH) and hyaluronic acid/poly-Llysine films PEI-(HA/PLL). Physico-chemical characterization was carried out by tensiometry, atomic force microscopy (AFM), and confocal microscopy. A biological study using human fibroblasts was carried out. Cell response was observed after 0, 2, 4 and 7 days in vitro using morphologic criteria (scanning electron microscopy), adhesion (fluorescence microscopy image analysis), and proliferation (Methyl Tetrazolium Test). Three methods of PMF sterilization were compared: ultraviolets (UV), ethylen oxyde (ETO), and autoclave. PSS/PAH films were characterized by the same physico-chemical tests and biological studies after the sterilization process. RESULTS: Results showed that PSS/PAH films were more biocompatible than HA/PLL films on both metallic biomaterials. The precursor based-layer study demonstrated that PEI was not biocompatible [3] . PMF sterilized by autoclave showed similar biocompatibility. Such results could not be found with UV and ETO. Nevertheless, some investigations have to be realized to prove that the structure of PMF was not perturbed after these processes. DISCUSSION & CONCLUSION: Based on our results, we decided to use PAH as the precursor based-layer (positively charged) and to stop using PEI on biomaterials made of Ti or NiTi. We are now planning to integrate bioactive molecules between layers of the PMF. To this aim other parameters need to be studied to characterize in vivo potential biomedical application, such as films aging in the oral cavity, as well as in salivary or fluored environment.

REFERENCES: [1] Decher G. Fuzzy nanoassemblies: toward layered polymeric multicomposites. Science 1997,277:1232-7. [2] Brunot, C, Grosgogeat B, Picart C, Lagneau C, Jaffrezic-Renault N, Ponsonnet L. Response of fibroblast activity and polyelectrolyte multilayer films coating titanium. Dent Mater 2008,24(8):1025-35.

ACKNOWLEDGEMENTS:

Timet ® Savoie society (Ugine, France) and AMF ® society (Lury sur Arnon, France) for providing samples.

Polyelectrolytes

  Multilayer Films on titanium and titanium alloys for biomedical applications C Brunot, L Mora-Ponsonnet, D Décoret, C Picart, B Beaugiraud, Brigitte Grosgogeat To cite this version: C Brunot, L Mora-Ponsonnet, D Décoret, C Picart, B Beaugiraud, et al.. Polyelectrolytes Multilayer Films on titanium and titanium alloys for biomedical applications. eCells and Materials Journal, 2014, eCM-2014, 27 (suppl.1), pp.17. hal-02004864

  Brunot 1,2 , L. Mora-Ponsonnet 3 , D. Décoret 1 , C. Picart 4 , B. Beaugiraud 5 , et B. Grosgogeat 1,6 Laboratoire des Matériaux et du Génie Physique, UMR CNRS 5628, MINATEC, 3 parvis Louis Néel, 38016 Grenoble. 5 Laboratoire de Tribologie et Dynamique des Systèmes, Ecole Centrale de Lyon, 36 avenue Guy de Collongue, 69134 Ecully Cedex. 6 Hospices Civils de Lyon, Service de Consultations et de Traitements Dentaires, 6-8 Place Depéret, 69365 Lyon Cedex 07.

	1 Laboratoire des Multimatériaux et Interfaces, UMR CNRS 5615, Université de Lyon, Equipe Biomatériaux
	et Interfaces Biologiques, UFR d'Odontologie, rue Guillaume Paradin, 69372 Lyon Cedex 08. 2 Faculté
	d'Odontologie, Université de Reims Champagne-Ardenne, 2 rue du Général Koenig, 51100 Reims.
	3 Laboratoire de Bio-Ingénierie de Polymères Cardiovasculaires, INSERM U 698, Institut Galilée, Université
	Paris13, 99 avenue JB Clément, 93430 Villetaneuse. 4

C.

L'


Alexandre Hollocou

Thomas Bonald

Marc Lelarge

Modularity-based Sparse Soft Graph Clustering

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Modularity is a quality function defined on node partitions that has been widely used to tackle the problem of node clustering in graphs [START_REF] Fortunato | Community structure in graphs[END_REF][START_REF] Newman | Finding and evaluating community structure in networks[END_REF]. Numerous algorithms have been proposed to solve the problem of maximizing this objective function [START_REF] Duch | Community detection in complex networks using extremal optimization[END_REF][START_REF] Guimera | Modularity from fluctuations in random graphs and complex networks[END_REF][START_REF] Newman | Fast algorithm for detecting community structure in networks[END_REF]. The Louvain algorithm [START_REF] Blondel | Fast unfolding of communities in large networks[END_REF], which is built in a number of graph analysis software packages, is arguably the most popular approach to solve this problem. A solution to the modularity optimization problem is a partition of the graph nodes. In particular, a node cannot belong to more than one cluster. Yet, in many real-life applications, it makes more sense to allow some elements to belong to multiple clusters. For instance, if we consider the problem of clustering the users of a social network such as Facebook or LinkedIn, we see that a person can belong to several social circles, such as her family, her colleagues and her friends, and we might be interested in recovering all these circles with a clustering algorithm, even if they do not form a partition. In order to overcome this limit of the classical definition of modularity, we introduce a relaxation of the problem to perform a so-called soft clustering that allows partial membership like fuzzy clustering methods [START_REF] Bezdek | Objective function clustering[END_REF][START_REF] Dunn | A fuzzy relative of the isodata process and its use in detecting compact well-separated clusters[END_REF]. A solution to this new problem gives for each node a degree of membership to each cluster, instead of returning a simple node partition. More precisely, we obtain a probability p ik ∈ [0, 1] for each node i to belong to a given cluster k, where the classical modularity-based techniques implicitly consider that this probability can be either 0 or 1.

The main contribution of this paper is to introduce an efficient algorithm to find an approximation of this soft version of the modularity maximization problem. This algorithm has four main advantages. First, the number of clusters does not need to be specified. Second, the algorithm is local, in the sense that each update of the membership information of a given node depends only on its direct neighbors in the graph. Third, the solution found by the algorithm is sparse, i.e. most of the membership probabilities p ik returned by our algorithm are equal to 0, which guarantees an efficient storage of the solution. Finally, an update performed by our algorithm reduces to an update performed by the Louvain algorithm as soon as the unique parameter t of our algorithm is larger than a threshold w/δ that depends on the graph.

The remainder of the paper is organized as follows. In section 2, we present the related work on the topic of modularity optimization. We introduce the relaxation of the modularity maximization problem in section 3. In section 4, we present our optimization method and obtain theoretical guarantees about its convergence. In section 5, we take benefits of the specificities of our optimization technique to propose an algorithm that is both local and memory efficient. In section 6, we study the ties of our approach to the Louvain algorithm.

Finally, in section 7, we present experimental results on synthetic and real-world data, and section 8 concludes the paper. Complete proofs of the main results of the paper are presented in the supplementary material.

In the rest of the paper, we consider that we are given a weighted and undirected graph G = (V, E, W ). We use V to denote the set of nodes, E the set of edges, and W = (W ij ) i,j∈V the adjacency matrix of the graph. If (i, j) ∈ E, W ij > 0 is the weight of edge (i, j), and if (i, j) / ∈ E, W ij = 0. We use w i to denote the weighted degree of node i, w i = j∈V W ij , and w the total weight of the graph w = i∈V w i . We use Nei(i) to denote the set of the neighbors of node i ∈ V , and j ∼ i as a notation for j ∈ Nei(i). Finally, we use n = |V | to denote the number of nodes, and m = |E| the number of edges of the graph.

Related work

Modularity optimization

The modularity function (1) has first been introduced by Newman and Girvan [START_REF] Newman | Finding and evaluating community structure in networks[END_REF] to measure the quality of a partition of graph nodes. For a given node partition P , modularity is defined as

Q(P ) = 1 w C∈P i,j∈C W ij - w i w j w . (1)
Modularity Q can be interpreted as the difference between the probability to pick an edge in G between two nodes of the same cluster C, i.e. an intra-cluster edge, and the probability to pick such an intra-cluster edge in a random graph generated with the so-called null model, where the probability that an edge (i, j) exists is w i w j /w.

Clustering approaches based on modularity try to solve the modularity maximization problem max

P ∈P Q(P ), (2)
where P is the set of the partitions of V .

We refer to the Newman and Girvan modularity as the hard modularity, to distinguish it from the soft modularity that we introduce in Section 3. The hard modularity optimization problem (2) cannot be solved exactly in real applications as it has been proven to be NP-hard [START_REF] Brandes | On finding graph clusterings with maximum modularity[END_REF]. However, many methods have been proposed to find good approximations of the modularity optimum in reasonable time. These techniques includes greedy algorithms [START_REF] Newman | Fast algorithm for detecting community structure in networks[END_REF], spectral approaches [START_REF] Newman | Finding community structure in networks using the eigenvectors of matrices[END_REF], and simulated annealing methods [START_REF] Guimera | Modularity from fluctuations in random graphs and complex networks[END_REF]. A thorough review of existing algorithms can be found in [START_REF] Fortunato | Community structure in graphs[END_REF]. One of the most popular algorithm for finding a good approximation of the modularity maximum is the Louvain algorithm [START_REF] Blondel | Fast unfolding of communities in large networks[END_REF], which is widely used in benchmarks for its ability to scale up to very large graphs (e.g. the Friendster social graph with 1.8 billions edges [START_REF] Prat-Pérez | High quality, scalable and parallel community detection for large real graphs[END_REF]). In this paper, we introduce an algorithm for optimizing the relaxation of the modularity maximization problem (4) that takes a parameter t, and we prove that each update performed by our algorithm reduces to the update performed by the Louvain algorithm if t is large enough. This constitutes a strong guarantee since the Louvain algorithm has proven its worth on multiple real-life applications.

Modularity relaxation

Relaxations of the hard modularity optimization problem have been tackled by few articles [START_REF] Chang | A unified framework for sampling, clustering and embedding data points in semi-metric spaces[END_REF][START_REF] Griechisch | Community detection by using the extended modularity[END_REF][START_REF] Havens | A soft modularity function for detecting fuzzy communities in social networks[END_REF][START_REF] Kawase | Additive approximation algorithms for modularity maximization[END_REF][START_REF] Lázár | Modularity measure of networks with overlapping communities[END_REF][START_REF] Lehmann | Deterministic modularity optimization[END_REF][START_REF] Nicosia | Extending the definition of modularity to directed graphs with overlapping communities[END_REF]. They all introduce membership matrices p ∈ R n×K , with K ≥ 1, where p ik ≥ 0 represents the degree of membership of node i ∈ V to the k th cluster, k ∈ 1, K . K is a given positive integer, that represents the maximum number of clusters. These relaxations all have similar forms and are very close to the one we introduce in Section 3. These works differ in the optimization techniques they use to solve the problem. In [START_REF] Chang | A unified framework for sampling, clustering and embedding data points in semi-metric spaces[END_REF], Chang et al. use a softmax parameterization for p, defining for each i ∈ V and k p ik = exp(θ ik )/( l exp(θ il )) where θ ∈ R n×K is the parameter to optimize. In [START_REF] Nicosia | Extending the definition of modularity to directed graphs with overlapping communities[END_REF], Nicosia et al. use a genetic algorithm to optimize a more general relaxation of the modularity problem. In [START_REF] Griechisch | Community detection by using the extended modularity[END_REF], Griechisch et al.

do not directly study the optimization of the relaxation problem, but they rely on an external quadratic solver. Finally, in [START_REF] Havens | A soft modularity function for detecting fuzzy communities in social networks[END_REF], Havens et al. use a spectral approach that does not directly solve the relaxation of the modularity problem, but where modularity is used as a selection criterion. The main limitation of these methods lies in the maximum number of clusters K that must be specified. We could get around this issue by taking large values for K, but all the approaches cited above do not scale well to large K. Indeed, the solutions p ∈ R n×K found by these methods are dense matrices. In other words, the number of parameters to store in memory and to optimize is in O(nK), which quickly becomes prohibitive for large values of K. For instance, in the approach of [START_REF] Chang | A unified framework for sampling, clustering and embedding data points in semi-metric spaces[END_REF], the matrix p is the densest possible matrix, i.e. all its coefficients are positive.

In the approach of [START_REF] Nicosia | Extending the definition of modularity to directed graphs with overlapping communities[END_REF], the genetic algorithm starts with dense random matrices of R n×K , and its hybridation and mutation mechanisms do not lead to sparser matrices, so that all the coefficients of the solution p found by this algorithm are positive with probability 1.

In this paper, we introduce an efficient algorithm to solve a relaxation of the modularity optimization problem with a membership matrix p ∈ R n×n . Thus, the maximum number of clusters does not need to be spec-ified. Besides, unlike the approaches presented above, the updates performed by our algorithm preserve the sparsity of the solution p and are local, in the sense that the membership vector p i• of a node i ∈ V is updated using only membership information from its neighbors in the graph. Thus, our algorithm can easily scale up to large datasets with large number of clusters, which is not the case of the algorithms listed above.

As proposed in [START_REF] Whang | Nonexhaustive, overlapping k-means[END_REF], the problem of soft graph clustering can also be tackled by a combination of spectral embedding [START_REF] Fowlkes | Spectral grouping using the nystrom method[END_REF][START_REF] Luxburg | A tutorial on spectral clustering[END_REF] and soft-clustering in the vector space, with algorithms such as the fuzzy C-means algorithm [START_REF] Dunn | A fuzzy relative of the isodata process and its use in detecting compact well-separated clusters[END_REF] or the NEO k-means algorithm [START_REF] Whang | Nonexhaustive, overlapping k-means[END_REF]. However, the number of clusters K has to be specified for these techniques. Besides, spectral embedding methods do not scale well to large size graphs [START_REF] Perozzi | Deepwalk: Online learning of social representations[END_REF]. For instance, they are unable to handle the subgraph of Wikipedia that we use in our experiments in Section 7. Finally, it is worth mentioning that several algorithms have been proposed for the related problem of overlapping community detection in networks [START_REF] Lancichinetti | Finding statistically significant communities in networks[END_REF][START_REF] Yang | Overlapping community detection at scale: a nonnegative matrix factorization approach[END_REF][START_REF] Zhang | Detecting overlapping communities in networks using spectral methods[END_REF]. Whereas these approaches return clusters C ⊂ V , soft-clustering methods give for each node the degrees of membership to all clusters, which represents a richer information. Note that some overlapping community detection algorithms also allow soft clustering [START_REF] Latouche | Overlapping stochastic block models[END_REF]. In the survey of [START_REF] Xie | Overlapping community detection in networks: The state-of-the-art and comparative study[END_REF], our algorithm would fall into the class of local expansion algorithms, where all algorithms have worst case complexity O(n 2 ) (except iCLD dealing with dynamic networks).

Soft modularity

Given an ordered partition P = (C 1 , . . . , C K ) of the nodes of V (i.e. a partition of V whose sets are given in a certain order), we define the membership matrix p associated with P as follows: ∀i ∈ V, ∀k ∈ 1, K , p ik = 1 if i ∈ C k , and 0 otherwise. It is easy to see that P is an optimal solution of the modularity maximization problem (2) if and only if its associated membership matrix p is an optimal solution of maximize

p∈Z n×n 1 w i,j∈V n k=0 W ij - w i w j w p ik p jk subject to ∀i ∈ V, n k=1 p ik = 1, ∀k ∈ 1, n , p ik ≥ 0. (3)
Given a membership matrix p, we use p i• to denote the vector that corresponds to the i th line of p and p •k to denote the vector that corresponds to its k th column. From the new formulation (3), we introduce the natural relaxation of the problem where coefficients of p can take real values maximize

p∈R n×n Q(p) subject to ∀i ∈ V, 1 T p i• = 1, p i• ≥ 0. (4)
where

Q(p) = 1 w i,j∈V W ij - wiwj w p i• T p j• .
Note that if p is a feasible solution of this problem, then p ik can be interpreted as a probability. It corresponds to the probability for a node i to belong to cluster k. So, by solving this relaxation of the modularity optimization problem we can obtain nodes that belong to multiple clusters unlike in the classical modularity maximization problem. We refer to a solution to this problem as a soft clustering of the nodes of G. We refer to this new objective function Q(p) as the soft modularity. Remark that, when

p i• ≥ 0, the constraint ∀i ∈ V, 1 T p i• = 1 is equivalent to the l 1 constraint p i• 1 = 1, where • 1 denotes the l 1 norm.

Alternating projected gradient descent

We can rewrite the relaxation of the modularity optimization problem as a classic minimization problem min p∈X J(p), where

J : p → -Q(p) is the cost func- tion and X = {p ∈ R n×n : ∀i, p i• ≥ 0, 1 T p i• = 1}
is the set of constraints. We define the normalized weighted Laplacian of the graph, as

L = I -D -1/2 W D -1/2
, where D is the diagonal matrix whose entries are the weighted degrees of the nodes: D = diag(w 1 , . . . , w n ).

The Laplacian matrix is symmetric and semi-definite positive [START_REF] Chung | Spectral graph theory[END_REF]. We use 0 ≤ λ 1 ≤ λ 2 ≤ . . . ≤ λ n to denote its eigenvalues.

Theorem 4.1. The function J is convex if and only if the second lowest eigenvalue λ 2 of the normalized Laplacian L verifies λ 2 ≥ 1.

Proof. We prove in the supplementary material that the hessian matrix H of J is congruent to a matrix whose eigenvalues are 0, λ 2 -1, . . . , λ n -1.

The condition λ 2 ≥ 1, which corresponds to a large spectral gap, is in general not satisfied. For instance, Lemma 1.7. in [START_REF] Chung | Spectral graph theory[END_REF] proves that λ 2 ≤ 1 if G is unweighted and not complete. Therefore, the loss function associated with our problem is in general non-convex. However, it is convex in p i• for a graph with no self-loop as shown in the following proposition. In the rest of the paper, we assume that the graph G does not contain any self-loop.

Proposition 4.2. If the graph does not contain any self loops, i.e. if W ii = 0 for all node i, the function p → J(p) is convex with respect to variable p i• for all i ∈ V .

Proof. If W ii = 0, the hessian matrix H i of J with respect to

p i• is H i = 2(w i /w) 2 I.
With this result, we can apply the methods of convex optimization with convex constraints to the problem min pi•∈Y J(p), with fixed p j• for j = i, where the set Y corresponds to the probability simplex of R n : Y = {q ∈ R n : q ≥ 0, 1 T q = 1}. We use π Y to denote the euclidean projection onto the probability simplex Y, π Y (x) = arg min y∈Y xy 2 . Then, using the projected gradient descent method [START_REF] Goldstein | Convex programming in hilbert space[END_REF][START_REF] Levitin | Constrained minimization methods[END_REF], we can define an update rule for variable p i• such that modularity is non-decreasing at each step. The gradient of J with respect to p i• is:

∇ i J = -2 w j∈V (W ij -w i w j /w)p j• . Theorem 4.3.
The soft modularity objective function Q(p) is non-decreasing under the update rule

p i• ← π Y   p i• + 2t w j∈V W ij - w i w j w p j•   ( 5
)
for all node i ∈ V if the step size t verifies t < (w/w i ) 2 . Q(p) is invariant under all these update rules iff p i• maximizes Q(p) with fixed p j• , j = i, for all node i ∈ V .
Proof. In the supplementary material, we prove that 2 where p + is the matrix p after the update, and

Q(p + ) ≥ Q(p) + t 1 -t wi w 2 G i (p)
G i (p) = (p i• - p + i• )/t.
We now consider the natural algorithm that cycles through the nodes of the graph to apply the update (5). Theorem 4.4. The algorithm converges to a local maximum of the soft modularity function p → Q(p) which is a fixed point of the updates of [START_REF] Chung | Spectral graph theory[END_REF].

Proof. This is a direct consequence of the analysis conducted in the proof of Theorem 4.3 provided in the supplementary material.

Soft clustering algorithm

In the previous section, we have presented updates rules that can be applied in an alternating fashion to find a local maximum of the soft modularity Q(p). However, there are three major issues with a naive implementation of this method. First, the gradient descent update for each node i is computationally expensive because the update rule (5) involves a sum over all nodes of V . Then, the size of a solution p is n 2 , which is prohibitive for large datasets. Finally, the computational cost of the projection π Y onto the probability simplex is classically in O(n log n) [START_REF] Duchi | Efficient projections onto the l 1-ball for learning in high dimensions[END_REF], which can be a limiting factor in practice.

In the present section, we present an efficient implementation of the algorithm that solves these three problems.

In particular, our implementation guarantees that the update step for node i only requires local computation, is memory efficient, and uses a fast mechanism for projection.

Local gradient descent step

The update of Theorem 4.3 relative to node i ∈ V can be decomposed into two steps:

1. pi• ← p i• + 2t w j∈V W ij - wiwj w p j• 2. p i• ← π Y ( pi• ).
At first sight, step 1 seems to depend on information from all nodes of V , and to require the computation of a sum over n terms. However, the gradient descent step can be written as

pi• ← p i• + 2t w j∼i W ij (p j• -p),
where p is weighted average of vector p j• : p = j∈V wj w p j• . The vector p ∈ R n can be stored and updated throughout the algorithm. Then, the computation of pi• is purely local, in the sense that it only requires information from neighbors of node i. Moreover, the sum to compute contains only d i terms, where d i = |Nei(i)| is the unweighted degree of node i. In most real-life graphs, the degree distribution follows a heavy tail distribution and most of nodes have a low degree d i .

Cluster membership representation

The cluster membership variable p is a n×n matrix. In a naive implementation, the memory cost of storing p is therefore in O(n 2 ). However, we will see below that the projection step can be written π Y ( pi• ) = ( pi• -θ1) + with θ ≥ 0, where [(x) + ] k = max(x k , 0). Thus, the projection acts as a thresholding step. We expect it to have the same effect as a Lasso regularization 1 and to lead to a sparse vector p i• . To take benefit of this sparsity, we only store the non-zero coefficients of p. If, for all node i ∈ V , the number of non-zero values in vector p i• is lower than a given L, then the memory cost of storing p is in O(nL). We will see, in our experiments in Section 7, that the average number of positive components that we observe for vectors p i• throughout the execution of our algorithm is lower than 2, and that the maximum number of positive components that we record is L = 65 for a graph with 731, 293 nodes. 1 The l1 proximal operator is x → sign(x)(|x| -λ)+.

Projection onto the probability simplex

The classic algorithm to perform the projection π Y ( pi• ) onto Y is described in [START_REF] Duchi | Efficient projections onto the l 1-ball for learning in high dimensions[END_REF]. It starts by sorting the vector pi• into µ:

µ 1 ≥ µ 2 ≥ • • • ≥ µ n . Then, it finds ρ = max j ∈ [n], µ j -1 j j r=1 µ r -1 > 0 , and defines θ = 1 ρ ( ρ i=1 µ i -1) so that π Y ( pi• ) = ( pi• -θ1) + .
This algorithm runs in O(n log n) because of the sorting step.

In our case, the complexity can be reduced to

O(L i log L i ), with L i = |supp i (p)|, where supp i (p) = {k ∈ 1, n , ∃j ∈ Nei(i) ∪ {i}, p ik > 0}
is the union of the supports of vectors p j• for j ∈ Nei(i) ∪ {i}. L i is upper-bounded by L(d i + 1) which is typically much smaller than the total number of nodes n. Proof. We prove this result in the supplementary material.

MODSOFT

Finally, the algorithm can be described as follows.

• Initialization: p ← I and p ← w/w. • One epoch:

For each node i ∈ V , -∀k ∈ supp i (p), pik ← p ik + t j∼i W ij (p jk -pk ) -p + i• ← project( pi• ) -p ← p + (w i /w)(p + i• -p i• ) and p i• ← p + i• .
where w is the vector (w i ) 1≤i≤n , project is the adapof the algorithm of [START_REF] Duchi | Efficient projections onto the l 1-ball for learning in high dimensions[END_REF] presented above, and t is the effective learning rate t = 2t/w. One epoch of the algorithm is typically repeated until the increase of modularity falls below a given threshold > 0, as in the Louvain algorithm. Note that we chose, in this implementation, to put each node in its own cluster at initialization, but our algorithm could also be initialized with the results of another algorithm (e.g. Louvain). We refer to this algorithm as MODSOFT (MODularity SOFT clustering).

Link with the Louvain algorithm

In this section, we show that if the learning rate t is larger than a graph-specific threshold, one update performed by our algorithm reduces to a local update performed by the Louvain algorithm. First, we observe that, if the membership matrix p verifies p ∈ {0, 1} n 2 , then the update rule (5) for a node i ∈ V becomes

p i• ← π Y 1 {i∈C k } + 2t w w i (C k ) -w i Vol(C k ) w k
where C k = {i ∈ V : p ik = 1} refers to the k th cluster defined by p, w i (C) denotes the degree of node i in cluster C ⊂ V , w i (C) = j∈C W ij , and Vol(C) denotes the volume of cluster C, Vol(C) = j∈C w j .

In the following, we assume that: ∀C, C ⊂ V, ∀i ∈ V ,

C = C ⇒ w i (C)-wiVol(C) w = w i (C )-wiVol(C ) w . (
C,C ⊂V i∈V C =C w i (C)-wiVol(C) w -w i (C )-wiVol(C ) w , then the update rule (5) for node i ∈ V reduces to p ik ← 1 if k = arg max l:j∈C l ,j∼i w i (C l ) -w i Vol(C l ) w
, and p ik ← 0 otherwise, where C k denotes the k th cluster defined by p.

Proof. We show in the supplementary material that assuming that the updated vector p i• has more than one positive component leads to a contradiction with t > w/δ.

In particular, this result shows that if p ∈ {0, 1} n 2 and t is large enough, then for each update of our algorithm, the updated membership matrix p + verifies p + ∈ {0, 1} n 2 , which corresponds to the maximum sparsity that can be achieved by a feasible solution p ∈ X .

The Louvain algorithm [START_REF] Blondel | Fast unfolding of communities in large networks[END_REF] is a popular heuristic to find a local maximum for the hard modularity maximization problem [START_REF] Blondel | Fast unfolding of communities in large networks[END_REF]. One epoch of the main routine of the Louvain algorithm considers successively each node i ∈ V , and (1) removes i from its current cluster, (2) applies an update rule, that we refer to as LouvainUpdate(i), which consists in transferring i to the cluster that leads to the largest increase in hard modularity. Proposition 6.2 gives an explicit formula to pick the cluster C k * to which node i is transferred by LouvainUpdate(i).

w i (C k ) -w i Vol(C k ) w
. Proof. We derive this result from the definition of the hard modularity in the supplementary material.

Using Proposition 6.1, we see that the update performed by LouvainUpdate(i) and the update rule ( 5) are equivalent if (H) is verified and t > w/δ.

7 Experimental results

Synthetic data

We consider a variant with overlaps of the popular Stochastic Block Model (SBM) [START_REF] Holland | Stochastic blockmodels: First steps[END_REF], defined as follows.

We have

V = C 1 ∪ • • • ∪ C K with |C 1 | = . . . = |C K | = c and |C k ∩ C k+1 | = o for all k.
For all nodes i and j such that i = j, the edge (i, j) exists with probability p in if i, j ∈ C k for some k, and with probability p out otherwise. In order to numerically evaluate the performance of our algorithm on this model, we first compute the relative difference between the modularity Q soft obtained with our algorithm and the modularity Q hard obtained with the Louvain algorithm, ∆Q = (Q soft -Q hard )/Q hard . We also compare the planted clusters of our model C k with the clusters Ĉk returned by our algorithm. We consider that these clusters are the non-empty sets Ĉk where Ĉk = {i ∈ V : p ik > 0}. Note that the clusters Ĉk returned by our algorithm can overlap just as the C k .

In order to compare the clusters C k and Ĉk , we use the average of the classic F1 score defined as the harmonic mean of the precision and the recall, which takes values between 0 and 1, and is equal to 1 if the estimated clusters Ĉk correspond exactly to the initial clusters C k . We compute the same two metrics for Chang's algorithm [START_REF] Chang | A unified framework for sampling, clustering and embedding data points in semi-metric spaces[END_REF]. Since the membership matrix returned by this algorithm is not sparse, we define Ĉk as Ĉk = {i ∈ V : p ik > α} for the F1 score evaluation. We present the results for α = 0.05, but other values of α leads to similar results.

For different values of the cluster size c, we run our algorithm, Chang's algorithm and the Louvain algorithm on 100 random instances of this model, with o = 2, K = 2, p in = 0.9, and p out = 0.1, We display the results in Figure 1. The average F1 scores for our algorithm and Chang's algorithm are close to 1 (> 0.99 on average), and are higher than those obtained by Louvain (0.94 on average). Besides, the relative difference to the modularity found by Louvain ∆Q for Chang's algorithm and our algorithm is always positive. However, we note that ∆Q and the differences between the F1 scores become negligible as the cluster size c increases. This can be explained by the fact that, as the cluster size c increases, the proportion of nodes with mixed membership, i.e. nodes that belong to multiple clusters C k , decreases. If this proportion is low, the misclassification of nodes with mixed membership has a low impact on the modularity score Q and on the F1 score.

We also represent in Figure 1 the sparsity of the solution in term of the number of non-zero coefficients of the matrix p, and the average execution time in milliseconds for the different algorithms. We observe that our algorithm provides performance comparable to Louvain (but with richer membership information). Besides, we see that MODSOFT runs more than 10 times faster and uses 5 times less memory than Chang's algorithm.

Real data

Wikipedia

We first consider a graph built from Wikipedia data, where nodes correspond to Wikipedia articles and edges correspond to hyperlinks between these articles (we consider that the edges are undirected and unweighted). Wikipedia articles can correspond to a wide variety of entities: persons, locations, organizations, concepts etc. In order to restrict ourselves to homogeneous entities, we choose to limit ourselves to Wikipedia pages that correspond to humans. For this purpose, we use Wikidata, a knowledge base that provides structured data about Wikipedia articles. In particular, Wikidata objects have an instance_of field that we use to determine if a Wikipedia article represents a human or not. The subgraph induced by the human entities has 731,293 nodes and 3,266,258 edges.

We run the Louvain algorithm, Chang's algorithm and MODSOFT on this subgraph. As described in Section Our algorithm leads to a higher modularity than the Louvain algorithm on this dataset. However, this increase represents a modularity increase ∆Q of slightly less than 1%. This can be explained by the fact that the nodes with mixed membership found by our algorithm represent less than 5% of the total number of nodes. Even if the performance increase brought by our algorithm in terms of modularity is marginal, we qualitatively observe that the results of our algorithm, and in particular the nodes with mixed membership identified by our approach, are particularly relevant.

In order to be able to graphically represent the results of our algorithm, we consider the subgraph induced by the top 100 nodes in term of degree. We display the results of our algorithm on this smaller graph in Figure 2. In particular, we observe that Napoleon, who became Emperor of the French in 1804, appears to belong to the same cluster as European leaders from the old regime, such as Charles V and Louis XIV, and, at the same time, to the same cluster as post-Frenchrevolution leaders such as Queen Victoria and Edward VII. We also see that Miles Davis, a famous American jazz trumpeter, is found to belong to the same cluster as classical music composers, such as Claude Debussy and Johan Sebastian Bach, but also to the same cluster as pop/rock musicians such as Jimi Hendrix and Michael Jackson. Finally, we observe that Albert Einstein who became politically involved during World War II, belongs to the same cluster as thinkers, intellectuals and scientists, such as Sigmond Freud and Plato, and to the same cluster as political leaders during World War II such as Adolf Hitler and Joseph Stalin.

Open Flights

Then, we consider a graph built from the Open Flights database that regroups open-source information on airports and airlines. The graph we consider is built as follows: the nodes correspond to the airports and the edges correspond to the airline routes, i.e. there is an edge between two airports if an airline operates services between these airports.

We apply our algorithm and the Louvain algorithm to this graph. During the execution of the algorithm, the number of positive components of vectors p i• remains lower than 11 and is on average 1.17. We measure a relative increase in modularity ∆Q of less that 1%, and yet we find that our algorithm brings insightful information about the dataset by identifying bridges between important zones of the world. In Figure 3, we display the results of our algorithm. We use plain colors to code for the cluster membership of pure nodes, i.e. nodes that belong to a unique cluster, and bigger Note that most of these nodes with mixed membership are located on the borders between large regions of the world. For instance, we remark that numbers of nodes are located on the frontier between Europe, Africa and Middle East, such as the Dubai airport in the United Arab Emirates. We also observe that airports with mixed membership are located all along the border between the United States on the one hand, and Canada and Alaska on the other hand.

SNAP datasets

We present further experiments on datasets from SNAP that come with ground-truth clusters [START_REF] Yang | Defining and evaluating network communities based on ground-truth[END_REF] (see Table 1). The Amazon dataset corresponds to a product co-purchasing network. Its nodes represent Amazon products and its edges correspond to frequently copurchased products. The ground-truth clusters are defined as the product categories. The DBLP dataset corresponds to a scientific collaboration network. The nodes represent the authors and the edges the coauthorship relations. The scientific conferences are used as ground-truth clusters. In the YouTube dataset, nodes represent users and edges connect users who have a friendship relation. User groups are used as ground-truth clusters. We use Louvain and the stateof-the-art algorithm OSLOM [START_REF] Lancichinetti | Finding statistically significant communities in networks[END_REF] as baselines for our experiments. OSLOM is known as one of the most effective overlapping community detection algorithm [START_REF] Xie | Overlapping community detection in networks: The state-of-the-art and comparative study[END_REF]. Note that Chang's algorithm does not run on these large graphs. We use a Cython implementation of the Modsoft algorithm that is publicly available 2 . In Table 2, we give the F 1 -scores and the execution times for Louvain (L), OSLOM (O) and Modsoft (M). The experiments were carried out on an AWS instance with 2 https://github.com/ahollocou/modsoft 64GB RAM and 16 vCPUs with Intel Xeon. A learning of t = 0.1 was used for Modsoft. Note that we did not include the results of OSLOM on the YouTube dataset because its the execution time exceeded 6 hours. We see that Modsoft outperforms OSLOM both in terms of F 1 -score (11% higher on average) and execution time (OSLOM is about 40 times slower). Note that, in our experiments, the objective function of Modsoft is within 1% of its optimal value after 3 to 5 epochs. Moreover, although Modsoft is slower than Louvain (about twice slower on average), it offers F 1 -scores that are at least twice as large as the ones provided by Louvain. This can be explained by the fact that the clusters returned by Louvain are necessary disjoint whereas the ground-truth clusters in all three datasets can overlap. In contrast, Modsoft is able to capture these mixed memberships. Besides, we observe in these experiments that, as expected, the larger t, the sparser the solution. For instance, on the Amazon dataset, the average number of posistive components per row is 5.58 for t = 0.5, 2.27 for t = 1, and 1.87 for t = 2.

Conclusion

We studied a relaxation of the popular modularity maximization problem, with the aim of performing soft clustering of graph nodes instead of hard partitioning. In order to efficiently solve this relaxation, we introduced a novel algorithm that is based on an alternating projected gradient descent. By diving into the specific form of the gradient descent updates, we were able to guarantee the locality of each algorithm step, and to make good use of the sparsity of the solutions. As a result, unlike existing methods, our approach is both local and memory efficient. Furthermore, we proved that our algorithm includes the main routine of the popular Louvain algorithm for t > w/δ. In this case, Modsoft outputs the sparsest possible solution, with only one component per row. We illustrated on real-life examples, that our algorithm has an execution time and a memory consumption comparable to the Louvain algorithm, but goes further than Louvain by identifying nodes with mixed memberships that represent important bridges between clusters, even if this does not always translate into a large modularity increase.

A generalization of standard modularity introducing a resolution parameter has been proposed by different authors [START_REF] Lambiotte | Laplacian dynamics and multiscale modular structure in networks[END_REF][START_REF] Reichardt | Statistical mechanics of community detection[END_REF] to perform multi-level clustering. In future works, we could add this resolution parameter to the soft modularity function in order to perform soft graph clustering at different scale. In the future, we would also like to compare our approach to a nonnegative matrix factorization method applied to the matrix W , which can be performed with a comparable alternating projected gradient descent [START_REF] Lin | Projected gradient methods for nonnegative matrix factorization[END_REF].

  Proceedings of the 22nd International Conference on Artificial Intelligence and Statistics (AISTATS) 2019, Naha, Okinawa, Japan. PMLR: Volume 89. Copyright 2019 by the author(s).

Proposition 5 . 1 .

 51 In order to compute π Y ( pi• ), we only need to sort the components k ∈ supp i (p) of pi• to determine ρ and θ. All components k / ∈ supp i (p) of π Y ( pi• ) are set to zero.

Figure 1 :

 1 Figure 1: Algorithm performance in function of the cluster size on the overlapping stochastic block model

Figure 2 :

 2 Figure 2: MODSOFT on a small subgraph of Wikipedia

Figure 3 :

 3 Figure 3: MODSOFT on the OpenFlights graph (red circles represent nodes with mixed-membership)

Table 1 :

 1 F 1 -scores and execution times in seconds on SNAP datasets

				Graph size
				|V |		|E|
		Amazon	334,863	925,872
		DBLP	317,080 1,049,866
		YouTube 1,134,890 2,987,624
		F 1 -score		Execution time (s)
		L	O	M	L	O	M
	Amazon	0.28 0.47 0.53 15 1038	32
	DBLP	0.14 0.35 0.38 17 1717	33
	YouTube 0.01	-0.15 45	-	98

Table 2 :

 2 F 1 -scores and execution times in seconds on SNAP datasets red dots to identify nodes with a mixed membership.


Khasim Cali

Emmanuel Scorsone

Krishna Persaud

ODORANT BINDING PROTEINS BASED SNIFFING DEVICE FOR DETECTION OF TOBACCO

Keywords: Odorant Binding Proteins, Surface Acoustic Wave Sensors, Sensor Array, Tobacco detection

Affinity constants for a selection of Odorant Binding Proteins (OBPs) from mammalian and insect sources were determined toward volatile signature compounds from tobacco in the liquid phase using a fluorescent probe and competitive ligand binding. The most sensitive proteins were selected and these were immobilised on to diamond coated Surface Acoustic Wave (SAW) sensors to construct an array of eight sensors to detect tobacco volatiles. The system was then tested for its ability to detect vapours of tobacco itself as well as other markers. The results were positive, a pattern of responses across the array was observed as each sensor gave a specific response towards different analytes. To determine the stability of these OBP biosensors over time, the sensors were repeatedly exposed to pulses of saturated D(+)-carvone vapour over three months. Little degradation was observed and these sensors were able to sensitively detect the target analytes after this period.

INTRODUCTION

Odorant Binding Proteins (OBPs) are small soluble proteins present in the olfactory systems of vertebrates and insects having affinity to a great number of organic molecules (odorants and pheromones). They are exceptionally stable to temperature and have a compact folding structure that makes these proteins resistant to proteolytic degradation. These characteristics are attractive to build robust, reliable and inexpensive biosensors for different applications.

Recently it has been shown that Odorant Binding Proteins (OBPs) can be used as Biosensors to detect variety of volatile organic compounds [START_REF] Cali | Development of mutant Odorant Binding Proteins (OBPs) for Security Applications: Detection of Explosives and Drugs[END_REF][START_REF] Cali | Odorant Binding Proteins Biosensors for detection of drugs[END_REF][START_REF] Manai | Grafting odorant binding proteins on diamond bio-MEMS[END_REF]. Recombinant porcine OBPs were chemically grafted on polycrystalline diamond surfaces for biosensor development; the resulting sensors were able to detect 2, 4-DNT in vapour form [START_REF] Manai | Grafting odorant binding proteins on diamond bio-MEMS[END_REF]. The aim of this work is to develop a Sniffing device based on Odorant Binding Proteins Biosensors for the detection of tobacco. Eleven OBPs were selected for this work. Two are porcine Wild Type OBP1 (pWTOBP1) and its mutant protein pF88W. Nine are from insect sources, namely wasp 1PDOMOBP1; six wild type proteins from mosquito Anopheles gambiae, AgamOBP1, AgamOBP4, AgamOBP5, AgamOBP19, AgamOBP21 and AgamOBP47; while the remaining two are mutants of AgamOBP1 developed using in silico mutagenesis and docking screening techniques, namely AgamOBP1_S82P and AgamOBP1_K93H [START_REF] Cali | Development of mutant Odorant Binding Proteins (OBPs) for Security Applications: Detection of Explosives and Drugs[END_REF][START_REF] Cali | Odorant Binding Proteins Biosensors for detection of drugs[END_REF].

METHODS

The SAW readout device was from Karlsruhe Institute of Technology (KIT), Germany. Intrinsic volatile signature compounds from tobacco were selected [START_REF] Coleman | Tobacco volatiles: Gas Chromatography, Encyclopedia of Separation Science[END_REF]. These were 2methylpropanal, 2-butanone, ethylpyrazine and acetone. Affinity constants of the screened OBPs toward the above analytes were initially determined in the liquid phase using a fluorescent probe and competitive ligand binding [START_REF] Cali | Development of mutant Odorant Binding Proteins (OBPs) for Security Applications: Detection of Explosives and Drugs[END_REF][START_REF] Cali | Odorant Binding Proteins Biosensors for detection of drugs[END_REF]. The best performing proteins were selected and these were immobilised on to diamond coated Surface Acoustic Wave (SAW) sensors to form an array of eight SAW Biosensors which was then mounted on to the SAW read-out platform (SAW DEVICE) [START_REF] Manai | Grafting odorant binding proteins on diamond bio-MEMS[END_REF]. The system was then tested for its ability to detect vapours of the selected markers plus three different brands of tobacco. The tobacco brands chosen were GOLDEN VIRGINIA (Imperial Brands/UK), PALL MALL (British American Tobacco Group/UK) and HOLBORN YELLOW (Japan Tobacco International/UK). For each sample, measurements were repeated fifty times over a period of 1 week. For stability studies the sensors were repeatedly exposed to pulses of saturated D (+)carvone vapour over three months.

RESULTS

Determination of Binding affinities in solution

The data in Figure 1 shows that OBPs can bind to these four tobacco volatiles and with differing affinity profiles for each of the analyte. The poorly responding proteins AgamOBP5, AgamOBP21, and AgamOBP47 were eliminated and the rest were selected for incorporation into an array of SAW Biosensors. Vapours of tobacco as well as marker volatiles compounds were presented to the device, and results for tobacco samples are presented in Figures 2 and3.

Stability of the biosensors over time

Despite being reported that Odorant Binding Proteins (OBPs) can be used as Biosensors to detect variety of volatile organic compounds [START_REF] Cali | Development of mutant Odorant Binding Proteins (OBPs) for Security Applications: Detection of Explosives and Drugs[END_REF][START_REF] Cali | Odorant Binding Proteins Biosensors for detection of drugs[END_REF][START_REF] Manai | Grafting odorant binding proteins on diamond bio-MEMS[END_REF], it was not known how stable these sensors are over time. The SAW biosensor array was exposed repeatedly to pulses of saturated D(+)-carvone vapour over three months (Figure 4). Carvone was used as a control analyte, because several OBPs are very sensitive to this compound [START_REF] Mulla | Capacitance-modulated transistor detects odorant binding protein chiral interactions[END_REF].

DISCUSSION

The data in Figure 2 and 3 shows that the sensors are sensitive to the tobacco vapour, each sensor displayed different response intensity to the same concentration of vapour, giving rise to a specific response pattern of the array towards the tobaccos. For example AgamOBP1 sensor is the best performing having highest sensitivity followed by AgamOBP19 while sensors pWTOBP1 and AgamOBP1_K93H have the lowest sensitivity.

The data in Figure 4 shows that the sensors retain their ability to detect analytes for more than three months. This is important and a key aspect of using these proteins as biosensors. Generally, the sensitivity of each sensor increases after sensor fabrication (Day 1 to Day 35) before achieving a stable sensitivity after some days (from day 56 to day 104). There are day to day variations in the intensity of the response -however this is largely due to the complexity of generating a consistent headspace for complex materials such as tobacco. However, what is important is that the sensors continue to operate over long periods of time and they do not fail with repeated exposure to volatile analytes. Also the data in Figure 2, 3 and 4 shows that the biosensors may be able to recognize the target analytes selectively, the response profiles between Tobaccos (Figures 2 and3) are quite different to that for carvone (Figure 4). For example while AgamOBP1 sensor is the best for tobacco, for carvone 1PDOMOBP1 and AgamOBP1_S82P sensors are the best performing compared to all other sensors.

CONCLUSIONS

A robust SAW Sniffing device made up of an array of Odorant Binding Proteins is reported. This opens up a new approach for detecting tobacco and other potential targets like drugs, explosives or other analytes of interest.

FUTURE WORK

The end goal is to integrate this SAW Biosensor device into a detection toolbox connected with the sampling system and supplied with operating electronics and pattern recognition software to be used for the detection of tobacco, and other potential targets for customs and security applications.

Figure 1 :Figure 2 :

 12 Figure 1: Affinity constants of the screened OBPs toward the selected signature volatiles of tobacco, determined in the liquid phase using a fluorescent probe and competitive ligand binding.3.2. Detection of tobacco by the SAW biosensor array

Figure 3 :

 3 Figure 3: Repeated response of the SAW biosensor array to the saturated tobaccos vapour. For each tobacco brand, fifty measurements were taken, as an example six repetitions are shown in this figure, where (a) GOLDEN VIRGINIA (SMOOTH YELLOW), (b) PALL MALL, and (c) HOLBORN YELLOW.

Figure 4 :

 4 Figure 4: Biosensor stability, measured by repeated responses to D(+)-carvone vapour measured by the SAW biosensor array over three months.

Acknowledgements:

The SNIFFER project European Community (FP7/2007-2013) n°285203 and project n°653323 -C-BORD (H2020).


T Iwanowska

M Kosinski

A Swiderski

D Syntfeld-Kazuch

Wolski

Comparison of prompt and delayed photofission neutron detection techniques using different types of radiation detectors

Pawel Sibczynski, Andrzej Dziedzic, Krystian Grodzicki, Joanna Iwanowska, Tymoteusz Kosiński, Michał Matusiak, Marek Moszyński, Lukasz Swiderski, Agnieszka Syntfeld-Każuch, Dariusz Wolski, Frédérick Carrel, Amélie Grabowski, Matthieu Hamel, Frederic Lainé, Adrien Sari, Alessandro Iovene, Carlo Tintori

Abstract-For several years, detection of various threats on country borders plays a significant role in the frame of Homeland Security applications. One of this threat is the illicit trafficking of nuclear materials (especially including Special Nuclear Material -SNM -235 U, 233 U or 239 Pu), which can be potentially used for production of nuclear weapon as well as radiological dispersal device (RDD) -known also as a "dirty bomb". In order to detect the potentially hidden nuclear material, systems using linear accelerators and a group of detectors are developed by several scientific groups around the world. Besides solutions focusing on detection of delayed γ-rays or neutrons, also the systems dedicated for prompt neutron detection were proposed. One of the possible prompt neutron detection technique is known as Threshold Activation Detection (TAD). This technique relies on activation of 19 F nuclei in the scintillator medium by fast neutrons and registration of high-energy β particles and γ-rays from the decay of reaction products (for example, 19 F(n,α) 16 N or 19 F(n,p) 19 O). Recent studies in the frame of the European Horizon 2020 C-BORD project showed that, despite the low 19 F(n,α) 16 N or 19 F(n,p) 19 O reaction cross-section, the method could be a good solution for detection of shielded nuclear material. A benchmark of the TAD technique based on fluorine detectors with reference method focused on delayed neutron detection with 3 He detectors will be presented in this paper. These experimental results were obtained using 9 MeV Varian Linatron M9 linear accelerator (LINAC).

Index Terms-photofission, prompt neutrons, BaF2, delayed neutrons, threshold activation detection, TAD, LINAC, special nuclear materials, SNM, scintillators.

I. INTRODUCTION

T HE aim of the study is to carry out a comparative study of nuclear material detection techniques basing on prompt and delayed neutron registration. The first technique, This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 653323. Disclaimer: This text reflects only the author's views and the Commission is not liable for any use that may be made of the information contained therein.

Pawel Sibczynski, Andrzej Dziedzic, Krystian Grodzicki, Joanna Iwanowska, Tymoteusz Kosiński, Michał Matusiak, Marek Moszyński, Lukasz Swiderski, Agnieszka Syntfeld-Każuch, Dariusz Wolski are with the National Centre for Nuclear Research (NCBJ), Otwock, 05-400 Poland (telephone: +48 22 273 1440, e-mail: pawel.sibczynski@ncbj.gov.pl, pawel.sibczynski@gmail.com) Frédérick Carrel, Amélie Grabowski, Matthieu Hamel, Frederic Lainé, Adrien Sari are with CEA, LIST, Sensors and Electronic Architectures Laboratory, F-91191 Gif-sur-Yvette, France.

Alessandro Iovene, Carlo Tintori are with the CAEN S.p.A., Via Vetraia 11, 55049 -Viareggio (LU) -Italy. applied for prompt neutron detection, is based on Threshold Activation Detection (TAD) using scintillators containing 19 F. The incident fast neutrons activate the 19 F nuclei, resulting in (n,α) or (n,p) reaction and emission of high energy β particles (up to 10.4 MeV) and γ-rays (in example 6.1 MeV). Potential scintillators of interest for this application are liquid fluorocarbon based EJ-313 [START_REF] Gozani | Neutron threshold activation detectors (TAD) for the detection of fissions[END_REF], [START_REF] Stevenson | Linac based photofission inspection system employing novel detection concepts[END_REF], [START_REF] Gozani | The investigation of fast neutron Threshold Activation Detectors (TAD)[END_REF] or BC-508, BaF 2 [START_REF] Sibczynski | Verification of threshold activation detection (TAD) technique in prompt fission neutron detection using scintillators containing 19 F[END_REF] or fluorine based plastic scintillators based on pentafluorostyrene [START_REF] Hamel | A fluorocarbon plastic scintillator for neutron detection: Proof of concept[END_REF]. The second technology is focused on the detection of delayed neutrons using 3 He detectors [START_REF] Loaiza | High-efficiency 3 He proportional counter for the detection of delayed neutrons[END_REF], having a very low sensitivity to γ-rays and high cross section for thermal neutrons of 5530 barns. In this paper we investigated the performance of 5" × 1" BaF 2 scintillator and compared with that for 3 He detector, Canberra 150NH100 model [START_REF] Carrel | New experimental results on the cumulative yields from thermal fission of and and from photofission of and induced by bremsstrahlung[END_REF]. Particularly, BaF 2 scintillators present many advantages comparing to fluorocarbon liquid scintillators, such as non-flammability and lack of hydrogen in the scintillator medium, decreasing the efficiency of threshold reaction due to neutron scattering. In addition, the risk of leakage of the liquid fluorocarbon scintillators is a significant disadvantage, promoting the BaF 2 scintillator as an interesting alternative.

II. RESULTS

A. Measurements of the BaF 2 response to the 252 Cf

The BaF 2 scintillator was coupled to an Electron Tubes ET9390 series 130 mm photomultiplier with silicone grease in order to improve light transportation to a photocathode. Finally, detector was mounted into laboratory-made housing designed and manufactured at NCBJ. The test site at NCBJ, where the detector was exposed to neutrons from a 252 Cf source is presented in Fig. 1.

Primarily, the response of the BaF 2 to neutrons was measured for 252 Cf neutron source at NCBJ, emitting 30000 neutrons per second in 4π. Source was shielded with 5 cm thick Pb brick in order to absorb high-energy prompt and delayed γ-rays. The 252 Cf was placed 15 cm from the front of the detector. Anode pulses from the photomultiplier were recorded by CAEN V1730 VME-type digitizer offering sampling rate of 500 MS/s and accuracy of 14 bits. Pulses were integrated by a DPP-PSD firmware implemented to the digitizer FPGA. The net spectrum recorded after 4 h was showed in Fig. 2. It can be noticed that fast neutron detection efficiency of the 5" × 1" BaF 2 based photodetector is significantly better than that for the detector presented in [START_REF] Gozani | The investigation of fast neutron Threshold Activation Detectors (TAD)[END_REF].

B. Tests with 9 MeV LINAC at SAPHIR facility

After examination of the detector efficiency at NCBJ, the performance of depleted uranium (DU) detection were investigated at SAPHIR facility in CEA LIST Saclay, France. The test site is shown on Fig. 3. We used the commercially produced Varian Linatron M9 LINAC, emitting bremsstrahlung radiation with endpoint at 9 MeV, in order to induce photofission in DU samples. Both BaF 2 and 3 He detectors were placed one meter from DU samples, which in turn were placed one meter from the LINAC tungsten conversion target. Measurements were performed for shielded and unshielded DU. A 5 cm thick Pb brick was used as the shielding. The Teflon in front of the BaF 2 detector was used as an additional TAD material containing 19 F, emitting 6.1 MeV γ-rays. During cargo inspection it is convenient to provide as low dose as reasonable possible due to potential exposition of illegal human transportation to radiation. Thus, the LINAC beam frequency was set to 40 Hz. Short time of inspection was introduced; irradiation lasted 10 s, then after 1 s of cooling time, the acquisition with the BaF 2 lasted 20 s. Spectra were acquired three times, then results were summed. Summarizing, the total time of acquisition lasted 93 s. In the case of 3 He detector, after single 10 s of irradiation, the measurement lasted 60 s.

Fig. 4 and Fig. 5 show the obtained net spectra for the BaF 2 and 3 He detector. The results expressed in the net number of counts were presented in Table I. The sensitivity of the TAD technique is definitely lower than the detection of delayed neutrons after irradiation with LINAC. The decrease of the Fig. 4. Net spectra obtained from the BaF 2 detector after three irradiations and counts summing. Data were binned with step of 500 keV. sensitivity of shielded DU detection is rather similar for these two techniques and varies between 20 -35%. However, the worldwide shortage of 3 He caused significant price increase of 3 He gas detectors. Thus, the BaF 2 detector can be considered as an alternative solution offered for considerable price.

III. CONCLUSION

The prompt neutron detection efficiency of the BaF 2 scintillator is significantly lower in comparison with the sensitivity of the 1 m long 3 He detector, used for delayed neutron detection. However, it is important to point out that the 3 He detector is much more expensive than 5" × 1" BaF 2 detector used in the present study. Thus, comparison tests with alternative delayed neutron detectors as well as detector geometry optimization will be continued.

Fig. 1 .

 1 Fig. 1. Test of the prototype BaF 2 detector designed at NCBJ with 252 Cf neutron source.

Fig. 2 .

 2 Fig. 2. The BaF 2 detector response to 252 Cf neutron source. The β continuum up to 10.4 MeV is clearly seen.

Fig. 3 .

 3 Fig. 3. Test site at SAPHIR facility, CEA LIST, France.

Fig. 5 .

 5 Fig. 5. The net counts registered with the 3 He detector in the time scale limit set to 60 s.

TABLE I NUMBER

 I OF NET COUNTS REGISTERED WITH THE BAF 2 AND 3 HE

		DETECTORS.
	Detector	Size	Counts
	BaF 2 + 5cm Pb	5" × 1"	189 ± 14
	BaF 2	5" × 1"	275 ± 17

3

He 150NH100 + 5cm Pb 1 m × 1" 3600 ± 60 3 He 150NH100 1 m × 1" 4450 ± 67

ACKNOWLEDGMENT

This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 653323. Disclaimer: This text reflects only the author's views and the Commission is not liable for any use that may be made of the information contained therein.