Jules Pitcho
email: jules.pitcho@uzh.ch

Massimo Sorella
email: massimo.sorella@epfl.ch

ALMOST EVERYWHERE NON-UNIQUENESS OF INTEGRAL CURVES FOR DIVERGENCE-FREE SOBOLEV VECTOR FIELDS

Keywords: Sobolev vector fields, generalized flows, continuity equation, ODE, integral curves 35A02 -35D30 -35Q49 -34A12

We construct divergence-free Sobolev vector fields in C([0, 1]; W 1,r (T d ; R d )) with r < d and d ≥ 2 which simultaneously admit any finite number of distinct positive solutions to the continuity equation. We then show that the vector fields we produce have at least as many integral curves starting from L d -a.e. point of T d as the number of distinct positive solutions to the continuity equation these vector fields admit. Our work uses convex integration techniques introduced in [4, 20] to study nonuniqueness for positive solutions of the continuity equation. We then infer non-uniqueness for integral curves from Ambrosio's superposition principle.

Introduction

In this paper we study positive solutions of the continuity equation

∂ t ρ + div(ρu) = 0, ρ(•, t) = ρ 0 (•) (1.1)
where u : [0, 1] × T d → R d is a prescribed vector field on the d-dimensional torus and ρ 0 : T d → R is the initial datum. Throughout this work, (1.1) will be understood in the sense of distributions which only requires that ρ and ρu be integrable. We then study integral curves of the vector field u.

In the smooth setting, the Cauchy-Lipschitz theory guarantees the existence of a unique flow X : [0, 1] × T d → T d of the vector field u satisfying ∂ t X(t, x) = u(t, X(t, x)), X(0, x) = x.

(1.

2)

The classical Liouville theorem then gives a representation of solutions of (1.1) in terms of the flow X of the vector field u through the formula ρ(t, •)L d = X(t, .) # (ρ 0 L d ).

(1.3)

For rough vector fields, the relationship between the continuity equation and the corresponding flow is an active field of research since the foundational work of DiPerna and Lions in [START_REF] Diperna | Ordinary differential equations, transport theory and Sobolev spaces[END_REF]. By means of a regularization scheme, they showed that if u ∈ L 1 ((0, 1); W 1,r (T d )) and div u ∈ L 1 ((0, 1) × T d ), then (1.1) is well-posed in the class L ∞ ((0, 1); L p (T d )), where p, r ≥ 1 satisfy the relation

1 p + 1 r ≤ 1.
In [START_REF] Ambrosio | Transport equation and cauchy problem for bv vector fields[END_REF], Ambrosio extended the work of DiPerna and Lions to the setting of BV vector fields. We now gather some useful definitions to investigate the relation of the ODE (1.2) and the PDE (1.1) in the non-smooth setting. Definition 1.1. Let u : (0, 1) × T d → R d be a Borel map. We say that γ ∈ AC([0, 1]; T d ) is an integral curve of u starting at x if γ(0) = x and γ ′ (t) = u(t, γ(t)) for a.e. t ∈ (0, 1).

The regular Lagrangian flow is then a suitable selection of integral curves of u by a compressibility condition (introduced in [START_REF] Ambrosio | Transport equation and cauchy problem for bv vector fields[END_REF][START_REF] Diperna | Ordinary differential equations, transport theory and Sobolev spaces[END_REF]). Definition 1.2 (Regular Lagrangian flow). Let u : (0, 1) × T d → R d be Borel. We say that a Borel map X : [0, 1] × R d → R d is a regular Lagrangian flow of u if (i) for L d -a.e. x ∈ T d , t → X(t, x) is integral curve of u with X(0, x) = x, (ii) there is a constant C > 0 such that for every t ∈ [0, 1], X(t, .)

# L d ≤ CL d .
The well-posedness of the regular Lagrangian flow for vector fields u in L 1 ((0, 1); W 1,1 (T d )) with the negative part of the divergence satisfying [div u] -∈ L 1 ((0, 1); L ∞ (T d )) was first derived from the wellposedness of (1.1) for bounded densities, and for such densities the formula (1.2) holds using as X the regular Lagrangian flow (see [START_REF] Diperna | Ordinary differential equations, transport theory and Sobolev spaces[END_REF] and see [START_REF] Ambrosio | Transport equation and cauchy problem for bv vector fields[END_REF] for the BV vector fields case). Later in [START_REF] Crippa | Estimates and regularity results for the DiPerna-Lions flow[END_REF], Crippa and De Lellis proved well-posedness of the regular Lagrangian flow without resorting to the PDE (1.1), but their approach only works for vector fields in L 1 ((0, 1); W 1,r (T d )) with r > 1. At any rate, the uniqueness of the regular Lagrangian flow does not imply L d -a.e. uniqueness of integral curves. Indeed, Bruè, Colombo and De Lellis recently produced divergence-free Sobolev vector fields -uniqueness of the regular Lagrangian flow associated to these vector fields therefore holds -for which almost everywhere uniqueness of integral curves fails (see [START_REF] Brué | Positive Solutions of Transport Equations and Classical Nonuniqueness of Characteristic curves[END_REF]Theorem 1.3]). We note that the case of continuous vector field is still open (question posed in [1, Section 2.3]) although in [START_REF] Caravenna | A directional lipschitz extension lemma, with applications to uniqueness and lagrangianity for the continuity equation[END_REF], Crippa and Caravenna proved almost everywhere uniqueness of the trajectories when u ∈ C((0, 1); W 1,r ) when r > d.

In this work we show that, for divergence-free Sobolev vector fields, uniqueness of integral curves of the ODE (1.2) can fail for a set of initial data with full measure. In fact, we show that the non-uniqueness for integral curves is even worse: for any natural number N we produce divergence-free Sobolev vector fields with at least N integral curves starting almost everywhere. We highlight that our result demonstrates the power of the selection principle of the regular Lagrangian flow ((ii) in Definition 1.2) for integral curves of Sobolev vector fields. Indeed, amongst at least N integral curves starting from L d -a.e. point of T d , the regular Lagrangian flow selects a single integral curve for L d -a.e. starting point.

Theorem 1.3. For every d, N ∈ N, d ≥ 2, r ∈ [1, d) and s < ∞ there is a divergence-free vector field u ∈ C((0, 1); W 1,r (T d ; R d ) ∩ L s ) such that the following holds for every Borel map v with u = v L d+1 -a.e.:

(NU) For L d -a.e. x ∈ T d there are at least N integral curves of v starting at x.

Ambrosio's superposition principle [START_REF]Transport equation and Cauchy problem for non-smooth vector fields[END_REF]Theorem 3.2] bridges the gap between positive solutions of the continuity equation (1.1) for a vector field u and the integral curves of u (as in Definition 1.1): it gives a way of representing positive solutions of the continuity equation in terms of integral curves of the vector field without any differentiability assumption, i.e. under more general assumptions than DiPerna-Lions theory. Using Ambrosio's superposition principle, we will derive Theorem 1.3 from a non-uniqueness result for positive solutions of (1.1), which in turn will be proved using a convex integration iterative procedure. The term convex integration is generic to designate iterative techniques by which wild solutions of PDEs are constructed. Such techniques were introduced in the study of the continuity equation in the groundbreaking work of Modena and Székelyhidi [START_REF] Modena | Non-uniqueness for the transport equation with Sobolev vector fields[END_REF][START_REF]Non-renormalized solutions to the continuity equation[END_REF] (see also [5, 6, 8-11, 13-15, 17, 18] for interesting results using the convex integration methods).

Theorem 1.4. Let d, N ∈ N with d ≥ 2. Let p ∈ (1, ∞), r ∈ [1, ∞] be such that 1 p + 1 r > 1 + 1 d ,
and denote by p ′ the dual exponent of p, i.e. 1/p + 1/p ′ = 1. Then there exists a divergence-free vector field In order to prove Theorem 1.4, we adapt the convex integration scheme for positive solutions of the continuity equation introduced in [START_REF] Brué | Positive Solutions of Transport Equations and Classical Nonuniqueness of Characteristic curves[END_REF]: our proof makes use of two new ideas. Firstly, we keep track of a fixed number of densities and one single vector field throughout the iteration scheme. Each density is perturbed using a distinct family of building blocks 1 . Each of these building blocks then interacts only with one term of the perturbation to the vector field (see the key identity (4.11)). Secondly, we localize in space the corrector parts of the perturbation to the densities (see (3.15)) which will be negative, in order to preserve the positivity of the solutions. We also note that to prove Theorem 1.4 in dimension d = 2 the ideas of [START_REF] Brué | Positive Solutions of Transport Equations and Classical Nonuniqueness of Characteristic curves[END_REF]Section 7] need to be adapted for technical reasons. This will be explained in Section 6.

u ∈ C 0 ([0, 1], W 1,r (T d ) ∩ L p ′ (T d ))

Preliminary lemmas

In this section, we gather some useful lemmas from [START_REF] Brué | Positive Solutions of Transport Equations and Classical Nonuniqueness of Characteristic curves[END_REF][START_REF] Modena | Non-uniqueness for the transport equation with Sobolev vector fields[END_REF]. We will write

T d for R d /Z d . Lemma 2.1. Let d, N ∈ N + . Then, there exist disjoint families, Λ i , of finite sets {ξ} ξ∈Λ i ⊆ ∂B 1 ∩ Q d for i = 1, .., N and smooth nonnegative coefficients a ξ (R) such that for every R ∈ ∂B 1 R = ξ∈Λ i a ξ (R)ξ
for any i = 1, .., N .

2.1.

Antidivergences. We recall that the operator ∇∆ -1 is an antidivergence when applied to smooth vector fields of zero mean. The following lemma proven in [START_REF] Modena | Non-uniqueness for the transport equation with Sobolev vector fields[END_REF]Lemma 2.3] and [START_REF] Modena | Convex integration solutions to the transport equation with full dimensional concentration[END_REF]Lemma 3.5] gives an improved antidivergence operator for functions with a particular structure. Lemma 2.2. (Cp. with [START_REF] Modena | Convex integration solutions to the transport equation with full dimensional concentration[END_REF]Lemma 3.5]) Let λ ∈ N and f, g : T d → R be smooth functions, and

g λ = g(λx). Assume that ´g = 0. Then if we set R(f g λ ) = f ∇∆ -1 g λ -∇∆ -1 (∇f • ∇∆ -1 g λ + ´f g λ ),
we have that div R(f g λ ) = f g λ -´f g λ and for some C := C(k, p)

D k R(f g λ ) L p ≤ Cλ k-1 f C k+1 g W k,p for every k ∈ N, p ∈ [1, ∞]. ( 2
f (x)g(λx) L p ≤ f L p g L p + C(p) √ d f C 1 g L p λ 1/p (2.2) and ˆf (x)g(λx) dx ≤ ˆf (x) g(λx) - ˆg dx + ˆf • ˆg ≤ √ d f C 1 g L 1 λ + ˆf • ˆg . (2.
3)

1
The term building block refers to smooth functions or vector fields which are fixed before iteration, and which are used to construct perturbations in a convex integration iterative scheme.

Building blocks.

The building blocks are the same as those of [START_REF] Brué | Positive Solutions of Transport Equations and Classical Nonuniqueness of Characteristic curves[END_REF]Section 4]. We recall them here for the convenience of the reader. Let 0 < ρ < 1 4 be a constant. We consider

ϕ ∈ C ∞ c (B ρ ) and ψ ∈ C ∞ c (B 2ρ ) which satisfy ˆϕ = 1, ϕ ≥ 0, ψ ≡ 1 on B ρ .
Given µ ≪ 1 we define the 1-periodic functions φµ (x) :=

k∈Z d µ d/p ϕ(µ(x + k)) ψµ (x) := k∈Z d µ d/p ′ ψ(µ(x + k)) .
Let ω : R d → R be a smooth 1-periodic function such that ω(x) = x • ξ ′ on B 2ρ (0). Given Λ as in Lemma 2.1, for any ξ ∈ Λ we chose ξ ′ ∈ ∂B 1 such that ξ • ξ ′ = 0 and we define

Ω µ ξ (x) := µ -1 ω(µ x)(ξ ⊗ ξ ′ -ξ ′ ⊗ ξ). Notice that div Ω µ
ξ is divergence-free since Ω µ ξ is skew-symmetric and div Ω µ ξ = ξ on supp( ψµ ) and supp( φµ ). For σ > 0 we set Wξ,µ,σ (t, x)

:= σ 1/p ′ div (Ω µ ξ ψµ )(x -µ d/p ′ σ 1/p ′ tξ) (2.4) Θξ,µ,σ (t, x) := σ 1/p φµ (x -µ d/p ′ σ 1/p ′ tξ) . ( 2.5)
Notice that Wξ,µ,σ is divergence-free since it is also the divergence of the skew-symmetric matrix Ω µ ξ ψµ . By construction we have

Wξ,µ,σ (t, x) = σ 1/p ′ ψµ ξ + Ω µ ξ • ∇ ψµ (x -µ d/p ′ σ 1/p ′ tξ)
, hence the following properties are easily verified. If we consider the translations

W ξ,µ,σ (t, x) = Wξ,µ,σ (t, x -v ξ ), Θξ,µ,σ (t, x) = Θξ,µ,σ (t, x -v ξ ),
we can prove the following result.

Lemma 2.3. Let d ≥ 3, Λ ⊂ ∂B 1 ∩ Q d be a finite set. Then there exists µ 0 > 0 such that the following holds.

There exist two families of functions

{Θ ξ,µ,σ } ξ,µ,σ ⊂ C ∞ (T d ), {W ξ,µ,σ } ξ,µ,σ ⊂ C ∞ (T d ; R d )
, where ξ ∈ Λ, σ, µ ∈ R such that for any µ ≥ µ 0 , σ > 0 we have

∂ t Θ ξ,µ,σ + div(W ξ,µ,σ Θ ξ,µ,σ ) = 0, (2.6) div W ξ,µ,σ = 0, ˆWξ,µ,σ = 0, (2.7) ˆWξ,µ,σ Θ ξ,µ,σ = σξ. (2.8)
For any k ∈ N and any s ∈ [1, ∞] one has

D k Θ ξ,µ,σ L s ≤ C(d, k, s)σ 1/p µ k+d(1/p-1/s) , ∂ k t Θ ξ,µ,σ L s ≤ C(d, k, s)σ 1+ k-1 p ′ µ k+d( k-1 p ′ +1-1 s ) (2.9) D k W ξ,µ,σ L s ≤ C(d, k, s)σ 1/p ′ µ k+d(1/p ′ -1/s) , ∂ k t W ξ,µ,σ L s ≤ C(d, k, s)σ k+1 p ′ µ k+d( k+1 p ′ -1 s ) . (2.10)
Finally, they have pairwise compact disjoint supports for any ξ = ξ ′ , namely

suppW ξ,µ,σ ∩ suppΘ ξ ′ ,µ,σ = suppW ξ,µ,σ ∩ suppW ξ ′ ,µ,σ = suppΘ ξ ′ ,µ,σ ∩ suppΘ ξ,µ,σ = ∅, (2.11)
for any ξ = ξ ′ .

The proof of the previous lemma follows combining [START_REF] Brué | Positive Solutions of Transport Equations and Classical Nonuniqueness of Characteristic curves[END_REF]Lemma 4.1 and Lemma 4.2].

Iteration scheme

The convex integration scheme to construct solutions of the continuity equation was first introduced in [START_REF] Modena | Non-uniqueness for the transport equation with Sobolev vector fields[END_REF]. The scheme was later adapted in [START_REF] Brué | Positive Solutions of Transport Equations and Classical Nonuniqueness of Characteristic curves[END_REF] to construct positive solutions of the continuity equation. To prove Theorem 1.3 we will adapt the convex integration scheme of [START_REF] Brué | Positive Solutions of Transport Equations and Classical Nonuniqueness of Characteristic curves[END_REF].

First, we define the notion of a family of a-open sets. This notion is useful because of the convolution step in the iteration scheme.

Definition 3.1. Let N ∈ N, a ∈ R + and {A i } 1≤i≤N be a finite family of open sets of T d . We say that the family {A i } 1≤i≤N is a-open if for any i = 1, . . . , N there exists a ball B i of radius 1 4aN such that B i ⊂ A i .
As in [START_REF] Modena | Non-uniqueness for the transport equation with Sobolev vector fields[END_REF] we consider the following system of equations in

[0, 1] × T d , where d ≥ 2,      ∂ t ρ q,i + div(ρ q,i u q ) = -div R q,i div u q = 0, (3.1)
where the indices are i, q ∈ N and 1 ≤ i ≤ N . We then fix three parameters a 0 , b > 0 and β > 0, to be chosen later only in terms of d, p, r, and for any choice of a > a 0 we define

λ 0 = a, λ q+1 = λ b q and δ q = λ -2β q .
The following proposition builds a converging sequence of functions with the inductive estimates

max t R q,i (t, •) L 1 ≤ δ q+1 (3.2) max t ( ρ q,i (t, •) C 1 + ∂ t ρ q,i (t, •) C 0 + u q (t, •) W 1,p ′ + u q (t, •) W 2,r + ∂ t u q (t, •) L 1 ) ≤ λ α q , ( 3.3)
for any i = 1, . . . , N , where α is yet another positive parameter which will be specified later.

Proposition 3.2. Let d, N ∈ N, d ≥ 3.
There exist α, b, a 0 , M > 5, 0 < β < (2b) -1 such that the following holds. For every family {A i } 1≤i≤N of a 0 -open in T d and for every a ≥ a 0 , if {(ρ q,i , u q , R q,i )} 1≤i≤N solve (3.1) and enjoy the estimates (3.2), (3.3), then there exist {(ρ q+1,i , u q+1 , R q+1,i )} 1≤i≤N which solve (3.1) and enjoy the estimates (3.2), (3.3) with q replaced by q + 1. Moreover, for any i = 1, . . . , N , the following hold:

(a) sup [0,1] [ (ρ q+1,i -ρ q,i )(t, •) p L p + (u q+1 -u q )(t, •) r W 1,r + (u q+1 -u q )(t, •) p ′ L p ′ ] ≤ M δ q+1 ; (b) the following properties inf [0,1]×(T d \Ai) ρ q,i ≥ 0, inf [0,1]×Ai ρ q,i ≥ c > 0, imply inf [0,1]×(T d \Ai) ρ q+1,i ≥ 0, inf [0,1]×Ai ρ q+1,i ≥ c -δ q+1 , (c) if for some t 0 > 0 we have that ρ q,i (t, •) = 1, R q,i (t, •) = 0 and u q (t, •) = 0 for every t ∈ [0, t 0 ], then ρ q+1,i (t, •) = 1, R q+1,i (t, •) = 0 and u q+1 (t, •) = 0 for every t ∈ [0, t 0 -λ -1-α q ], (d) if for some t 0 > 0 we have that supp ρ q,i (t, •) ⊂ B i , R q,i (t, •) = 0 and u q (t, •) = 0 for every t ∈ [t 0 , 1], then supp ρ q+1,i (t, •) ⊂ B i,λ -1-α q , R q+1,i (t, •) = 0 and u q+1 (t, •) = 0 for every t ∈ [t 0 + λ -1-α q , 1], where B i,λ -1-α q := {x ∈ T d : d(x, B i ) < λ -1-α q }.
Remark 3.3. We highlight that the constant a 0 in the proposition above does not depend on the sets A i but only on the number N . Therefore, when we apply this proposition (to prove Theorem 1.4), we choose the sets A i after having fixed a 0 .

To prove Proposition 3.2 we use a convex integration scheme similar to the one in [4, Proposition 2.1]. However, the end products of our scheme are different from those of [START_REF] Brué | Positive Solutions of Transport Equations and Classical Nonuniqueness of Characteristic curves[END_REF]. Indeed, we seek to produce a single vector field u and N densities ρ i with mutually disjoint compact supports for some time such that (ρ i , u) weakly solves (1.1). Accordingly, we modified the iterative proposition of [START_REF] Brué | Positive Solutions of Transport Equations and Classical Nonuniqueness of Characteristic curves[END_REF] in two essential ways: we index N distinct densities ρ q,i by the parameter i = 1, . . . , N ; we have refined the control inf ρ q+1,i over the subregion T d \ A i thanks to (b). The former is achieved by taking N disjoint families of building blocks {Λ i } 1≤i≤N , and the latter by localizing to A i the corrector part of the perturbation to ρ q,i .

3.1. Choice of the parameters. The choice of parameter is the same of [4, Section 5.1]. We define first the constant

γ := 1 + 1 p min d p , d p ′ , -1 -d 1 p ′ - 1 r -1 > 0,
Notice that, up to enlarging r, we can assume that the quantity in the previous line is less than 1/2, namely that γ > 2. Hence we set α

:= 4 + γ(d + 1), b := max{p, p ′ }(3(1 + α)(d + 2) + 2), (3.4)
and

β := 1 2b min p, p ′ , r, 1 b + 1 = 1 2b(b + 1) . (3.5)
Finally, we choose a 0 and M sufficiently large (possibly depending on all previously fixed parameters) to absorb numerical constants in the inequalities. We set

ℓ := λ -1-α q , ( 3.6
)

µ q+1 := λ γ q+1 . ( 3.7)
3.2. Convolution. The convolution step is the same of [4, Section 5.2]. We just write here the definitions. We first perform a convolution of ρ q and u q to have estimates on more than one derivative of these objects and of the corresponding error. Let φ ∈ C ∞ c (B 1 ) be a standard convolution kernel in space-time, ℓ as in (3.6) and define ρ ℓ,i := ρ q,i * φ ℓ , u ℓ := u q * φ ℓ , R ℓ,i := R q,i * φ ℓ . We observe that (ρ ℓ,i , u ℓ , R ℓ,i + (ρ q,i u q ) ℓρ ℓ,i u ℓ ) solves system (3.1) for any i = 1, . . . , N and by (3.2), (3.5) enjoys the following estimates

R ℓ,i L 1 ≤ δ q+1 , (3.8) ρ ℓ,i -ρ q,i L p ≤ ℓ ρ q,i C 1 ≤ ℓλ α q ≤ λ -1 q ≤ δ 1/p q+1 2 , (3.9) u ℓ -u q L p ′ ≤ Cℓλ α q ≤ Cδ 1/p ′ q+1 , u ℓ -u q W 1,r ≤ Cℓλ α q ≤ Cδ 1/r
q+1 . Indeed note that by (3.5)

ℓλ α q = λ -1 q = δ 1 2bβ q+1 ≤ δ max{1/p,1/p ′ ,1/r} q+1 . Next observe that ∂ S t ρ ℓ,i C 0 + ρ ℓ,i C S + u ℓ W 1+S,r + ∂ S t u ℓ W 1,r ≤ C(S)ℓ -S+1 ( ρ q,i C 1 + u q W 2,r ) ≤ C(S)ℓ -S+1 λ α q
for every S ∈ N \ {0} and for every i = 1, . . . , N . Using the Sobolev embedding W d,r ⊂ W d,1 ⊂C 0 we then conclude

∂ S t u ℓ C 0 + u ℓ C S ≤ C(S)ℓ -S-d+2 λ α q
. By Young's inequality we estimate the higher derivatives of R ℓ,i in terms of R q,i L 1 to get

R ℓ,i C S + ∂ S t R ℓ,i C 0 ≤ D S ρ ℓ,i L ∞ R q,i L 1 ≤ C(S)ℓ -S-d ≤ C(S)λ (1+α)(d+S) q (3.10)
for every L ∈ N and i = 1, . . . , N . Finally, thanks to [4, Lemma 5.1] for the last part of the error we have

(ρ q,i u q ) ℓ -ρ ℓ,i u ℓ L 1 ≤ Cℓ 2 λ 2α q ≤ 1 4 δ q+2 , (3.11)
where we have assumed that a is sufficiently large.

3.3. Definition of the perturbation. Let µ q+1 > 0 be as in (3.7) and let

χ ∈ C ∞ c (-3 4 , 3 4 ) such that n∈Z χ(τ -n) = 1 for every τ ∈ R. Let χ ∈ C ∞ c (-4 5 , 4 5 ) be a nonnegative function satisfying χ = 1 on [-3 4 , 3 4 ]. Notice that n∈Z χ(τ -n) ∈ [1, 2] and χ • χ = χ. Fix a parameter κ = 40p ′ /δ q+2 and consider 2N disjoint families {Λ 1 i } 1≤i≤N , {Λ 2 i } 1≤i≤N as in Lemma 2.1. Next, for n ∈ N, define [n]
to be 1 or 2 depending on the congruence class of n. Finally, we take our building blocks according to Lemma 2.3 with Λ = ∪ 2 i=1 ∪ N j=1 Λ i j and observe that their spatial supports are disjoint. We define the new density and vector field by adding to ρ ℓ and u ℓ principal terms and correctors, namely we set

ρ q+1,i := ρ ℓ,i + θ (p) q+1,i + θ (c) q+1,i , u q+1 := u ℓ + N i=1 (w (p) q+1,i + w (c) q+1,i ) .
The principal perturbations are given, respectively, by

w (p) q+1,i (t, x) = n≥12 χ(κ|R ℓ,i (t, x)| -n) ξ∈Λ [n] i W ξ,µq+1,n/κ (λ q+1 t, λ q+1 x), (3.12) θ (p) q+1,i (t, x) = n≥12 χ(κ|R ℓ,i (t, x)| -n) ξ∈Λ [n] i a ξ R ℓ,i (t, x) |R ℓ,i (t, x)| Θ ξ,µq+1,n/κ (λ q+1 t, λ q+1 x) , (3.13)
where it is understood that the terms in the second sum vanish at points where R ℓ,i vanishes. In the definition of w

(p)
q+1,i and θ

(p)
q+1,i the first sum runs for n in the range

12 ≤ n ≤ Cℓ -d δ -1 q+2 ≤ Cλ d(1+α)+2βb 2 q ≤ Cλ d(1+α)+1 q ≤ λ d(1+α)+2 q , ( 3.14)
where the last holds providing

a 0 ≥ C. Indeed χ(κ|R ℓ,i (t, x)| -n) = 0 if n ≥ 20δ -1 q+2 R ℓ,i C 0 + 1

and by (3.10) we obtain an upper bound for n.

The aim of the corrector term for the density is to ensure that the overall perturbation has zero average. So we set

θ (c) q+1,i (t, x) := -g i (x) ˆθ(p) q+1,i (t, x) dx, ( 3.15)
where

{g i } 1≤i≤N ⊂ C ∞ (T d ) such that ´Td g i = 1
, g i ≥ 0 and supp(g i ) is compactly contained on A i and g i C S (λ 0 N ) d+S (where means inequality up to a geometric constant depending only on d). Here we used the property that the family {A i } 1≤i≤N is λ 0 -open. We observe that the functions {g i } 1≤i≤N do not depend on q, they only depend only on the fixed open sets {A i } 1≤i≤N . The aim of the corrector term for the vector field is to ensure that the overall perturbation has zero divergence. Thanks to (2.7), we can apply Lemma 2.2 to define

w (c) q+1,i := - n≥12 ξ∈Λ [n] i R ∇χ(κ|R ℓ,i (t, x)| -n) • W ξ,µq+1,n/κ (λ q+1 t, λ q+1 x)
Moreover, since W ξ,µq+1,n/κ is divergence-free, the argument inside R has 0 average for every t ≥ 0, and so w (c) q+1,i is indeed well defined. Notice finally that the perturbations equals 0 on every time interval where R ℓ,i vanishes identically for any i = 1, . . . , N .

Proof of Proposition 3.2

For the sake of readability, the quantifier "for every i = 1, . . . , N " will be implicit in the rest of this paper. Before coming to the main arguments, we recall [4, Lemma 6.1] for the "slowly varying coefficients".

∂ m t χ(κ|R ℓ,i | -n) C S + ∂ m t χ(κ|R ℓ,i | -n) C S ≤ C(m, S)δ -2(S+m) q+2 ℓ -(S+m)(1+d) ≤ C(m, S)λ (S+m)(d+2)(1+α) q ∂ m t (a ξ ( R ℓ,i |R ℓ,i | )) C S ≤ C(m, S)δ -S-m q+2 ℓ -(S+m)(1+d) ≤ C(m, S)λ (S+m)(d+2)(1+α) q on {χ(κ|R ℓ,i | -n) > 0}.
4

q+1,i L p ≤ n≥12 ξ∈Λ [n] i χ(κ|R ℓ,i (t, x)| -n)a ξ R ℓ,i (t,x) |R ℓ,i (t,x)| L p Θ ξ,µq+1,n/κ (λ q+1 t, λ q+1 x) L p + 1 λ 1/p q+1 n≥12 ξ∈Λ [n] χ(κ|R ℓ,i (t, x)| -n)a ξ R ℓ,i (t,x) |R ℓ,i (t,x)| C 1 Θ ξ,µq+1,n/κ (λ q+1 t, λ q+1 x) L p ≤ C n≥12 (n/k) 1/p χ(κ|R ℓ,i (t, x)| -n) L p + Cλ -1/p q+1 δ 1/p q+2 λ (d+2)(1+α)+(1+1/p)((d(1+α)+1)) q ≤ C R ℓ,i 1/p L 1 + Cλ -1/p q+1 δ 1/p q+2 λ 3(d+2)(1+α) q ≤ Cδ 1/p q+1 , (4.1)
provided that in the second last inequality we use

(d + 2)(1 + α) + (1 + 1/p)((d(1 + α) + 1)) ≤ 3(1 + α)(d + 2)
and in the last inequality we use (3.8).

Next , we use Θ ξ,µq+1,n/κ L 1 ≤ ( n κ )

q+1,i (t, x)| ≤ √ dλ -1 q+1 n≥12 ξ∈Λ [n] g i L ∞ χ(κ|R ℓ,i | -n)a ξ R ℓ,i |R ℓ,i | C 1 Θ ξ,µq+1,n/κ L 1 + n≥12 ξ∈Λ [n] g i L ∞ χ(κ|R ℓ,i | -n)a ξ R ℓ,i |R ℓ,i | L 1 Θ ξ,µq+1,n/κ L 1 ≤ g i L ∞ λ -1 q+1 µ -d/p ′ q+1 λ 3(1+α)(d+2) q + g i L ∞ µ -d/p ′ q+1 n≥12 (n/κ) 1/p χ(κ|R ℓ,i | -n)a ξ R ℓ,i |R ℓ,i | L 1 ≤ g i L ∞ µ -d/p ′ q+1 + g i L ∞ µ -d/p ′ q+1 R ℓ,i 1/p L 1 (N λ 0 ) d λ -1 q+1 ≤ δ q+1 2 , ( 4.2)
where in the second to last inequality, we enlarge a to absorb the constant N . Now, recall that θ

(p) q+1 is nonnegative by definition. Therefore, inf [0,1]×T d [θ (p) q+1,i + θ (c) q+1,i ] ≥ inf [0,1]×T d θ (c) q+1,i ≥ - δ q+1 2 . Also inf [0,1]×(T d \Ai) [θ (p) q+1,i + θ (c) q+1,i ] ≥ 0.
Since ρ ℓ,i is nonnegative whenever ρ q,i is nonnegative, by (3.9) we get property (b) of Proposition 3.2.

4.2. Estimate on w q+1,i L p ′ and Dw q+1,i L r . Exactly with the same computation as in (4.1), replacing p with p ′ , we have that

w (p) q+1,i L p ′ ≤ C R ℓ,i 1/p ′ L 1 + Cλ -1/p ′ q+1 δ 1/p ′ q+2 λ 3(d+2)(1+α) q ≤ Cδ 1/p ′ q+1 (4.3)
Concerning the corrector term w

(c) q+1 , we use (2.10) (precisely W ξ,µq+1,n/κ L p ′ ≤ ( n κ ) 1/p ′ ≤ λ (d+2)(1+α) q ) Lemma 2.2 and (3.14) to get w (c) q+1,i L p ′ ≤ 1 λ q+1 n≥12 ξ∈Λ [n] i χ(κ|R ℓ,i | -n) C 2 W ξ,µq+1,n/κ L p ′ ≤ Cλ -1 q+1 n≥12 λ 2(d+2)(1+α) q (n/κ) 1/p ′ ≤ Cλ -1 q+1 δ 1/p ′ q+2 λ 4(d+2)(1+α) q ≤ δ 1/p ′ q+2 ≤ δ 1/p ′ q+1 . (4.4)
Computing the gradient of w (p) q+1 and combining Lemma 4.1 with (2.10) we have

Dw (p) q+1,i L r ≤ n≥12 ξ∈Λ [n] i χ(κ|R ℓ,i | -n) C 1 W ξ,µq+1,n/κ L r + n≥12 ξ∈Λ [n] i λ q+1 DW ξ,µq+1,n/κ L r (4.5) ≤ Cδ 1/p ′ q+2 λ 3(1+α)(d+2)+2+bγd(1/p ′ -1/r) q + Cδ 1/p ′ q+2 λ b+3(1+α)(d+2)+bγ(1+d(1/p ′ -1/r)) q (4.6) ≤ δ 1/p ′ q+2 ≤ δ 1/r q+1 . ( 4.7)
Concerning the corrector, by Lemma 2.2 and similar computations as above,

Dw (c) q+1,i L r ≤ C n≥12 ξ∈Λ [n] i χ(κ|R ℓ,i | -n) C 3 W ξ,µq+1,n/κ L r (4.8) ≤ Cδ 1/p ′ q+2 λ 5(1+α)(d+2) q µ d( 1 r -1 p ′ ) q+1 ≤ δ 1/p ′ q+2 λ 5(1+α)(d+2)-b(1+1/p) q ≤ δ 1/r q+1 .
(4.9) 4.3. Definition of the new error R q+1,i . This part is the similar to [START_REF] Brué | Positive Solutions of Transport Equations and Classical Nonuniqueness of Characteristic curves[END_REF], however here θ (c) q+1,i is not constant in space and we have to adapt the argument accordingly: we will pick up a new error term which we call R space .

By definition the new error R q+1,i must satisfy

-div R q+1,i = ∂ t ρ q+1,i + div(ρ q+1,i u q+1 ) = div(θ (p) q+1,i w (p) q+1 -R ℓ,i ) + ∂ t θ (p) q+1,i + ∂ t θ (c) q+1,i + div(θ (p) q+1,i u ℓ + ρ ℓ,i w q+1 + θ (p) q+1,i w (c) q+1 ) + div((u ℓ + w q+1 )θ (c) q+1,i ) + div((ρ q,i u q ) ℓ -ρ ℓ,i u ℓ ) (4.10)
In the second equality above we have used that (ρ ℓ , u ℓ , R ℓ,i + (ρ q,i u q ) ℓρ ℓ,i u ℓ ) solves (3.1).

We now decompose

∂ t θ (p) q+1,i = n≥12 ξ∈Λ [n] i χ(κ|R ℓ,i | -n)a ξ R ℓ,i |R ℓ,i | ∂ t Θ ξ,µq+1,n/κ (λ q+1 t, λ q+1 x) + n≥12 ξ∈Λ [n] i ∂ t χ(κ|R ℓ,i | -n)a ξ R ℓ,i |R ℓ,i | Θ ξ,µq+1,n/κ (λ q+1 t, λ q+1 x) =: (∂ t θ (p) q+1,i ) 1 + (∂ t θ (p) q+1,i ) 2 ,
We now observe that

θ (p) q+1,i N k=1 w (p) q+1,k = θ (p) q+1,i w (p) q+1,i = n≥12 ξ∈Λ [n] i χ(κ|R ℓ,i | -n)a ξ R ℓ,i |R ℓ,i | (Θ ξ,µq+1,n/κ W ξ,µq+1,n/κ )(λ q+1 t, λ q+1 x), (4.11)
where the first equality holds because Λ k i ∩ Λ p j = ∅ for any p, k and i = j; the second equality follows from (2.11) and the definitions of χ and χ. Also, Θ ξ,µq+1,n/κ and W ξ,µq+1,n/κ solve the transport equation (2.6). These observations in conjunction with Lemma 2.1 yield the cancellation of the error R ℓ,i up to lower order terms div(θ

(p) q+1,i N k=1 w (p) q+1,k ) + (∂ t θ (p) q+1,i ) 1 -div R ℓ,i = n≥12 ξ∈Λ [n] i ∇ χ(κ|R ℓ,i | -n)a ξ R ℓ,i |R ℓ,i | (Θ ξ,µq+1,n/κ W ξ,µq+1,n/κ )(λ q+1 t, λ q+1 x) -div R ℓ,i + n≥12 ξ∈Λ [n] i χ(κ|R ℓ,i | -n)a ξ R ℓ,i |R ℓ,i | λ q+1 ∂ t Θ ξ,µq+1,n/κ + div(Θ ξ,µq+1,n/κ W ξ,µq+1,n/κ ) (λ q+1 t, λ q+1 x) = n≥12 ξ∈Λ [n] i ∇ χ(κ|R ℓ,i | -n)a ξ R ℓ,i |R ℓ,i | (Θ ξ,µq+1,n/κ W ξ,µq+1,n/κ )(λ q+1 t, λ q+1 x) - n κ ξ + n≥12 ξ∈Λ [n] i ∇ χ(κ|R ℓ,i | -n)a ξ R ℓ,i |R ℓ,i | n κ ξ -div R ℓ,i = n≥12 ξ∈Λ [n] i ∇ χ(κ|R ℓ,i | -n)a ξ R ℓ,i |R ℓ,i | (Θ ξ,µq+1,n/κ W ξ,µq+1,n/κ )(λ q+1 t, λ q+1 x) - n κ ξ + div( Rℓ,i -R ℓ,i ), (4.12)
where Rℓ,i :

= n≥12 χ(κ|R ℓ,i | -n) R ℓ,i |R ℓ,i | n k .
We have We can now define R q+1,i which satisfies (4.10) as

|R ℓ,i -Rℓ,i | ≤ 11 n=-1 χ(κ|R ℓ,i | -n)R ℓ,i + n≥12 χ(κ|R ℓ,i | -n) R ℓ,i |R ℓ,i | n k -R ℓ,i ≤ 13 κ + n≥12 χ(κ|R ℓ,i | -n) |R ℓ,i | - n κ ( 4
-R q+1,i :=R quadr i + ( Rℓ,i -R ℓ,i ) + R time i + R space i + θ (p) q+1,i u ℓ + ρ ℓ,i w q+1 + θ (p) q+1,i w (c) q+1 + [(ρ q,i u q ) ℓ -ρ ℓ,i u ℓ ], (4.14)
where

R quadr i := n≥12 ξ∈Λ [n] i R ∇ χ(κ|R ℓ,i | -n)a ξ R ℓ,i |R ℓ,i | • (Θ ξ,µq+1,n/κ W ξ,µq+1,n/κ )(λ q+1 t, λ q+1 x) - n κ ξ , (4.15)
R time i := ∇∆ -1 ((∂ t θ (p) q+1,i ) 2 + ∂ t θ (c) q+1,i + m i ), (4.16
)

m i := n≥12 ξ∈Λ [n] i ˆ∇ χ(κ|R ℓ,i | -n)a ξ R ℓ,i |R ℓ,i | (Θ ξ,µq+1,n/κ W ξ,µq+1,n/κ )(λ q+1 t, λ q+1 x) - n κ ξ dx, R space i := (u ℓ + w q+1 )θ (c) q+1,i . ( 4
q+1,i ) 2 + ∂ t θ (c) q+1,i + m i has 0 mean, so that R time i is well defined. 4.4. Estimate on R q+1,i L 1 .
Recall that the estimate on (ρ q,i u q ) ℓρ ℓ,i u ℓ L 1 has been already established in (3.11). By the property (2.1) of the antidivergence operator R, Lemma 4.1 and (3.14) we have

R quadr i L 1 ≤ C λ q+1 n≥12 ξ∈Λ [n] i χ(κ|R ℓ,i | -n)a ξ R ℓ,i |R ℓ,i | C 2 Θ ξ,µq+1,n/κ W ξ,µq+1,n/κ L 1 ≤ Cδ q+2 λ 4(1+α)(d+2)+2 q λ q+1 ≤ δ q+2 20 .
To estimate the terms which are linear with respect to the fast variables, we take advantage of the concentration parameter µ q+1 . First of all, by Calderon-Zygmund estimates we get

R time i L 1 ≤ C (∂ t θ (p) q+1,i ) 2 + ∂ t θ (c) q+1,i -m i L 1 ≤ (∂ t θ (p) q+1,i ) 2 L 1 + |∂ t θ (c) q+1,i | + |m i |. Next, notice that (∂ t θ (p) q+1,i ) 2 L 1 ≤ C n≥12 ξ∈Λ [n] ∂ t χ(κ|R ℓ,i | -n)a ξ R ℓ,i |R ℓ,i | C 0 Θ ξ,µq+1,n/κ L 1 (4.18) ≤ Cδ 1/p q+2 λ 3(1+α)(d+2) q µ -d/p ′ q+1 ≤ δ q+2 20λ q+1 . ( 4
|∂ t θ (c) q+1,i | + |m i | ≤ g i L ∞ n≥12 ξ∈Λ [n] ˆχ(κ|R ℓ,i | -n)a ξ R ℓ,i |R ℓ,i | ∂ t Θ ξ,µq+1,n/κ (λ q+1 t, λ q+1 x) dx + |m i | + g i L ∞ (∂ t θ (p) q+1,i ) 2 L 1 ≤ g i L ∞
n≥12 ξ∈Λ [n] ˆχ(κ|R

ℓ,i | -n)a ξ R ℓ,i |R ℓ,i | div (Θ ξ,µq+1,n/κ W ξ,µq+1,n/κ )(λ q+1 t, λ q+1 x) dx + |m i | + δ q+2 20 = 2 g i L ∞ n≥12 ξ∈Λ [n] ˆ∇ χ(κ|R ℓ,i | -n)a ξ R ℓ,i |R ℓ,i | • (Θ ξ,µq+1,n/κ W ξ,µq+1,n/κ )(λ q+1 t, λ q+1 x) - n k ξ dx + δ q+2 20 ≤ δ q+2 20 + 2 g i L ∞ √ d λ q+1 n≥12 ξ∈Λ [n] χ(κ|R ℓ,i | -n)a ξ R ℓ,i |R ℓ,i | C 2 Θ ξ,µq+1,n/κ W ξ,µq+1,n/κ L 1 ≤ δ q+2 20 + C g i L ∞ λ -1 q+1 δ q+2 λ 4(1+α)(d+2) q ≤ 1 10 δ q+2 ,
and using also (4.2) and (3.3) , we estimate the error

R space i R space i L 1 ≤ (u ℓ + w q+1 )θ (c) q+1,i L 1 ≤ ( u ℓ L 1 + w q+1 L 1 ) θ (c) q+1,i L ∞ ≤ 2λ α q (N λ 0 ) d λ -1 q+1 ≤ δ q+2 20 ,
where the last inequality holds up to enlarging a 0 depending on N . We also have that

θ (p) q+1,i u ℓ + ρ ℓ,i w (p) q+1 L 1 ≤ θ (p) q+1,i L 1 u ℓ L ∞ + ρ ℓ,i L ∞ w (p) q+1 L 1 ≤ n≥12 ξ∈Λ [n] i χ(κ|R ℓ,i | -n)a ξ R ℓ,i |R ℓ,i | L ∞ Θ ξ,µq+1,n/κ L 1 u ℓ L ∞ + ρ ℓ,i L ∞ χ(κ|R ℓ,i | -n) L ∞ W ξ,µq+1,n/κ L 1 ≤ Cδ 1/p q+2 λ 2(1+α)(d+2) q µ -d/p ′ q+1 + Cδ 1/p ′ q+2 λ 2(1+α)(d+2) q µ -d/p q+1 ≤ δ q+2 20 .
In the last inequality we used 2βb 2 ≤ 1, the definition of γ, and b(1

+ 1/p) ≥ 2(1 + α)(d + 2) + 1.
Finally, from (4.1) and (4.4)

(ρ ℓ,i + θ (p) q+1,i )w (c) q+1 L 1 ≤ ( ρ ℓ,i C 1 + θ (p) q+1,i L p ) w (c) q+1 L p ′ ≤ Cλ 4(1+α)(d+2)+α q λ -1 q+1 ≤ 1 20 δ q+2 . (4.20)
4.5. Estimates on higher derivatives.

ρ q+1,i C 1 ≤ ρ ℓ,i C 1 + θ q+1,i C 1 ≤ ρ q,i C 1 + n≥10 ξ∈Λ [n] i χ(κ|R ℓ,i | -n)a ξ R ℓ,i |R ℓ,i | C 1 Θ ξ,µq+1,n/κ (λ q+1 x) C 1 ≤ Cλ α q + Cλ 3(1+α)(d+2) q λ q+1 µ 1+d/p q+1 ≤ λ α q+1 . (4.21)
An entirely similar estimate is valid for ∂ t ρ q+1,i C 0 and the one for

u ℓ + w (p) q+1,i W 2,r is analogous. Concerning w (c)
q+1,i , we use Lemma 2.2 and (3.14)

w (c) q+1,i W 2,r ≤ n≥12 ξ∈Λ [n] i λ q+1 χ(κ|R ℓ,i | -n) C 4 W ξ,µq+1,n/κ W 2,r ≤ Cλ 6(1+α)(d+2) q λ 2 q+1 µ 2+d(1/p ′ -1/r) q+1 ≤ λ α q+1 . It remains just to estimate ∂ t u q+1 L 1 ≤ ∂ t u ℓ L 1 + N i=1 ∂ t w (p) q+1,i L 1 + ∂ t w (c) q+1,i L 1 .
From (2.10) and Lemma 4.1

∂ t w (p) q+1,i L 1 ≤ n≥12 ξ∈Λ [n] i λ q+1 ∂ t W ξ,µq+1,n/κ L 1 + ∂ t χ(κ|R ℓ,i | -n) L ∞ W ξ,µq+1,κ/n L 1 Cδ 2/p ′ q+2 λ (1+2/p ′ )(d(1+α)+1) q λ q+1 µ 1+γ(1+d(2/p ′ -1)) q+1 ≤ λ 2+γ(d+1) q+1 ≤ λ α q+1 .
A similar computation is valid for

∂ t w (c) q+1,i L 1 .
5. Proof of main results

Proof of Theorem 1.4.

Proof of Theorem 1.4 assuming Proposition 3.2. Let α, b, a 0 , M > 5, β > 0 be fixed as in Proposition 3.2.

Let a ≥ a 0 be chosen such that

+∞ q=0 δ 1/p q+1 < 1 32M , +∞ q=0 λ -1-α q < 1 32N
,

Let {φ i } 1≤i≤N ⊂ C ∞ (T d ) be nonnegative functions with mutually disjoint compact supports such that ´Td φ i (x)dx = 1, {x ∈ T d : φ i (x) > 1} contains a ball of radius 1
4N , and d(supp φ i , supp φ j ) ≥ 1/4N for i = j. We also require that φ i C S ≤ (100N ) d+S for any S ∈ N. Let χ : [0, 1] → [0, 1] be a smooth function such that χ ≡ 0 on [0, 2/5], χ ≡ 1 on [3/5, 1] with χ ′ is compactly supported on (2/5, 3/5), and

∂ t χ L ∞ ≤ 20. Define ρ 0,i (t, x) := (1 -χ(t)) + χ(t)φ i (λ 0 x) and set u 0 ≡ 0. We also set R 0,i (t) := -∇∆ -1 ∂ t ρ 0,i (t) + div(ρ 0,i (t)u 0 (t)) = -∇∆ -1 ∂ t ρ 0,i (t) = -∂ t χ∇∆ -1 (φ i (λ 0 •) -1) .
We then have N starting triples {(ρ 0,i , u 0 , R 0,i )} 1≤i≤N for our iteration scheme which enjoy (3.1) with q = 0 for any i = 1, . . . , N . Moreover, thanks to Lemma 2.2, we have R 0,i L 1 ≤ Cλ -1 0 . Thus (3.2) is satisfied because 2β < 1 (here we have taken λ 0 = a 0 sufficiently large to absorb the constant C). Next, we have ∂ t ρ 0,i C 0 + ρ 0,i C 1 ≤ Cλ 0 . Since u 0 ≡ 0 and α > 1 we conclude that (3.3) is satisfied as well.

Finally we observe that the family of sets

A i := {x ∈ T d : φ(λ 0 x) > 1} for i = 1, . . . , N form a a 0 -open family.
We can recursively apply Proposition 3.2 to obtain a family of sequences {(ρ q,i , u q , R q,i ) q∈N } 1≤i≤N of smooth solutions to (3.1) and such that

• the sequences {ρ q,i } q∈N is Cauchy in C(L p ) and we denote by ρ i its limit for any i = 1, . . . , N ,

• the sequence of divergence-free {u q } q∈N is Cauchy in C(L p ′ ∩ W 1,r ) and we denote by u its limit (whose divergence understood in the sense of distribution vanishes). Thanks to property (3.2) we get that (u, ρ i ) solve the continuity equation for any i = 1, . . . , N . Property (b) and inf Ai ρ 0,i (t,

•) ≥ 1 also yield inf Ai ρ i (t, •) ≥ 1 - +∞ q=0 δ 1/p q+1 ≥ 1 2 .
This implies that A i ⊂ supp(ρ i (t, •)) for any t ∈ [0, 1] and any i = 1, • • • , N . Thus supp(ρ i (t, •)) has non-empty interior, and inf

T d \Ai ρ i ≥ 0.
So ρ i are nonnegative.

Finally, since ρ 0,i (t, •) ≡ 1 for t ∈ [0, 2/5] and

+∞ q=0 λ -1-α q < 1 λ 1+α 0 < 1 15
, and by property (c) and (d) of Proposition 3.2, we get that ρ i (t, •) ≡ 1 for t ∈ [0, 1/3] for any i = 1, . . . , N . Also, by property (d), and since

+∞ q=0 λ -1-α q < 1
λ0N and d(supp ρ 0,i (1, •), supp ρ 0,j (1, •)) ≥ 1 λ0N for i = j, we must have that supp(ρ i ) ∩ supp(ρ j ) is negligible for i = j.

Proof of Theorem 1.3.

Proof of Theorem 1.3 assuming Theorem 1.4. Let {ρ 0 } 1≤i≤N ⊂ C t L p be nonnegative densities and u ∈ C t (L p ′ ∩ W 1,r ) a divergence-free vector field given by Theorem 1.4. Then (ρ i , u) solves (1.1) for i = 1, . . . , N . Thanks to the Ambrosio's superposition principle (see [START_REF]Transport equation and Cauchy problem for non-smooth vector fields[END_REF]Theorem 3.2]), each nonnegative L 1 ([0, 1] × T d ) solution is transported by a generalized flow η i of the vector field u. More precisely,

η i ∈ M + (AC([0, 1]; T d ) × T d
) is concentrated on pairs (γ, x) such that γ an integral curve of u starting from x, and we have ρ i (x, t)L d = (e t ) ♯ η i for every t ∈ [0, 1].

Observe that the family of probability measures {η i } 1≤i≤N does not depend on the pointwise representative of u. Indeed, given two pointwise representative v and w of u (u and v are two Borel maps such that v = u = w L d+1 -a.e.), by Fubini and by the superposition principle, we have for each integer 1

≤ i ≤ N ˆAC([0,1],T d )×T d ˆ1 0 |v(γ(s)) -w(γ(s))|ds dη i (γ, x) = ˆ1 0 ˆTd |v(y) -w(y)|ρ i (s, x)dL d (x) ds = 0.
(5.1)

Thus, η i is concentrated on integral curves of v if and only if η i is concentrated on integral curves of w.

By the superposition principle we have

ˆTd ψ(x)ρ i (1, x)dL d (x) = ˆTd ˆAC([0,1],T d ) ψ(γ(1))dη i x (γ)dL d (x), (5.2)
and for every ψ ∈ C(T d ). Therefore, for L d -a.e. x ∈ T d and η i x -a.e. γ ∈ AC([0, 1]; T d ), we have γ(1) ∈ A i . Since A i ∩ A j = ∅ for i = j, it follows that for L d -a.e. x ∈ T d , the measures {η i x } 1≤i≤N have mutually disjoint supports. Therefore, for L d -a.e. x ∈ T d , there are at least N integral curves starting from x.

Dimension d = 2

The two dimensional case (i.e. for d = 2) is slightly more technical. We can no longer use Lemma 4.2 of [START_REF] Brué | Positive Solutions of Transport Equations and Classical Nonuniqueness of Characteristic curves[END_REF] to translate in space the tubes supporting the building blocks and thereby make these tubes disjoint. In [START_REF] Brué | Positive Solutions of Transport Equations and Classical Nonuniqueness of Characteristic curves[END_REF], the authors found a way around this issue. They are able make the building blocks in the case d = 2 disjoint. They take advantage of the presence of a single error to argue that only building blocks with comparable speeds -that is building blocks for which the speed ratio is of order ∼ 10 -2 need to have disjoint supports. Indeed in [START_REF] Brué | Positive Solutions of Transport Equations and Classical Nonuniqueness of Characteristic curves[END_REF], the speeds at the inductive step q ∈ N of the convex integration scheme are

w n = µ d/p ′ q+1 ( n κ ) 1/p ′ for n = 1, ..., λ d(1+α)+2 q
. The supports of the building blocks are then translated in space suitably, the speeds {w n } n=1,..,λ d(1+α)+2 q are approximated by {v n } n=1,..,λ d(1+α)+2 q , and at the price of a small error the authors obtain building blocks satisfying

W ξ,µq+1,vn • Θ ξ ′ ,µq+1,vm (t, x) = 0 for any (x, t) ∈ T 2 × R + , ( 6.1)
for any ξ = ξ ′ , |n -m| ≤ 1, n, m = 1, .., λ d(1+α)+2 q . However, they deal with only two distinct families Λ 1 and Λ 2 of directions for the building blocks in and a single error R q at each step q of the iteration, whereas in our setting there are 2N distinct families of directions {Λ i } 1≤i≤2N for the building blocks and N errors {R q,i } 1≤i≤N . We therefore need that any building block with direction in Λ i has disjoint support with any other building block with direction in Λ j for any i = j because in the convex integration scheme we need the key identity (4.11) to hold. More precisely we will need W ξ,µq+1,vn • Θ ξ ′ ,µq+1,vm (t, x) = 0 for any (x, t) ∈ T 2 × R + , (6.2)

whenever ξ = ξ ′ ∈ Λ, n, m = 1, .., λ d(1+α)+2 q
. This identity is achievable because the speed ratios of the building blocks are at most of order λ d(1+α)+2 q , typically a very small number compared to µ q+1 in the iterative proposition (see Section 3.1). So we will prove that we can find ∼ λ d(1+α)+2 q balls of radius ∼ µ -1 q+1 which are moving with speed ratio at most ∼ λ d(1+α)+2 q and which don't intersect at any time. We will proceed similarly to [4, Section 7], although our argument differs in some parts for reasons which were outlined above.

Lemma 6.1. Let ξ, ξ ′ ∈ S 1 ∩ Q 2 be two distinct vectors and let w = A N < λ d(1+α)+2 q
where A and N are positive, coprime integers such that N < λ d(1+α)+2 q

. Then there exists C = C(ξ, ξ ′ ) such that for any ε > 0

L 1 ([0, 1] \ {s : d T 2 (tξ, (tw + s)ξ ′ ) ≥ ε ∀ t ≥ 0}) < CN ελ d(1+α)+2 q . Proof. Let ε > 0. Set T int := {(t, t ′ ) : ξt = ξ ′ t ′ on T 2 }
and observe that T int ⊂ Q 2 since the matrix with columns ξ and ξ ′ is invertible with rational coefficients. Moreover T int is an additive discrete subgroup of R 2 , hence it is a free group of rank k ∈ {0, 1, 2}. Denoting by T and T ′ the period of, respectively, t → ξt and t → ξ ′ t one has that (T, 0), (0, T ′ ) ∈ T int . This implies that the rank of T int is two, hence we can find two generators (t 1 , t ′ 1 ), (t 2 , t ′ 2 ) ∈ T int . Let us finally introduce A := {ξt ∈ T 2 : (t, s) ∈ T int for some s ∈ R} to denote the set of points in T 2 where the supports of the curves t → tξ and t → tξ ′ intersect. Let s ∈ [0, 1] be such that d T 2 (tξ, (t + s)wξ ′ ) < ε for some t ≥ 0. There exists q ∈ A such that d T 2 (tξ, q) ≤ cε, where c = c(ξ, ξ ′ ) > 1, hence up to modifying t we can assume that tξ =: q ∈ A and d T 2 (q, (tw + s)ξ ′ ) ≤ 3cελ d(1+α)+2 q . Since tξ ∈ A there exists t ′ such that (t, t ′ ) ∈ T int and, exploiting the fact that (t 1 , t ′ 1 ), (t 2 , t ′ 2 ) ∈ T int are generators, we can find

k 1 , k 2 ∈ Z such that t = k 1 t 1 + k 2 t 2 and t ′ = k 1 t ′ 1 + k 2 t ′ 2 .
The following identity holds on T

2 (tw + s)ξ ′ = t ′ ξ ′ -t ′ ξ ′ + (tw + s)ξ ′ = q -(k 1 t ′ 1 + k 2 t ′ 2 )ξ ′ + ((k 1 t 1 + k 2 t 2 )w + s)ξ ′ = q + (k 1 (t 1 w -t ′ 1 ) + k 2 (t 2 w -t ′ 2 ) + s)ξ ′ therefore d T 2 ((k 1 (t 1 w-t ′ 1 )+k 2 (t 2 w-t ′ 2 )+s)ξ ′ , 0) ≤ 3cελ d(1+α)+2 q this implies that -s ∈ B 3cελ d(1+α)+2 q ((k 1 (wt 1 - t ′ 1 ) + k 2 (wt 2 -t ′ 2 ))) + ZT ′ . Notice now that the set E := {k 1 (wt 1 -t ′ 1 ) + k 2 (wt 2 -t ′ 2 ) : k 1 , k 2 ∈ Z} is discrete, so any two neighbouring points in E are at least a distance c ≥ c ′ (ξ, ξ ′ )N -1 > 0 from each other, and E + ZT ′ = E. In particular L 1 ([0, 1]\{s : d T 2 (tξ, (tw + s)ξ ′ ) ≥ ε ∀ t ≥ 0}) ≤ L 1 [0, 1] ∩ r∈E B 3cελ d(1+α)+2 q (r) (6.3) ≤ 2 c ′ (ξ, ξ ′ )N -1 3cελ d(1+α)+2 q ≤ 6c c ′ (ξ, ξ ′ ) ε(λ d(1+α)+2 q ) 2 ,
where in the last we used the inequality N < λ d(1+α)+2 q

.

We now need a number theory lemma, it is just a property on real numbers, but we state it for a sequence of real numbers, since we will apply it for a sequence.

Lemma 6.2. Let {α

n } n∈N ⊂ R + such that α n ≤ λ d(1+α)+2 q
, then there exists {ṽ n } n∈N ⊂ Q such that the following holds:

• ṽn = a n + pn qn , with p n , q n ∈ N, • a n = ⌊α n ⌋ ≤ λ α q , • q n , p n ≤ λ d(1+α)+2 q , • 0 ≤ α n -ṽn ≤ 2 λ d(1+α)+2 q , for any n ∈ N.
Proof. Fix α n , we define a n := ⌊α n ⌋ and α n := α na n ∈ [0, 1). We want to approximate α n with dyadic numbers. We define ℓ := max{N ∈ N : 2 N ≤ λ d(1+α)+2 q }. Since the dyadic intervals are such that

2 ℓ -1 i=0 i 2 ℓ , i + 1 2 ℓ = [0, 1)
there exists i = 0, .., 2 ℓ such that α n ∈ i 2 ℓ , i+1 2 ℓ , defining p n = i and q n = 2 ℓ , we get the thesis. Proposition 6.3. Consider a finite number of disjoint sets Λ j i for i = 1, .., N , j = 1, 2, as in Lemma 2.1 and their union Λ :=

N i=1 2 j=1 Λ j i ⊂ R 2 . Let {w n } n=1,..,λ d(1+α)+2 q ⊂ R satisfy w n = a n + p n q n
where a n , q n , p n are positive integers and they are less or equal than λ d(1+α)+2 q

. Then there exists a constant c 0 := c 0 ( C, Λ) > 0 with the following property: for every ξ ∈ Λ and n ∈ N there exists a ξ,n ∈ [0, 1] such that the family of curves

x ξ,n (t) := (w n t + a ξ,n )ξ with ξ ∈ Λ, n = 1, .., λ d(1+α)+2 q (6.4) satisfies d T 2 (x ξ,n (t), x ξ ′ ,m (t)) ≥ c 0 (λ d(1+α)+2 q
) 4 for every t ≥ 0, when ξ = ξ ′ . (6.5)

Proof. We fix c 0 such that Cc 0 |Λ| < 1. We define the following sets A ξ,ξ ′ ,n,m := s ∈ [0, 1] : d T 2 ((w n t + s)ξ, (w m t + s)ξ ′ ) ≥ c 0 (λ

d(1+α)+2 q ) 4
for ξ, ξ ′ ∈ Λ and n, m = 1, .., λ d(1+α)+2 q . We define

A := λ d(1+α)+2 q n,m=1 ξ =ξ ′ ∈Λ A ξ,ξ ′ ,n,m
and the thesis will follow by proving that A is not empty. We claim that L 1 (A) > 0. Using Lemma 6.1 we notice that the measure of the complement of the set A ξ,ξ ′ ,n,m satisfies L 1 (A c ξ,ξ ′ ,n,m ) ≤

Cc0(λ d(1+α)+2 q ) 2 (λ d(1+α)+2 q ) 4 = Cc0 (λ d(1+α)+2 q ) 2 . Then L 1 (A c ) = L 1   λ d(1+α)+2 q n,m=1 ξ =ξ ′ ∈Λ A c ξ,ξ ′ ,n,m   ≤ λ d(1+α)+2 q n,m=1 ξ =ξ ′ ∈Λ L 1 (A c ξ,ξ ′ ,n,m ) ≤ (λ d(1+α)+2 q ) 2 |Λ| Cc 0 (λ d(1+α)+2 q ) 2 ≤ Cc 0 |Λ| < 1,
and so A is not empty. , where the sequence {ṽ n } n=1,..,λ d(1+α)+2 q is given by Lemma 6.2 applied to the sequence α n = n k 1/p ′ and v n = ṽp ′ n (notice that the assumption α n ≤ λ d(1+α)+2 q is satisfied thanks to the bound (3.14)). We apply Proposition 6.3 to {w n } n=1,..,λ d(1+α)+2 q (notice that the assumptions are satisfied in view of Lemma 6.2) obtaining the family {a ξ,n : ξ ∈ Λ, n = 1, .., λ d(1+α)+2 q }. Finally, starting from the building blocks introduced in Section 2.3, we define W ξ,µq+1,vn (t, x) := Wξ,µq+1,vn (t, xa ξ,n ξ), Θ ξ,µq+1,vn (t, x) := Θξ,µq+1,vn (t, xa ξ,n ξ), (6.6) for any n = 1, ..., λ d(1+α)+2 q and ξ ∈ Λ. We now show that W ξ,µq+1,vn • Θ ξ ′ ,µq+1,vm (t, x) = 0 for any (x, t) ∈ T 2 × R + , (6. ) 4 . (6.9) (6.9) follows from our choice of µ q+1 = λ bγ q , because γ > 1 and b > 4(d(1 + α) + 2).

q+1,i = n≥12 ξ∈Λ [n] i χ(κ|R ℓ,i | -n)a ξ R ℓ,i
|R ℓ,i | (Θ ξ,µq+1,vn W ξ,µq+1,vn )(λ q+1 t, λ q+1 x), (6.10) as a consequence of (6.2), (3.14), the fact that χ • χ = χ and χ(κ|R ℓ,i |n) • χ(κ|R ℓ,i |m) = 0 when |n -m| > 1.

Since the average of Θ ξ,µq+1,vn W ξ,µq+1,vn which appears from the forth line of formula (4.12), in the definition of R quadr and in m is now v n ξ rather than n/κξ, the definition of Rℓ,i should now be replaced by Rℓ,i :=

Lemma 4 . 1 .

 41 For m ∈ N, S ∈ N \ {0} and n ≥ 2 we have

6. 1 .

 1 Disjointness of the supports. Set wn := µ d/p ′ q+1 v 1/p ′ n

7 )

 7 for any ξ = ξ ′ ∈ Λ, n, m = 1, .., λ d(1+α)+2 q . Indeed for any fixed t ≥ 0 one has the inclusionssupp W ξ,µq+1,vn (t, •) ⊂ B 2ρµ -1 q+1 (tw n ξ + a ξ,n ξ), supp Θ ξ ′ ,µq+1,vm (t, •) ⊂ B ρµ -1 q+1 (tw m ξ ′ + a ξ ′ ,m ξ ′ ), (6.8) hence we just need to check that B 2ρµ -1 q+1 (tw n ξ + a ξ,n ξ) ∩ B ρµ -1 q+1 (tw m ξ ′ + a ξ ′ ,m ξ ′ ) = ∅. Proposition 6.3 guarantees d T 2 (tw n ξ + a ξ,n ξ, tw m ξ ′ + a ξ ′ ,m ξ ′ )

6. 2 .

 2 Proof of the Proposition 3.2 in the case d=2. The estimates up to Section 4.2 are done in the same way by observing that ṽn = v 1/p ′ n and (n/κ) 1/p ′ are comparable up to a factor 2. In Section 4.3, we computed the product θ(p) q+1,i w (p)q+1 in (4.12) with which we were able to compensate the old error R ℓ,i (for i = 1, . . . , N ). Now this product has the form θ

  n≥12 χ(κ|R ℓ,i |n) R ℓ,i |R ℓ,i | v n ,and the obvious modification takes place for the definition of R quadr and m. Observing that|v n -n k | ≤ p ′ n k |v 1/p ′ n -( n k ) 1/p ′ | ≤ p ′ λ d(1+α)+2 q δq+2 40p ′ 2 λ d(1+α)+2 q = δq+220 the estimate (4.13) now works analogously to give |R ℓ,i -Rℓ,i | ≤ 16 20 δ q+2 . The rest of the estimates work as in Sections 4.3, 4.4 and 4.5.

  .1. Estimate on θ q+1,i L p and on inf T d θ q+1,i . We apply the improved Hölder inequality of (2.2),

	Lemma 4.1 and (3.14) to get
	θ	(p)

  .17) Property (d) is now clear from the definition of R q+1,i and the definition of ρ q+1,i . Notice that R quadr i is well defined since by (2.8) the function (Θ ξ,µq+1,n/κ W ξ,µq+1,n/κ )(λ q+1 t, λ q+1 x) -n κ ξ has 0 mean. From the second equality in (4.10) and since the average of (∂ t θ ) 1 is m i by integration by parts, we deduce that (∂ t θ

	(p)
	q+1,i (p)

this hypothesis is in [4,Lemma 7.2] where it is required that w < 10

Acknowledgements. MS has been supported by the SNSF Grant 182565. The authors wish to thank Maria Colombo for bringing the problem of non-uniqueness of integral curves to their attention and for useful suggestions.


Benyahia Samia

Meftah Boudjelal

Olivier Lézoray

Skin lesion classification using convolutional neural networks based on Multi-Features Extraction

Keywords: Feature extraction, Classification, Skin lesion, Convolutional Neural Networks

In the recent era, deep learning has become a crucial technique for the detection of various forms of skin lesions. Indeed, Convolutional neural networks (CNN) have became the state-of-the-art choice for feature extraction. In this paper, we investigate the efficiency of three state-of-the-art pre-trained convolutional neural networks (CNN) architectures as feature extractors along with four machine learning classifiers to perform the classification of skin lesions on the PH2 dataset. In this research, we find out that a DenseNet201 combined with Cubic SVM achieved the best results in accuracy: 99% and 95% for 2 and 3 classes, respectively. The results also show that the suggested method is competitive with other approaches on the PH2 dataset.

Introduction

Cancer is considered one of the leading causes of death in the world. It is estimated that the number of people diagnosed with cancer will double in the next few decades. Fortunately, some cancers have a high chance of cure through early detection, appropriate treatment, and care in the early stages. Skin lesions are an abnormal change in the tissue either on the surface of the skin or under the skin. Skin lesion usually grows in an irregular way beyond their usual boundaries as compared to the surrounding tissue. Is is primarily caused by excessive exposure to ultraviolet radiation. Skin lesions can be classified into two categories: benign skin tumors such as nevus, or malignant tumors such as melanoma that is the least common but the most harmful form of skin cancer. Until the last few years, computer-aided diagnosis was used for the early detection of skin cancer from dermoscopy images. It was based on handcrafted features extraction, such as statistical pixel-level features, shape features, texture features, and relational features from the images in order to train classical machine learning models for distinguishing the lesion from the surrounding healthy skin. Currently, such a computer-aided diagnosis is still a challenging task. Most of researchers are now investigating deep learning techniques. It is expected that they will in the next future be able to reach a performance similar to those of dermatologists, directly from dermoscopic images. Indeed, in recent years there have been many successful applications of machine learning techniques (ANN, SVM, KNN) as well as deep learning approaches and in particular of Convolutional Neural Networks (CNNs) such as AlexNet, VGGNet, ResNet, DenseNet, GoogleNet, Inception, EffcientNets. Several contributions have been proposed recently in terms of new models that have made significant improvements in the detection and the classifications of skin lesions. Hopefully, such new techniques will help to improve patient survival rates.

In this paper, we investigate the efficiency of various commonly pre-trained convolutional neural networks architectures as feature extractors with various machine learning classifiers to perform the classification for dermoscopic images from the PH2 dataset.

The rest of the paper is structured as follows. Section 2 presents some studies on the application of CNN for melanoma diagnosis on the PH2 dataset. Section 3 describes the dataset used, the various pretrained CNN used as feature extractor, the machine learning classifiers, and the evaluation metrics. Section 4 describes expriments and discusses of the obtained results. Last Section provides the conclusion of this study.

Related works

Several automated recognition methods for skin lesion images have been proposed in the last decade. In particular, recent classification approaches of skin lesions have been dominated by CNN approaches. In this section, we review some studies on the application of CNN for melanoma diagnosis on the PH2 dataset.

Ozkan et al. [START_REF] Ozkan | Skin Lesion Classification using Machine Learning Algorithms[END_REF] proposed a study on skin lesions classification based on dermoscopic images to classify images of the PH2 datasets into three classes: normal, abnormal, and melanoma. They used four different machine learning classifiers: ANN, SVM, KNN, and Decision Tree. The achieved accuracies were 92.50%, 89.50%, 82.00%, and 90.00% for ANN, SVM, KNN, and DT, respectively.

Ghasem et al. [START_REF] Ghasem Shakourian | A hierarchical structure based on Stacking approach for skin lesion classification[END_REF] proposed two-hybrid approaches to combine four heterogeneous classifiers KNN, SVM, ENN, and MLP. Their first approach was the Structure-Based on Stacking (SBS), while the other was the Hierarchical Structure Based on Stacking (HSBS). The authors have considered a preprocessing step, a segmentation step for analyzing the lesion area and have applied different feature extraction methods based on the shape, color, and texture. Finally, a classification step is used either SBS or HSBS by combining the different classifiers. The evaluation was based on the PH2 dermoscopic images with different selected features. The achieved accuracy was 96.7% for the HSBS method and 98.5% for SBS.

Ann et al. [START_REF] Salido | Using Deep Learning to Detect Melanoma in Dermoscopy Images[END_REF] used a deep learning approach to classify skin lesions in dermoscopic images. First, the authors preprocess the images to remove unwanted artifacts, such as hairs, using morphological operators and an inpainting algorithm. They classify the images using a CNN AlexNet architecture, and finally, tested the classifier using both preprocessed and unprocessed images from the PH2 dataset. The obtained accuracy for the two classes was 93%, while the accuracy for the three classes was 67.5% with preprocessed images.

Singh et al. [START_REF] Singh | Designing a Retrieval-Based Diagnostic Aid using Effective Features to Classify Skin Lesion in Dermoscopic Images[END_REF] proposed an approach for the skin lesions classification using a segmentation step by a thresholding method and the ABCD rule for feature extraction. Finally, a SVM classifier is used to make a decision. The proposed system achieved an accuracy of 92.5% on the PH2 data set.

Filali et al. [START_REF] Filali | Efficient fusion of handcrafted and pre-trained CNNs features to classify melanoma skin cancer[END_REF] proposed a skin lesion classification system based on a fusion of handcrafted features (shape, skeleton, color, and texture) and features extracted from deep learning architectures (VggNet16, ResNet18, AlexNet, and GoogLeNet). They used a Support Vector Machine (SVM) classifier to make a final decision. Their approach achieved an accuracy of 98% on the PH2 dataset.

Sanket et al. [START_REF] Sanket | Skin Cancer Classification using Machine Learning for Dermoscopy Image 1457[END_REF] proposed a skin lesion diagnosis system. The authors begin with a processing step using median and Wiener filters to remove noise, followed by a segmentation step using a watershed and morphological filters. After that, a feature extracting step is performed using the Grey Level Co-occurrence Matrix (GLCM), color and geometrical features. Finally, the classification is performed using KNN, SVM, and an ensemble method. They achieved an accuracy of 92% on the PH2 dataset.

Khalid et al. [START_REF] Khalid | Skin Cancer Classification using Deep Learning and Transfer Learning[END_REF] proposed a skin lesions classification framework based on fine-tuning a pre-trained deep learning network AlexNet. They replaced the last layer with a softmax to classify the lesion on three different classes from the PH2 dataset and achieved an accuracy of 98.61%.

Methodology and Materials

The process we propose is made up of four successive parts. The first component is the input of dermoscopic images. The second is the feature extraction part with three commonly pre-trained convolutional neural networks (CNN) architectures. From the obtained features, two databases are created: learning and test databases. The third component of the system is the learning step by a set of classifiers where each one is performed individually. Finally, the last is the validation step from the test dataset.

In addition, as shown in Figure 1, there are two kinds of extracted features with the skin lesions dataset:

1. Extracted features using the original dataset. 2. Extracted features with an augmented dataset using data augmentation techniques (will be detailed later).

PH2 dataset

The Portuguese dermatological service of the Pedro Hispano Hospital and the University of Porto have collaborated to establish the PH2 database [START_REF] Mendonan | PH2 -A dermoscopic image database for research and benchmarking[END_REF]. There are 200 dermoscopic images in the PH2 database, with 80 common nevus, 80 atypical nevus, and 40 melanomas. These are 8-bit RGB color images with a 768 x 560 pixels. The images can be obtained online from Hospital Pedro Hispano (https://www.fc.up.pt/addi/ph2%20database.html) and some are shown in Figure 2.

Convolutional Neural Networks

The convolutional neural networks (CNN) are a special type of neural networks, that possesses several convolutional layers to extract learned features directly from images in a hierarchical manner, and they can be used for classification tasks. Recently, CNN have been used to improve the performance in many different applications with various architectures such as LeNet, AlexNet, VGGNet, ResNet, DenseNet, GoogLeNet among others. Such CNN architectures can be used as pre-trained architectures that were initially trained on large datasets. In this section, we present a brief overview of such pre-trained CNNs that can be used as extractors of features for skin lesion images.

ResNet ResNet [START_REF] He | Deep Residual Learning for Image Recognition[END_REF] EfficientNet EfficientNet architecture [START_REF] Tan | EfficientNet: Rethinking Model Scaling for Convolutional Neural Networks[END_REF] is based on the principle of scaling the different dimensions such as depth, width, and image resolution of the network at the same time uniformly by using a fixed compound coefficient. Effi-cientNet is considered as a family of eight different CNN models: EfficientNet-B0 to B7. The EfficientNet-B0 represents the baseline version with an input size of 224 × 224, and it is based on the inverted bottleneck residual and squeeze-andexcitation blocks.

Classifiers

Accurate automated classification of a skin lesion in its early stages saves effort, time, and human life. The purpose of our study will be to find the most accurate classifier from several machines learning classifiers for skin lesion classification. This section provides a summary of four machine learning classifiers that will be employed for the classification of skin lesions: Artificial neural network, Support Vector Machines, K-Nearest Neighbor, and Random Forest.

Artificial Neural Network ANN are described as compound systems made up of at least two layers of neurons, an input layer, and an output layer, usually including hidden layers [START_REF] Livingstone | Artificial Neural Networks: Methods and Applications[END_REF]. Each layer contains a large number of artificial neurons that constitute an interconnected network in a weighted way. Each neuron in the network receives digital information as signals from neighboring neurons, and each of these values is assigned a particular "weight" representative of the strength of the connection.

Support Vector Machines SVMs algorithm plots the feature's value at a point in a high-dimensional space [START_REF] Vapnik | Statistical learning theory[END_REF]. Then, it performs classification by finding the hyperplane, that determines a an hyperplane that can classify the data and differentiates the space into two zones. The set of points near the hyperplane is referred to as the Support Vectors.

K-Nearest Neighbor KNN is a nonparametric method associated with only one parameter that represents the number K of nearest neighbors and a training data set [START_REF] Larose | Discovering Knowledge in Data: an Introduction to Data Mining[END_REF]. The principle of the KNN model consists of choosing the K data points closest to the point under study in order to predict its class.

Random Forest RF [START_REF] Breiman | Random Forests[END_REF] consists in constructing and training multitude of decision trees in parallel on slightly different data subsets with random variables.

Evaluation metrics

In order to evaluate the performance of our classifiers and to compare the different results of different scenarios, we used four evaluation metrics. These measures are accuracy, sensitivity, specificity, and precision.

Accuracy = T P + T N T P + F P + T N + F N (1)
Sensitivity = T P T P + F N (2) Specif icity = T N T N + F P (3)
P recision = T P T P + F P (4)

Where T P , T N , F P , F N , and refer to true positive, true negative, false positive and false negative respectively.

Experimental results and discussion

In our study for a melanoma diagnosis, two scenarios of experiments are proposed using the PH2 dataset:

-Classification of skin lesions into two types of lesions: melanoma or nonmelanoma -Classification of skin lesions into three types: melanoma, atypical nevus, or common nevus For both scenarios DenseNet201, ResNet50, EfficientB0 pre-trained CNN models are used as feature extractors, and ANN, SVM, KNN, RF classifiers used for classification. The latter are compared one to the other on the PH2 dataset. As previsouly mentioned, the PH2 dataset contains a total of 200 images that we divided into two parts.

-80% of the dataset (160 images) is used for training -20% of the dataset (40 images) is used for testing the effectiveness of the models.

For the second scenario concerning the classification into three classes, a data augmentation techniques are included with the following augmentations: rotations, zooming, shearing, flips (top-bottom, left right), skew-left-right, contrast enhancement. All experiments were performed with the same training dataset and tested with the same test set.

Table 1 and Figure 3 depict the performance results obtained for classifying PH2 dataset in two classes: melanoma or non-melanoma with DenseNet201, ResNet50, EfficientB0 CNN extractor and ANN, SVM, KNN, RF classifiers.

Table 2 and Figure 4 depict the performance results obtained for classifying PH2 dataset in 3 classes: melanoma, atypical nevus or common nevus without augmentation with DenseNet201, ResNet50, EfficientB0 CNN extractor and ANN, SVM, KNN, RF classifiers. Table 3 and Figure 5 depicts the performance results obtained for classifying PH2 dataset in 3 classes: melanoma, atypical nevus or common nevus with aug- Through the evaluation and results obtained from Table 1 to Table 3, we find that the DenseNet201 model achieved the best results as compared to other CNN models. On the other side, the cubic SVM classifier achieved the best results as compared to others classifiers ANN, KNN, RF. Through Figure 3 to Figure 5, we notice that data augmentation considerably improves the results of classification for three classes. From Table 1, Table 2 and Table 3, the statistics show that SVM classifiers perform the highest accuracy of 99% for melanoma and nonmelanoma detection while 95% for melanoma, atypical nevus, or common nevus with augmentation combined with DenseNet architecture. Table 4 and Table 5 depict the comparison of the proposed approach with the state-of-the-art models in terms of accuracy for the PH2 dataset. As it can be seen our proposed approach is competitive with the state-of-the-art.

Conclusion

In order to take advantage of deep learning models' ability to extract features from skin lesion images, several experiments were conducted with convolutional neural networks (CNN) architectures combined with various classifiers. In this work, we have proposed to use DenseNet201, ResnNet50, EfficientB0 pre-trained CNN architectures as feature extractors and the Artificial neural network (ANN), Support Vector Machines (SVM), K-Nearest Neighbor (KNN), Random Forest (RF) as classifiers to evaluate the classification of skin lesions from PH2 datasets with 2 or 3 classes. The results found show that the DenseNet201 model combined with the SVM classifier gives the better score for two classes with an accuracy of 99%, moreover for three classes with 95%, which is competitive with the actual state-of-the-art.

Fig. 1 .

 1 Fig. 1. Proposed method's flowchart.

Fig. 2 .

 2 Fig. 2. Example of PH2 skin images.

Fig. 3 .

 3 Fig. 3. Accuracy according to DenseNet201, ResnNet50, and EfficientB0 with four classifiers for melanoma or non-melanoma .

Fig. 4 .

 4 Fig. 4. Accuracy according to DenseNet201, ResnNet50, and EfficientB0 with four classifiers for melanoma, atypical nevus or common nevus without augmentation.

Fig. 5 .

 5 Fig. 5. Accuracy according to DenseNet201, ResnNet50, and EfficientB0 with four classifiers for melanoma, atypical nevus or common nevus with augmentation.

  introduces residual blocks that include skip-connection between layers where each layer feeds into the next layer and directly into the layers about 23 hops away. ResNet architecture contains a 3x3 convolution filter, global average pooling, and max-pooling layers, residual blocks, batch normalization layers followed by a fully connected layer, and softmax for the classification. ResNet architecture has several variations with a different number of layers such as ResNet18, ResNet34, ResNet50, ResNet101, and ResNet152.DenseNet DenseNet architecture[START_REF] Huang | Densely Connected Convolutional Networks[END_REF] uses shortcut connections composed of a dense block linked by transition layers. Each dense block contains convolution layers where each layer is connected to all preceding other layers. All the feature maps from previous layers are passed to all subsequent layers of the same block. DenseNet architecture contains a convolution filter, global average pooling, maxpooling layers, transition layers, dense blocks followed by fully connected layers, and softmax for the classification. DenseNet architecture proposes several vari-

ations with a different number of layers such as DensNet-121, DensNet-169, DensNet-201, and DensNet-246.

Table 1 .

 1 Summary of results for melanoma or non-melanoma

			DenseNet201			ResNet50	EfficientB0
	SVM	99	99	99	99	97.5 83.33 83.33 98.68 92.5 65.31 65.31 72.36
	KNN	97.5 98.64 98.64 87.5 95	81.98 91.98 91.98 97.5 98.64 98.64 87.5
	ANN	95	96.96 96.96 88.88 95	97.36 97.36 75.00 87.50 78.12 78.12 81.16
	RF	99	99	99	99	97.5 83.33 83.33 98.68 95	97.29 97.29 80

Acc.% Sens.% Spec.%Prec.% Acc.% Sens.% Spec.%Prec.% Acc.% Sens.% Spec.%Prec.%

Table 2 .

 2 Summary of results for melanoma, atypical nevus or common nevus without augmentation Acc.% Sens.% Spec.%Prec.% Acc.% Sens.% Spec.%Prec.% Acc.% Sens.% Spec.%Prec.% SVM 62.50 62.77 78.38 73.19 60.00 60.20 77.16 71.44 62.50 40.09 76.82 60.45 KNN 70.00 70.37 81.48 75.91 60.00 57.06 75.59 61.11 62.50 41.14 76.20 62.19 ANN 70 72.80 83.33 72.87 60.00 28.57 74.04 34.43 50.00 30.51 72.17 57.44 RF 57.50 58.69 75.30 69.59 70.00 44.75 82.09 76.50 70.00 47.78 82.38 67.35

	DenseNet201	ResNet50	EfficientB0

Table 3 .

 3 Summary of results for melanoma, atypical nevus or common nevus with augmentation.

			DenseNet201	ResNet50	EfficientB0
	SVM	95	96.67 97.22 96.29 87.50 91.25 92.77 91.02 90.00 92.91 94.16 92.59
	KNN	87.50 91.25 93.51 86.31 82.50 80	90.37 77.14 82.50 87.50 89.72 87.30
	ANN	92.50 87.77 93.88 92.50 95.00 96.29 96.00 96.66 85.00 88.88 90.90 89.16
	RF	85.00 88.75 91.57 83.57 82.50 87.50 90.92 79.39 82.50 87.50 90.04 82.47
	mentation with DenseNet201, ResNet50, EfficientB0 CNN extractor and ANN,
	SVM, KNN, RF classifiers.

Acc.% Sens.% Spec.%Prec.% Acc.% Sens.% Spec.%Prec.% Acc.% Sens.% Spec.%Prec.%

Table 4 .

 4 Results of the proposed approach compared to various approaches for PH2 dataset

	Authors	No. of classes Method	Accuracy
	Ghasem et al. [2]	2	SBS	98.50%
	Ghasem et al. [2]	2	HSBS	96.70%
	Ann et al. [3]	2	AlexNet	93.00%
	Filali et al. [5]	2	SVM	98.00%
	Sanket et al. [6]	2	SVM	92.00%
	Proposed work	2	DenseNet201+SVM 99.00%

Table 5 .

 5 Results of the proposed approach compared to various approaches for PH2 dataset

	Authors	No. of classes Method	Accuracy
	Ozkan et al. [1]	3	MLP	92.50%
	Ann et al. [3]	3	AlexNet	67.50%
	Singh et al. [4]	3	SVM	92.50%
	Khalid et al. [7]	3	AlexNet	98.61%
	Proposed work	3	DenseNet201+SVM 95.00%

Acknowledgement

This work was completed as part of the Hubert Curien Partnership (PHC) TAS-SILI cooperation program between France and Algeria under the project code 19MDU212.