Benoît Perthame
email: perthame@dma.ens.fr

Stéphane Génieys
email: genieys@math.univ-lyon1.fr

Dynamics of Nonlocal Fisher concentration points: a nonlinear analysis of Turing patterns

Keywords: Adaptive evolution, Redistributed resources, Turing instability, Nonlocal Fisher equation, Dirac concentrations, Hamilton-Jacobi equation

The so-called 'nonlocal Fisher' model takes into account an influence neighborhood for inhibition in classical Fisher ecological invasion. In this area, it has been introduced to represent front propagation with redistributed resources. More recently it has also been proposed as the simpler model exhibiting Turing instability and the biological interpretation refers to adaptive evolution. One aspect of the present paper is to propose a nonlinear analysis of these Turing patterns.

More precisely, we introduce a rescaled equation in order to take into account rare mutations (small diffusion). We analyze in which circumstances such a model exhibits stable patterns, among them the Dirac concentrations (that are interpreted as morphs in adaptive dynamics) are remarkable. We use a change of variables, similar to the phase in WKB method, that describes more accurately the phenomenon and leads to a constrained Hamilton-Jacobi equation. It allows us to interpret several features of the patterns, as the weights of the Dirac concentration points, the asymmetry variable regulating their velocities and other relevant quantities.

Model system; the Nonlocal Fisher equation

We consider the periodic solution n ε to the Nonlocal Fisher equation

           ∂nε ∂t (t) = ε∆n ε + 1 ε n ε (1 -Φ * n ε ), 0 ≤ x ≤ 1, n ε (t, 0) = n ε (t, 1), ∂nε(t,0) ∂x = ∂nε(t,1) ∂x , n ε (t = 0, x) = n 0 ε (x) ≥ 0, (1)
with an initial data n 0 ε exhibiting a Dirac mass behavior as we describe it later in more details. The convolution kernel Φ, among other forthcoming properties, has bounded support,

Φ ≥ 0, Φ = 1, Φ(x) = 0, for x / ∈ [-b, b] (b < 1 2 ). (2)
This is a rescaled version of the competition model introduced by Génieys, Volpert and Auger [START_REF] Génieys | Adaptive dynamics: modeling Darwin's divergence principle[END_REF]. The authors interpret the variable x as a physiological trait, the diffusion term models mutations along with usual genetic population laws (and the scaling means rare mutations with a long range effect). The convolution term mimics competition for resources between individuals whose traits are close enough. But the Nonlocal Fisher equation was also introduced as an ecological model; in the correct range of parameters, it describes an invasion front for a sysytem with redistributed resources, see [START_REF] Gourley | Travelling front of a nonlocal Fisher equation[END_REF] and the references therein.

The interpretation and the results in [START_REF] Génieys | Adaptive dynamics: modeling Darwin's divergence principle[END_REF] make a direct link between evolution theory and Turing's instability. Indeed, this model exhibits both Turing's instability (under specific properties on the kernel Φ that we discuss in Section 5.1) and a phenomenon of evolutionary branching: a monomorphic population (i.e. an initial data with a single Dirac mass) may become dimorphic in order to lower the competition. It is an illustration of Darwin's divergence principle, which is a mechanism of branching due to the competition and can arise even in a constant homogeneous environment [START_REF] Darwin | On the origin of species by means of natural selection[END_REF]. This principle is illustrated numerically by using the cellular automata approach ( [START_REF] Atamas | Self-organization in computer simulated selective systems[END_REF] for instance) or differential systems ( [START_REF] Diekmann | The dynamics of adaptation : an illuminating example and a Hamilton-Jacobi approach[END_REF][START_REF] Geritz | Evolutionary singular strategies and the adaptive growth and branching of the evolutionary tree[END_REF][START_REF] Perthame | Transport equations in biology[END_REF]). It is mathematically investigated by using a probabilistic approach in [START_REF] Champagnat | Unifying evolutionary dynamics: from individual stochastic processes to macroscopic models[END_REF] and game theory in [START_REF] Hofbauer | Evolutionary games and population dynamics[END_REF]. Its relevance in physiological contexts is discussed in [START_REF] Kupiec | Ni Dieu ni gène. Pour une autre théorie de l'hérédité[END_REF].

The re-scaling we have introduced allows us to identify two length scales: ε for the activation and b for the inhibition, typical of Turing's instability [START_REF] Turing | The chemical basis of morphogenesis[END_REF][START_REF] Meinhardt | Models of biological pattern formation[END_REF]. Even though the model is very simple it allows us to recover a general principle in the topic. For long rang activation and short range inhibition, i.e. b ≪ ε, we expect front propagation. As b → 0 we recover the classical fisher equation which is well-known to exhibit a traveling wave (see [START_REF] Murray | [END_REF][START_REF] Berestycki | Reaction-Diffusion Equations and Propagation Phenomena[END_REF] for instance). For b small, traveling waves in a nonlocal Fisher equation is studied in [START_REF] Gourley | Travelling front of a nonlocal Fisher equation[END_REF].

In the other regime short range activation and long range inhibition, we expect pattern formations. This arises from an instability that occurs indeed for certain competition kernels Φ characterized by the fact that their Fourier transform is 'negative enough', as proved in [START_REF] Génieys | Adaptive dynamics: modeling Darwin's divergence principle[END_REF] (see also Section 5.1). See also the appendix of the present paper for an estimation of the leading wave length. Even though Turing's instability, and thus pattern formation, are more traditional for systems of parabolic PDEs, notice that convolution models are also commonly used in this area, [START_REF] Lefever | On the origin of tiger bush[END_REF][START_REF] Desvillettes | Infinite Dimensional Reaction-Diffusion for Population Dynamics[END_REF].

Our purpose is to study, in this situation when ε vanishes, wether one can expect the asymptotic expansion

n ε (t, x) ≈ I i=1 ρ i (t) δ(x -x i (t)).
(

) 3
This is indeed the same asymptotic regime than the formalism for adaptive dynamics introduced in [START_REF] Diekmann | The dynamics of adaptation : an illuminating example and a Hamilton-Jacobi approach[END_REF] and also used in [START_REF] Barles | Concentrations and constrained Hamilton-Jacobi equations arising in adaptive dynamics[END_REF][START_REF] Carrillo | Adaptive dynamics via Hamilton-Jacobi approach and entropy methods for a juvenile-adult model[END_REF] for various models arising in population biology. In this language ( [START_REF] Geritz | Evolutionary singular strategies and the adaptive growth and branching of the evolutionary tree[END_REF][START_REF] Diekmann | Beginner's guide to adaptive dynamics[END_REF]), this asymptotic can be interpreted as polymorphism and I denotes the number of different traits that can be represent. The interesting feature here is that this number I is not fixed in advance by opposition to the chemostat problem in [START_REF] Diekmann | The dynamics of adaptation : an illuminating example and a Hamilton-Jacobi approach[END_REF].

The main questions we wish to answer here are • Is this expansion (3) always true when Turing's instability occurs?

• Can one predict the densities ρ i of each morph in terms of the competition kernel Φ? • Can one predict the 'velocity' of the concentration points x i (t) in terms of the asymmetry of the competition kernel Φ? • For steady states situations (with symmetric competition kernel Φ) can one predict the spacing between the concentration points x i (t), and is it given by the wave length of the linearized equation?

• What is the shape of n ε in this regime?

This paper uses the Hamilton-Jacobi formalism that was introduced in [START_REF] Diekmann | The dynamics of adaptation : an illuminating example and a Hamilton-Jacobi approach[END_REF] to predict quantitatively the answers to these questions and explain some apparent paradoxes: one of them is that very 'asymmetric' competition kernels can lead to slower velocities of the points x i (t) than seemingly 'less asymmetric' kernels. Also, several rules have been derived, still validated on specific examples, for Turing instabilities based on the linearized eigenvalues. Here we recover for instance the rule, already explained, on relative ranges for activators and inhibitors related to waves or patterns. However the linearized eigenvector do not necessarily give a good hint on the wave length of the patterns.

The outline is as follows. We first present the asymptotic theory based on Hamilton-Jacobi formalism. Then we analyze the steady states and answer to most of the questions mentioned above. After that, we give the main parameter that determines the dynamics of the Dirac concentration points. We conclude with several mathematical proofs that we gather in a final section.

Asymptotic analysis

In this section, we recall the method introduced in [START_REF] Diekmann | The dynamics of adaptation : an illuminating example and a Hamilton-Jacobi approach[END_REF] which allows to analyze the limit ε → 0 and study the concentration effects. It is based on a classical idea which consists in introducing a 'phase' along with the WKB method for oscillations. Here the parabolic aspect leads us to a real phase as in 'front propagation' analysis introduced in [START_REF] Evans | A PDE approach to geometric optics for certain semilinear parabolic equations[END_REF][START_REF] Barles | Wavefront propagation for reaction diffusion systems of PDE[END_REF].

A simple example that motivates the forthcoming analysis is the 'zero temperature maxwellian' 2ε) . It converges to a Dirac mass at x = 0, but it is easier to pass to the limit in the quantity

n ε = 1 √ 2πε e -|x| 2 /(
ε ln(n ε ) = -|x| 2 2 - ε 2 ln(2πε) ---→ ε→0 - |x| 2 2 ≤ 0.
In the same way, one can handle several Dirac masses. We consider for example

n ε = 1 √ 2πε I i=1 ρ i θ i e -|x-x i | 2 /(2εθ i ) ---→ ε→0 I i=1 ρ i (t) δ(x -x i (t))
and

ε ln(n ε ) ---→ ε→0 -min i |x -x i | 2 2θ i ≤ 0.
With this in mind, we set, following [START_REF] Diekmann | The dynamics of adaptation : an illuminating example and a Hamilton-Jacobi approach[END_REF][START_REF] Barles | Concentrations and constrained Hamilton-Jacobi equations arising in adaptive dynamics[END_REF][START_REF] Carrillo | Adaptive dynamics via Hamilton-Jacobi approach and entropy methods for a juvenile-adult model[END_REF],

n ε (t, x) = e ϕε(t,x)/ε , n 0 ε (x) = e ϕ 0 ε (x)/ε
with the property ϕ 0 ε (x) ---→ ε→0 ϕ 0 (x) ≤ 0, max ϕ 0 = 0.

Then the points where ϕ(t) vanishes are the Dirac concentrations but their weights are not known by this construction. Inserting this ansatz in equation ( 1), we find

   ∂ϕε ∂t (t) = ε∆ϕ ε + |∇ϕ ε | 2 + (1 -Φ * n ε ), ϕ ε (t = 0, x) = ϕ 0 ε (x), (4)
From Theorem 5.3 (see the technical results in the last section), we know that, after extraction of a subsequence if necessary, ϕ ε → ϕ (uniformly, locally in time).

The mass control from below and from above, proved in Theorem 5.1 gives us a first property of the limit ϕ. Namely max

x ϕ(t, x) = 0 ∀t ≥ 0. (5)
Indeed, would this max be negative, then n ε would go extinct which is not true; would this max be positive, then by continuity this would be true in an interval and n ε would blow-up contradicting the finite mass property.

Passing to the limit as ε vanishes in (4), we find the Hamilton-Jacobi equation (H.-J. eq. in short)

   ∂ϕ ∂t (t) = |∇ϕ| 2 + r(t, x), ϕ(t = 0, x) = ϕ 0 (x), (6)
and

r(t, x) = lim ε→0 1 -Φ * n ε , (7)
notice that when Φ ∈ C(R) this is a uniform limit and r(t, •) ∈ C(0, 1) (continuous functions) for all times. It is a standard point that the limiting solution ϕ is not smooth (it is only Lipschitz continuous) and satisfies the H.-J. eq. in the viscosity sense of Crandall-Lions (see [START_REF] Crandall | User's guide to viscosity solutions of second order partial differential equations[END_REF][START_REF] Fleming | Controlled Markov processes and viscosity solutions[END_REF]). There is a specific difficulty here because the time dependency of r in time is just L ∞ but this is handled also with arguments which are standard by now and we refer the reader to [START_REF] Barles | Concentrations and constrained Hamilton-Jacobi equations arising in adaptive dynamics[END_REF][START_REF] Perthame | Transport equations in biology[END_REF] for uniqueness and further mathematical analysis of this problem.

Assuming n ε converges to the sum of Dirac masses as expressed in (3), we find in ( 5) and ( 7):

ϕ(t, x i (t)) = 0, i = 1, ..., I, (8)
r(t, x) ≡ 1 - I i=1 ρ i (t) Φ(x -x i (t)). (9)
We recall the interpretation that ϕ satisfies the constrained H.-J. eq. ( 5), [START_REF] Champagnat | Unifying evolutionary dynamics: from individual stochastic processes to macroscopic models[END_REF]. To fulfill the constraint on the max imposes to introduce Lagrange multipliers, i.e., free parameters that adapts the right-hand side to these constraints, and these are the weights ρ i (t). Therefore the system (5), ( 6), ( 9) can be seen as a coupled problem. As we know it from [START_REF] Diekmann | The dynamics of adaptation : an illuminating example and a Hamilton-Jacobi approach[END_REF] there might be several solutions when the max in ( 5) is attained at multiple points and one expects that an additional information has to be added in this case (and it is an open question to find a convenient criteria). Only in the case where a single maximum point is possible, one can prove uniqueness ( [START_REF] Barles | Concentrations and constrained Hamilton-Jacobi equations arising in adaptive dynamics[END_REF]).

From the relations ( 8) and ( 5), we deduce that ∂ ∂t ϕ(t, x i (t)) = ∇ϕ(t, x i (t)) = 0. Therefore we find the fundamental relation

r(t, x i (t)) = 0 1 ≤ i ≤ I. ( 10
)
This formula is standard in the theory of adaptive dynamics; r is called the invasion exponent. The quantity r x also appears (see the canonical equation [START_REF] Lefever | On the origin of tiger bush[END_REF] for instance) and is called the selection gradient and the points where it vanishes are the so-called singular points.

Steady state

We begin with the simpler situation when the potential Φ is even. Then, the steady states solutions to (1) are reached as it will be explained in Section 4 (but for 'asymmetric potentials' we cannot hope for steady states). Then, we have the alternative to study the time asymptotic patterns that satisfy equation

   -ε∆n ε = 1 ε n ε (1 -Φ * n ε ), 0 ≤ x ≤ 1, n ε (0) = n ε (1), ∂nε(0) ∂x = ∂nε(1) ∂x . (11)
Then, the limiting equation is simpler and reads

-|∇ϕ| 2 = r(t, x) = 1 - I j=1 ρ j Φ(x -x j ) =    ≤ 0 f or x = x i , = 0 f or x = x i . ( 12
) max x ϕ(x) = 0 = ϕ(x i ), i = 1, ..., I. (13)
We analyze now the numerical results obtained with different even competition kernels in view of this asymptotic equation. We have distinguished two cases, both satisfy the Turing instability criteria, but have different behaviors in terms of pattern formation. We refer to them as well-or ill-differentiated.

Well-differentiated competition

The simplest situation is when we choose a competition kernel exhibiting stronger competition for slightly different traits x (a behavior that is largely accepted and was observed already in [START_REF] Darwin | On the origin of species by means of natural selection[END_REF]). More precisely, we take for real numbers a > 0, b > 0,

Φ(x) = Φ(0) (1 + a |x| b ) 1 {|x|≤b} , Φ = 1. ( 14
)
On this class of kernels, we have checked numerically that Dirac concentrations are always obtained, and we can discuss the asymptotic prediction made in Section 2.

• (Location of the concentration points) Firstly, we observ that the patterns, even though they concentrate as Dirac masses for small ε, depend upon the initial data. Figure 1 exhibits four or five peaks depending on the initial data. Secondly, we can also check for the example a = 2 (but we have tested successfully other values which are not shown here) that the Dirac locations corresponds to the points where ∇ϕ vanishes; see [START_REF] Atamas | Self-organization in computer simulated selective systems[END_REF] with the well differentiated even potential ( 14) with b = .15, a = 2, this corresponds to the third line of Table 1. We represent the density n ε (left) and the phase ϕ ε (right). The numerics is performed with 3000 points, ε = 2.10 -4 . n ε and ϕ ε in Figure 2.

• (Validity of the H.-J. equation) Thirdly, we check the validity of the asymptotic limit ( 5), [START_REF] Champagnat | Unifying evolutionary dynamics: from individual stochastic processes to macroscopic models[END_REF]. The Figure 3 depicts the potential Φ * n ε compared to the computed value of 1 + |∇ϕ ε | 2 . The undershoots are just computational effects when computing the numerical derivative of the solution ϕ ε .

• (Weights of the Dirac masses) Observe that, in the limiting expression Φ * n ε = I j=1 ρ j (t) Φ(x -

x j (t))dis, some overlapping is necessary for b = .15 because the values Φ * n ε (x) = 0 are unstable. But it occurs on a sufficiently small set in such a way that at the points x i there is no overlaping. But also, from [START_REF] Evans | A PDE approach to geometric optics for certain semilinear parabolic equations[END_REF] at the points x i , we have ∇ϕ(x i ) = 0 (maximum point of ϕ). Therefore, we predict that the weights ρ i should satisfy ρ i Φ(0) = 1. The Table 1 shows the computed values ρ i and compare them to 1/Φ(0). A very good agreement is achieved (to the expense of a very fine grid).

Ill-differentiated competition

A very interesting case is that initially proposed in [START_REF] Génieys | Adaptive dynamics: modeling Darwin's divergence principle[END_REF] where the competition kernel is given by

Φ(x) = Φ(0) 1 {|x|≤b} , Φ = 1. (15)
It is proved in [START_REF] Génieys | Adaptive dynamics: modeling Darwin's divergence principle[END_REF] that it exhibits Turing's instability (see also Section 5.1 for an explicit computation in our framework). It appears, from numerical simulations, that several stable steady states can be produced depending upon the initial data and not all of them converge to Dirac masses as ε → 0. Here, we only focus on patterns that remain smooth as ε vanishes. Figure 4 shows results indicating that, in place of (3), one should rather expect

n ε (t, x) ≈ I i=1 u(x -x i ),
where the profile u is a continuous even function with bounded support such that u = 1. Denoting by c the length of this support, numerical simulations suggest that

c + b = 1 I . ( 16
)
This means that the number of patterns is maximizing space occupation with the property that the potential terms Φ * u(x -x i ) do not overlap on the support of n ε , i.e. there is a single i such that Φ * u(x -x i ) = 0. In the case of Dirac concentrations, as mentioned earlier some overlapping, away from the x i 's, is also necessary to avoid the state Φ * u(x) = 0 but there is no free parameter c that allows for such an equality since c = 0.

Our claim here is that such a steady state is possible only if c < b. Indeed similar to the case of u(y)dy does not depend on x, which is compatible with equation ( 17).

Hence we have obtained 0 < c < b which combined with ( 16) yields

1 2b < I < 1 b . ( 18
)
Since I ∈ N, only a finite number of solutions are possible. In the situation of figure 4, we still take b = .15 and the numerics is in accordance with the conclusion that steady-states solutions may exist with I = 4, 5 or 6. These three possibilities occur, still depending upon the initial data (as in Figure 1). Notice that for the well-differentiated potentials already presented in Subsection (3.1) we have seen that stable steady states are Dirac masses. Hence the free parameter c is not available to justify equation [START_REF] Gourley | Travelling front of a nonlocal Fisher equation[END_REF]. However the numerical simulations of Figure 1 suggest that the possible numbers of peaks are still determined by equation [START_REF] Kupiec | Ni Dieu ni gène. Pour une autre théorie de l'hérédité[END_REF].

Asymmetry and dynamic of concentration points

We now study the more general case of 'asymmetric' potentials that allow for motion of the concentration points. The question here is to determine which 'asymmetry' criteria controls this motion. A simple and somewhat natural direction would be to measure it by the 'oddness' of the kernel Φ using the quantity A[Φ] = xΦ(x)dx. Our purpose here is to show that the asymptotic theory of Section 2 gives the answer, and it turns out that this 'oddness' quantity A[Φ] is not the correct one.

Following [START_REF] Diekmann | The dynamics of adaptation : an illuminating example and a Hamilton-Jacobi approach[END_REF][START_REF] Perthame | Transport equations in biology[END_REF], the H.-J. eq. ( 6)-( 7) contains information on the dynamic of the concentration points x i (t) defined in the asymptotic formula (3) because they satisfy φ(t, x i (t)) = 0. Namely, we recover the velocities with the formula

ẋi (t) = -D 2 ϕ(t, x i (t)) -1 ∂r ∂x (t, x i (t)). ( 19
)
and we know from the definition [START_REF] Crandall | User's guide to viscosity solutions of second order partial differential equations[END_REF], that

∂r ∂x (t, x) = - I j=1 ρ j (t) Φ ′ (x -x j (t)). ( 20
)
As a consequence, the correct measure of asymmetry for the speed of the concentration points x i (t) is given by the quantity Φ ′ (0). These points move only when Φ ′ (0) = 0 (and thus Φ is not even). We check numerically that x i (t) moves according to the value of Φ ′ (0). We have used the three following potentials

Φ 1 (x) = a 1 (x + .5b) + 1 {|x|≤b} , Φ = 1, Φ 2 (x) = a 2 (x + b) + 1 {|x|≤b} , Φ = 1, Φ 3 (x) = 1 {-.5b≤x≤b} , Φ = 1.
We have

Φ ′ 1 (0) > Φ ′ 2 (0) > Φ ′ 3 (0)
= 0, and thus we can expect that the concentration points move faster for Φ 1 than for Φ 2 and that they do not move with Φ 3 . This behavior is indeed obtained as depicted in Figure 5.

Mathematical proofs and numerical algorithms

In this section we gather several mathematical results that have been used throughout the paper. They mostly concern a priori estimates. We also describe, for the sake of completeness, the numerical algorithms that we have used.

Turing's instability and estimation of the typical wave length

For the sake of completeness we reproduce here the computation for deriving the Turing instability condition in model (1), following [START_REF] Génieys | Adaptive dynamics: modeling Darwin's divergence principle[END_REF].

Linearizing the equation (1) around the steady state n = 1, we obtain

   ∂r ∂t (t) -ε∆r = -1 ε Φ * r, 0 ≤ x ≤ 1, r(t, •) 1-periodic.
As usual one tries to find solutions with exponential growth, which means eigenvectors

   λr(x) -ε∆r = -1 ε Φ * r, 0 ≤ x ≤ 1, r(•) 1-periodic, λ > 0. (21)
It is natural to decompose a possible solution r in Fourier series

r(x) = n∈Z r(n)e 2iπnx , r(n) = 1 0 r(x)e -2iπnx dx.
Then, equation [START_REF] Murray | [END_REF] becomes

r(n) λ + ε(2πn) 2 + 1 ε Φ(n) = 0.
Consequently, the eigenvalues are given by

λ = -ε(2πn) 2 - 1 ε Φ(n). ( 22
)
Turing instability occurs when some Fourier coefficient of Φ is negative and the leading positive eigenvalue is given for ε small, by

λ ε ≈ - 1 ε min n Φ(n) = - 1 ε Φ(n 0 ).
The corresponding eigenvector is e 2iπn 0 x , which exhibits a wave length L 0 = 1 |n 0 | . We can specify the example Φ = 1 2b 1 {|x|≤b} , with b < 1 2 and we find that

Φ(n) = b -b e -2iπnx dx = 1 2πnb
sin(2πnb).

The smallest frequency such that instability occurs, Φ(n 0 ) = 0, is given here by n 0 = E( 1 2b ) + 1 and thus L 0 = 1/[E( 1 2b ) + 1] (which is close to 2b when b is small). This is also the largest possible length of patterns as expressed in [START_REF] Kupiec | Ni Dieu ni gène. Pour une autre théorie de l'hérédité[END_REF] but in practice a different wave length is obtained (see Section 3).

Mass control for the Nonlocal Fisher Eq.

Theorem 5.1 Assume (2). Then the solution to the Nonlocal Fisher equation (1) has saturated growth and non-extinction. More precisely the total mass M ε (t) = n ε (t, x)dx satisfies

min M ε (0), 1 Φ ∞ ≤ M ε (t) ≤ max M ε (0), 1 Φ m c 2 ,
where c and Φ m are any numbers such that Φ ≥ Φ m on (-2c, 2c). Moreover, as t → ∞

lim inf M ε (t) ≥ 1 Φ ∞ , lim sup M ε (t) ≤ 1 Φ m c 2 .
Remark 5.2 In view of the numerical experiments, one expects that L p norms are unbounded for p > 1 and ε small in the case of a well-differentiated kernel Φ.

Proof. (i) Non-extinction. Integrating in x the equation [START_REF] Atamas | Self-organization in computer simulated selective systems[END_REF], we obtain the relation

d dt M ε (t) = 1 ε 2 M ε (t) -n ε Φ * n ε , (23)
and we have

Φ * n ε ≤ Φ ∞ M ε (t).
Therefore, since n ε ≥ 0, we arrive to the differential inequality

d dt M ε (t) ≥ M ε (t) ε 2 (1 -M ε (t) Φ ∞ ) ,
which proves the lower bound.

(ii) Limited growth. Consider a value c > 0 as in the statement of Theorem 5.1, and an interval I of length c where

I n ε dx ≥ c M ε (t).
Then, still with the notations in the statement of Theorem 5.1, we have the lower bounds

n ε Φ * n ε ≥ Φ m |x-y|≤2c n ε (t, x) n ε (t, y) dxdy ≥ Φ m I×I n ε (t, x) n ε (t, y) dxdy ≥ Φ m c 2 M ε (t).
Therefore, we obtain from ( 23)

d dt M ε (t) ≤ M ε (t) ε 2 1 -c 2 M ε (t)Φ m ,
and again the upper bound follows directly.

5.3

Estimates for the H.-J. Eq.

Theorem 5.3 We assume (2), Φ ∈ W 1,∞ (R) and M ε (0) = M 0 < ∞, and that ϕ 0 ε and ∇ϕ 0 ε are bounded. Then, for all time t > 0, the solution to equation (4) satisfies the a priori bounds

ϕ ε (t, •) L ∞ (0,1) + ∂ ∂x ϕ ε (t, •) L ∞ (0,1) ≤ C(t),
and, after extraction of a subsequence, ϕ ε → ϕ ∈ W 1,∞ (0, T ) × (0, 1) (uniformly locally in time) and ( 6) holds (in the viscosity sense).

Proof. Differentiating equation ( 4), we obtain the equation for w

= ∇ϕ ε , ∂w ∂t (t) = ε∆w + 2∇ϕ ε • ∇w -∇Φ * n ε . (24)
Here, we can upper bound

|∇Φ * n ε (t)| ≤ ∇Φ ∞ M (t),
where M (t) is an upper bound of the total mass M (s) for times 0 ≤ s ≤ t which is controlled thanks to Theorem 5.1. Therefore, by the maximum principle, we have

|w(t, x)| ≤ max x∈(0,1) |w(t = 0, x)| + ∇Φ ∞ M (t)t.
This uniform estimate on the gradient gives a uniform estimate on ϕ ε in L ∞ (0, T ) × (0, 1) when coming back to the equation ( 4) because it boils down to heat equation with a bounded right-hand side. And thus, the L ∞ bounds on ϕ ε and ∂ ∂t ϕ ε , as stated in Theorem 5.3, are proved.

Then we can obtain a local uniform L 2 estimate on ∂ϕε ∂t . To do that, we multiply for instance the equation ( 4) by ∂ϕε ∂t , and integrate. We obtain, denoting by C an absolute bound for

|∇ϕ ε | 2 +1-Φ * n ε 1 0 | ∂ϕε(t) ∂t | 2 dx ≤ ε 1 0 ∂ϕε(t) ∂t ∆ϕ ε (t) + C 1 0 | ∂ϕε(t) ∂t |dx ≤ -ε 2 d dt 1 0 |∇ϕ ε (t)| 2 dx + C 1 0 | ∂ϕε(t)
∂t |dx. Therefore, when integrating in time and using Cauchy-Schwarz inequality for the second term on the right, we find

T 0 1 0 | ∂ϕ ε (t) ∂t | 2 dx ≤ ε 2 1 0 |∇ϕ ε (t = 0)| 2 dx + C √ T T 0 1 0 | ∂ϕ ε (t) ∂t | 2 dx 1/2 .
This proves a uniform bound on the quantity

T 0 1 0 | ∂ϕε(t) ∂t | 2 dx.
Together with the bound on the x derivative, this proves the compactness of ϕ ε for the uniform topology.

Passing to the limit in viscosity sense and almost everywhere is standard then. The constraint max ϕ(t, •) = 0 follows from the mass constraint (see [START_REF] Diekmann | The dynamics of adaptation : an illuminating example and a Hamilton-Jacobi approach[END_REF][START_REF] Barles | Concentrations and constrained Hamilton-Jacobi equations arising in adaptive dynamics[END_REF][START_REF] Perthame | Transport equations in biology[END_REF]).

Numerical methods

We have used two different numerical schemes that we describe now. Both of them are very simple and based on finite differences.

In the first scheme, the discretization of equation ( 1) is based on a time splitting of the reaction term n(1 -Φ * n) and the differential term. Considering a number N of points in x, we approximate the solution by a discrete vector n k i for i = 1, ..., N and k the label for discrete time. We set ∆x = 1/N the space stepping and ∆t the time stepping.

In order to avoid strong limitations on the time step, we use an exact resolution of the reaction term

n k+1/2 i = n k i exp ∆t ε (1 -Φ * n k ) ,
where the convolution Φ * n is computed according to

Φ * n ≈ J M j=-J M Φ j n i-j ,
where J M = b * N . We use a three point implicit (or explicit when ε is small enough) scheme for the differential term

n k+1 i = n k+1/2 i + ε ∆t 2∆x 2 [n k+1 i+1 + n k+1 i-1 -2n k+1 i ].
In the second scheme, we use an explicit scheme for the reaction as well as for the diffusion. The convolution is computed by using the trapeze formula:

Φ * n ≈ J M -1 j=0
(Φ j n i-j + Φ j+1 n i-j-1 )/2 + (Φ -j n i+j + Φ -j-1 n i+j+1 )/2. This gives a better accuracy and allows to check that the number of discretization points is enough to resolve the singularities.

Finally, notice that for accuracy reasons it is necessary to ensure some kind of CFL condition

∆t ε ≪ 1. ( 25
)
This is imposed in order to ensure that the reaction term is well resolved. In practice we can choose ∆t ε = .1. As a consequence, when ε is small, the CFL condition for the diffusion term, namely ε ∆t 2∆x 2 ≤ 1, can be achieved. Then the explicit scheme for the diffusion is stable.

Figure 1 :

 1 Figure 1: Well-differentiated potential (14) with a = .015, b = .15, ε = 10 -4 . Initial conditions are: 4 peaks evenly spaced (left), 5 peaks evenly spaced (right).

Figure 2 :

 2 Figure2: (Symmetric potential) Numerical solution of the problem[START_REF] Atamas | Self-organization in computer simulated selective systems[END_REF] with the well differentiated even potential (14) with b = .15, a = 2, this corresponds to the third line of Table1. We represent the density n ε (left) and the phase ϕ ε (right). The numerics is performed with 3000 points, ε = 2.10 -4 .

Figure 3 :

 3 Figure 3: (Symmetric potential) Same as Figure 2. We depict the potential Φ * n ε (left) and a centered evaluation of 1 + |∇ϕ ε | 2 (right).

  a

Figure 4 :

 4 Figure 4: Ill-differentiated potential (15) with b = .15, ε = 10 -4 . Initial conditions are: 4 peaks evenly spaced (left), 5 peaks evenly spaced (right).

Φ-c/ 2 u 2 -c/ 2

 222 Dirac masses, for x ∈ [x i -c/2, x i + c/2], |∇ϕ| 2 = r(t, x) = 1 -I j=1 * u(x -x j ) = 0 because n ε remains positive and thus ϕ vanishes. Since the terms Φ * u(x -x j ) are non-overlapping on [x i -c/2, x i + c/2], this means that Φ * u(x) ≡ 1, for x ∈ [-c/2, c/2]. (17) But for x ∈ [-c/2, c/2], we can compute Φ * u(x) = min(x+b,c/2) max(x-b,-c/2) u(y)dy. If one had c > b, for x large enough (say x ∈ [b-c/2, c/2]), we would also have Φ * u(x) = which depends on x and contradicts equation (17). Similarly, for x small enough (x ∈ [-c/2, -b+c/2]), Φ * u(x) = x+b (y)dy which contradicts equation (17). On the other hand if c < b, Φ * u(x) = c/

Figure 5 :

 5 Figure 5: (Asymmetric potential) Dynamic of the concentration points with various asymmetric potentials. In the upper left the value Φ ′ (0) is larger than in the upper right (a 1 and a 2 respectively). The lower pictures are obtained with an asymmetric potential satisfying Φ ′ (0) = 0 (potential Φ 3 ) with 3000 points on the left and 6000 points on the right. The abscissae are x and the ordinates are time.

Table 1 :

 1 Comparison between the predicted values of the weights of the Dirac mass and the weights obtained by numerical simulations. We have taken the parameter b = .15 and left: ε = 2 10 -4 , 3000 discretization points, right: ε = 5 10 -5 , 10000 discretization points,.

		1/Φ(0)	numerical weight	a	1/Φ(0)	numerical weight
	0	.3	.300	0	.3	.300
	1	.45	.442	1	.45	.447
	2	.6	.587	2	.6	.594
	3	.75	.724	3	.75	.747


Nicolas Totaro

Sandra Forget

Jean-Louis Guyader

iPTF methods: How Green's identity and FEM solver can be used for acoustic inverse methods

Keywords: no, 43, 35, Sx, 43, 60, Pt

may come

iPTF methods: How Green's identity and FEM solver can be used for acoustic inverse methods

Introduction

Source identication methods like Near-eld Acoustic Holography [START_REF] Williams | Sound source reconstructions using microphone array[END_REF][START_REF] Maynard | Neareld acoustic holography: I. Theory of generalized holography and the development of NAH[END_REF] aim to locate a source of noise or to reconstruct its vibratory eld. The diversity of applications in this eld of research has caused to develop several dierent methods with their own domain of validity (acoustic environment, frequency range, shape of the source, etc.) and their own experimental setup (synchronous or sequential measurements, microphones antennas, etc.). Despite the number of existing methods, there are still some issues for industrial applications. Measurements in non-anechoic conditions, reconstruction of elds on 3D sources, eld separation are some of the diculties that still slow down the industrial (c) European Acoustics Association applications. The present paper deals with a way to develop new source identication methods based on the use of Green's identity. Green's identity is a well-known technique to solve problems in acoustics. Two dierent developments will be shown here to illustrate the possibilities of the Green's identity in the framework of inverse methods. These developments, called u-iPTF [START_REF] Aucejo | Identication of source velocities on 3D structures in non-anechoic environments: Theoretical background and experimental validation of the inverse patch transfer functions method[END_REF][START_REF] Totaro | Sound elds separation and reconstruction of irregularly shaped sources[END_REF] (inverse Patch Transfer Functions with uniform boundary conditions) and m-iPTF (inverse Patch Transfer Functions with mixed boundary conditions) have several advantages: (i) the properties of the Green's identity ensure the deconnement of the source from its acoustic environment (ii) the judicious choice of the associate problem opens a lot of development possibilities (iii) the use of a numerical solver permits to treat irregularly shaped sources (iv) the source is fully characterized (pressure, normal velocity, normal intensity elds, acoustic power, radiation eciency). These examples of developments using Green's identity will be presented and illustrated by a numerical example and experimental applications.

Green's identity in acoustics

Let's take the classical case presented in Fig. 1 of a baed structure radiating in any acoustic environment. In this acoustic environment, some other (stationnary) sources may contribute to the noise measured at point N. To compute the pressure at point N, one can use the Green's identity. This well known method is extensively used in acoustics to compute the pressure anywhere in an acoustic volume knowing pressure and particle velocity elds at boundaries. In the case of Fig. 1 there is no strictly speaking acoustic volume, just an (semi-innite) acoustic environment. In fact, in Green's identity, there is no need of a "real" acoustic volume, a "virtual" volume suits well. Let's dene a "virtual" acoustic volume Ω containing the point N and the vibrating surface Σ and excluding other sources as presented in Fig. 2. Obviously, the pressure and particle velocity elds on the "virtual" surface Σ are due to the contribution of the radiating surface Σ and to the contribution of all other sources, reections, absorption outside the volume Ω. Let's consider that the pressure and the particle velocity are continuously known on surfaces Σ (vibrating surface), Σ (physically rigid surface) and Σ ("virtual" surface).

Then, the Green's identity writes

Ω (∆p(N ) + k * 2 p(N ))Φ(N )dN = Ω (∆Φ(N ) + k * 2 Φ(N ))p(N )dN + Σ p(Q) ∂Φ(Q) ∂n -Φ(Q) ∂p(Q) ∂n dQ + Σ p(Q ) ∂Φ(Q ) ∂n -Φ(Q ) ∂p(Q ) ∂n dQ + Σ p(Q ) ∂Φ(Q ) ∂n -Φ(Q ) ∂p(Q ) ∂n dQ , (1)
Figure 2. Denition of a "virtual" volume Ω and its boundary surfaces Σ, Σ and Σ for the application of Green's identity.

where k * = k(1 -jη) is the complex wavenumber, with η is the damping loss factor and n is the outward normal. The boundary conditions are prescribed by the pressure and particle velocity on the surfaces Σ, Σ and Σ . In addition, the Euler relation gives

∂p(Q) ∂n = -jωρ 0 V n (Q) ∀Q ∈ Σ, (2)
∂p(Q ) ∂n = 0 ∀Q ∈ Σ , (3)
∂p(Q ) ∂n = -jωρ 0 V n (Q ) ∀Q ∈ Σ , (4)
where V n (Q) is the normal velocity at point Q, ω is the angular frequency and ρ 0 is density of air. In addition, if there is no source inside the volume Ω, then the pressure p(N ) is given by the Helmholtz equation

∆p(N ) + k * 2 p(N ) = 0 ∀N ∈ Ω. (5)
Considering Eqs. ( 1) to ( 5), one can write

- Ω (∆Φ(N ) + k * 2 Φ(N ))p(N )dN = + Σ p(Q) ∂Φ(Q) ∂n + jωρ 0 V n (Q)Φ(Q)dQ + Σ p(Q ) ∂Φ(Q ) ∂n + Σ p(Q ) ∂Φ(Q ) ∂n dQ + Σ jωρ 0 V n (Q )Φ(Q )dQ . (6)
The main interest of Green's identity is that the function Φ can be arbitrarily chosen provided that it is twice dierentiable. Let's choose Φ as a mode shape φ n (N ) of the equivalent nite volume Ω such as

∆φ n (N ) + k 2 n φ n (N ) = 0 ∀N ∈ Ω, (7)
where k n is the natural wavenumber of mode n of the virtual acoustic volume Ω. Mode shapes φ n (N ) constitute an orthonormal basis of functions on which the pressure at point N can be decomposed

p(N ) = ∞ n=0 a n φ n (N ) (8)
Introducing Eqs. ( 7) and (8) in Eq. ( 6) and using the orthonormal property of mode shapes, one obtains

p(N ) = ∞ n=1 φ n (N ) Λ n (k * 2 -k 2 n ) C n , (9)
where

C n = Σ p(Q) ∂φ n (Q) ∂n dQ +jωρ 0 Σ V n (Q)φ n (Q)dQ + Σ p(Q ) ∂φ n (Q ) ∂n dQ + Σ p(Q ) ∂φ n (Q ) ∂n dQ +jωρ 0 Σ V n (Q )φ n (Q )dQ . (10)
Eqs. ( 9) and (10) are valid whatever the boundary conditions of the "virtual" mode shapes φ n (N ). Depending on the chosen boundary conditions, some terms in Eq. (10) will be suppressed. Then, it has been here proven that the pressure at point N can be expressed using a "virtual" problem for which the boundary conditions can be arbitrarily chosen.

Acoustic radiation of an oil pan

A baed oil pan, shown in Fig. 3 is excited by a point force and radiates in an semi-innite acoustic environment.

Reference computation

A classical computation based on the use of innite elements has been done using ACTRAN software. An example of a pressure map is given in Fig. 4 Let's consider the modes shapes φ n (N ) with Neumann boundary conditions on surfaces Σ, Σ and Σ as

∂φ n (Q) ∂n = 0 ∀Q ∈ Σ ∪ Σ ∪ Σ (11)
In that case, C n simplies to

C n = jωρ 0 Σ V n (Q)φ n (Q)dQ +jωρ 0 Σ V n (Q )φ n (Q )dQ , (12)
and the pressure p(N ) at any point of the domain Ω is computed knowing the particle velocities on surfaces Σ and Σ only. One can see in Eq. ( 12) that the choice of a Neumann condition on surface Σ allows not to need to know information on this surface. This is actually due to the fact that the chosen boundary condition for the associate problem (the mode shapes) is the same as the real problem.

Green's function with mixed boundary condition

A second possibility is to consider the boundary condition on the surface Σ as a Dirichlet boundary condition. In that case,

∂φ n (Q) ∂n = 0 ∀Q ∈ Σ ∪ Σ (13)
and

φ n (Q ) = 0 ∀Q ∈ Σ . (14)
Thus, C n simplies to

C n = jωρ 0 Σ V n (Q)φ n (Q)dQ + Σ p(Q ) ∂φ n (Q ) ∂n dQ , (15)
and the pressure p(N ) at any point in the domain Ω is computed knowing the particle velocity on surface Σ and pressure on surface Σ . However, in this formulation, one has to know the spatial derivative of the mode shapes expressed in terms of pressure. These are proportional to the modes shapes expressed in terms of particle velocity as

∂φ n (Q) ∂n = -jω n ρ 0 χ n (Q), (16)
where χ n (Q) is the mode shapes of the virtual acoustic volume expressed in particle velocity. These mode shapes can be obtained with a numerical solver as ACTRAN.

Comparion of computations

The aim is here to verify if Eq. ( 12) (uniform conditions) and Eq. (15) (mixed conditions) permit to compute the pressure at point N by comparing the results to reference calculation presented in Fig. 4. For that purpose, the integrals are replaced by sums (discretisation of the surfaces) and the pressure and/or particle velocities on surfaces Σ and Σ are taken from reference calculation. The modes shapes φ n (N ) and its derivative ∂φn ∂n are computed using a numerical solver (ACTRAN) up to 5kHz on the virtual volume presented in Fig. 5. The pressures p(N ) obtained with the three computations are compared in Fig. 6 which demonstrates that Green's identity can be used to compute the pressure in a domain whatever the chosen boundary conditions. The slight discrepancies between curves are due to the discretisation of the problem and to the convergence of the modal sums. It is important to underline that neither the Neumann nor the Dirichlet boundary conditions are representative of the real conditions. This does not prevent the calculation. However, the convergence of the calculation seems to be quicker with Dirichlet condition.

Source reconstruction using

Green's identity

In both equations ( 12) and (15) the normal velocity of surface Σ appears. Let's consider now that this veloc- ity eld is unknown. If the pressure p(N ) is measured at several points, these equations can be inverted and the velocity eld can be identied. This is the principle of the method called inverse Patch Transfer Function [START_REF] Aucejo | Identication of source velocities on 3D structures in non-anechoic environments: Theoretical background and experimental validation of the inverse patch transfer functions method[END_REF][START_REF] Totaro | Sound elds separation and reconstruction of irregularly shaped sources[END_REF]. This method corresponds here to the inversion of the Green's identity with Neumann boundary conditions. In this paper, it will be called u-iPTF. The other possibility, also presented here, will be called m-iPTF.

Uniform BC: u-iPTF

If the vibrating surface Σ is discretized into N s elementary surfaces called patches and the "virtual" sur-face Σ is discretized into N m patches, then, Eqs. ( 9) and ( 12 where pi and Vj are the space averaged pressure and particle velocity on patches i and j of surface Σ and Vl is velocity of patch l on surface Σ. Z ij is the acoustic impedance between patches i and j. In a matrix form, one obtains

pi = Z ij Vj + Z il Vl (18)
The velocities space averaged on patches l on the vibrating surface Σ can then be deduced using

Vl = Z -1 il pi -Z ij Vj (19)
The velocity eld on the surface of the source can be obtained measuring pressure and particle velocity on a virtual surface Σ . These colocated measurements can be done using pU probe. 9) and (15), the pressure at one point N in the domain can be expressed as

p(N ) = Nm j=1 Y * N j pj + Ns l=1 Z * N l Vl (20)
If the pressure is measured at several points in the cavity, Eq. ( 20) can be written in a matrix form

p i = Y * ij pj + Z * il Vl
(21) Finally, the velocity eld on the source is obtained by

Vl = Z * -1 il p i -Y * ij pj (22)
Compared to Eq. ( 19), the velocity eld Vl is obtained measuring the pressure on the virtual surface Σ and at several points in the domain Ω. In addition, impedance Z * il and the admittance Y * ij are computed with mode shapes with mixed BC while Z il and Z ij are obtained with Neumann BC.

Application on the oil pan

In both cases, Eqs. ( 19) and ( 22), the problem is reduced to the classical form Ax = b. Both resolution are ill-posed and need regularization. In the following, the GCV method [START_REF] Leclère | Acoustic imaging using under-determined inverse approaches: Frequency limitations and optimal regularization[END_REF] will be used to nd the best regularization parameter of the classical Tikhonov regularization technique. Fig. 7 shows the mean square velocity of the vibrating surface Σ obtained with reference calculation and identication methods and Fig. 8 presents the comparison between the reference and the velocity elds identied by u-iPTF and m-iPTF at 700Hz. This application shows that these approaches, based on the use of the Green's identity coupled to FEM solver, can lead to powerful inverse methods. These methods are intrinsically independent of the acoustic environment and can handle sources with complex geometries.

Experimental validation

An experimental validation of the described procedure has been performed with both approaches. However, as the results using u-iPTF have already been exposed in [START_REF] Totaro | Sound elds separation and reconstruction of irregularly shaped sources[END_REF], only the results obtained using m-iPTF will be presented here. The oil pan was xed on a wall and excited by a shaker. It was radiating in a semi-anechoic chamber. Some reference measurements have been done using a pU probe near the oil pan (on a rectangular box close to the surface of the oil pan). An example of these measurements can be seen in Fig. 10 (left).

Then, to apply m-iPTF, the pressures have been measured on a virtual surface (a bigger rectangular box) and inside the virtual volume. The result obtained at 976Hz is shown in Fig. 10 (right). The comparison in Fig. 10 demonstrates that the source eld reconstruction is possible in real conditions. In addition, contrary to reachable measurements, the identication is computed on the real surface of the structure. Finally, the source is completely characterized as pressure and intensity elds can be deduced from identied velocity eld and the direct formulation, Eq. ( 19) or (22). As a demonstration, the acoustic power of the source, measured with pU probe or identied, is plotted in Fig. 11. Both curves agree. As a consequence, the acoustic power of the source is obtained measuring pressures only thanks to m-iPTF.

CONCLUSIONS

The Green's identity is a technique extensively used in acoustics to solve vibro-acoustic problems. In this paper, it is shown how the Green's identity and FEM solver can be used to develop source eld reconstruction methods. Two developments are presented here. The rst one, the u-iPTF, permits to reconstruct the velocity eld of a source with complex geometry measuring colocated pressure and particle velocity on a virtual surface surrounding the source. The second one, the m-iPTF, only needs pressure measurements on a virtual surface and inside the virtual volume. Both of them need regularization as the problem is ill-conditionned. However, both of them allow source elds reconstructions in non-anechoic room on complex geometry. In addition, the source is then fully described by velocity, pressure, intensity elds and acoustic power and radiation eciency.

Figure 1 .

 1 Figure 1. A baed structure radiates in an acoustic environment.

Figure 3 .

 3 Figure 3. Oil pan.

  at 2000Hz.

Figure 4 .

 4 Figure 4. Map of the pressure eld computed with innite elements at 2000Hz.

3. 2 .

 2 Green's function with Neumann boundary conditions

Figure 5 .

 5 Figure 5. Virtual volume modeled by Finite Elements.

Figure 6 .

 6 Figure 6. Pressure at point N as a function of frequency. Black solid line: reference computation using Innite Elements; Blue dashed-dotted line: Green's function with Neumann BC; Red dotted line: Green's function with mixed BC.

Figure 7 .

 7 Figure 7. Mean Square Velocity as function of frequency. Black solid line: Reference; Blue dashed-dotted line: u-iPTF; Red dotted line: m-iPTF

4. 2 .

 2 Mixed BC: m-iPTF Using Eqs. (

Figure 9 .

 9 Figure 9. Experimental setup. Left: oil pan and line of microphones; right: electrodynamic shaker.

Figure 10 .

 10 Figure 10. Amplitude of the velocity eld at 976Hz. Left: measured with pU probe; Right: identication using m-iPTF.

Figure 8 .

 8 Figure 8. Velocity eld of the vibrating surface Σ at 700Hz. Left: Reference; Center: u-iPTF; Right: m-iPTF

Acknowledgement

This work was performed within the framework of the Labex CeLyA of Université de Lyon, operated by the French National Research Agency (ANR-10-LABX-0060/ ANR-11-IDEX-0007).


S Forget

N Totaro

J-L Guyader

M Schaeffer

SOURCE FIELDS RECONSTRUCTION ON A 3D STRUCTURE IN NOISY ENVIRONMENT

Nowadays, noise source identification remains a current challenge in industry to be able to act at the origin of noise. Since decades, many identification methods have been developed but most of them are not suitable for industrial applications. In this context, the presented method named Mixed -inverse Patch Transfer Functions (M-iPTF), could be promising as it is able to reconstruct acoustic fields on irregularly shaped sources, only needs pressure measurements and can be performed in a noisy environment. The present paper describes the theoretical background of the M-iPTF method, based on the application of the Green's identity on a virtual cavity surrounding the source, and its application results on a baffled oil pan excited with an electro-dynamic shaker, with and without the presence of a secondary correlated source.

INTRODUCTION

Since early 1990s, many noise source identification methods have been developed and some have matured and have been used throughout many industries [START_REF] Lanslots | Selecting sound source localization techniques for industrial applications[END_REF]. Even though all these techniques have the same goal, that is to say to know from where the sound originates, it does not exist one method that prevails on other. A choice has to be made depending on the source object, the operating conditions and the available sensors and facilities. Usually, these methods are classified into three categories: Near-field Acoustic Holography (NAH) approaches, Beamforming techniques, and inverse methods.

The NAH, firstly introduced by Maynard et al [START_REF] Julian D Maynard | Nearfield acoustic holography: I. theory of generalized holography and the development of nah[END_REF], uses the 2D Fourier transform of the Helmholtz equation and pressure measurements from a microphones array to reconstruct the source fields. The main drawbacks are : (i) the back propagation can only be applied to simple geometries (plane, cylindrical or spherical [START_REF] Earl | Fourier acoustics: sound radiation and nearfield acoustical holography[END_REF]), (ii) measurements have to be done in the near field to prevent boundary effects and (iii) the source has to be placed in free field conditions.

The beamforming technique [START_REF] Billingsley | The acoustic telescope[END_REF] uses a microphones array placed in the far field and process suitable time delays on signals depending on the microphones position, to determine by summation the source location. However, it is only a qualitative approach and the method presents some limitations at low frequencies.

In the third category, a promising approach, the inverse Patch Transfer Functions (iPTF) developed in [START_REF] Mathieu Aucejo | Identification of source velocities on 3d structures in non-anechoic environments: Theoretical background and experimental validation of the inverse patch transfer functions method[END_REF][START_REF] Vigoureux | Inverse patch transfer functions method as a tool for source field identification[END_REF], could address the problems mentioned above. Indeed by defining a closed acoustic volume surrounding the source and by using the Green's identity, the iPTF method can separate the contribution of the source from what coming from its environment. In this sense, it is a separation method which allows to be applied in a non anechoic environement. Additionally, it uses the computation of eigen modes of this cavity by FEM to handle cavities and sources with complex shapes [START_REF] Totaro | Sound fields separation and reconstruction of irregularly shaped sources[END_REF]. The back propagation can thus be made directly on the source surface. The challenge is the experimental process. It requires to measure particle velocities and pressures on the virtual surface surrounding the cavity. These can be addressed by using a p-U probe but its quite expensive cost could be a barrier to an industrial development.

The proposed method relies on the same concept than the iPTF method, i.e the definition of a virtual cavity surrounding the source and the application of the Green's identity. The difference comes from the boundary conditions taken for this cavity. Indeed, as they are independent from the real boundary conditions of the problem, they can be chosen arbitrarily. For the classical iPTF method, named Homogeneous-iPTF (H-iPTF), all surfaces of the cavity satisfy the homogenous Neumann condition. The developed method, named Mixed-iPTF (M-iPTF), uses mixed boundary conditions for the virtual cavity so that the measured data change. Here, instead of measuring pressures and particle velocities, the M-iPTF method only requires pressure measurements. This enables to consider a possible industrial measurement process by using a simple microphones array. Obviously, the M-iPTF method keeps the source separation property, which does not require free field conditions and the back propagation can be made directly on a complex 3D source.

MIXED -INVERSE PATCH TRANSFER FUNCTIONS

Theoretical background

Let's consider the vibro-acoustic problem presented in Figure 1. It consists in a vibrating surface Σ (source denoted S) radiating in any acoustic environment. In the present example, the source S is located near a rigid wall Σ . To insulate the source S from its acoustic environment, which may contain secondary stationnary sources (denoted S in Figure 1) or objects, ones defines a virtual surface Σ , arbitrarily chosen, surrounding the source to identify. It is important to notice here that this virtual surface has no physical meaning.

A "virtual" closed acoustic volume (denoted Ω in Figure 1) is then delimited by surfaces Σ, Σ and Σ . The Green's identity can be used to compute the pressure at any point in the volume knowing quantities (pressure and/or particle velocity) on boundaries. In its general form, the Green's identity can be written as :

Ω [Φ∆Ψ -Ψ∆Φ] dΩ - Γ [Φ ∂Ψ ∂n -Ψ ∂Φ ∂n ] dΓ = 0 (1)
The functions Φ and Ψ of Equation ( 1) are replaced by p(M ) and φ n (M ) in Equation ( 2), which corresponds respectively to the pressure at any point M in the cavity and the mode shapes of this cavity expressed in pressure. This choice is arbitrary provided that both functions are twice differentiable. In that case the pressure and pressure eigenmodes of the cavity allow to simplify the left-hand side of the Green's identity in Equation ( 2), thanks to the use of the Helmholtz-Kirchhoff equations ( 3) and [START_REF] Billingsley | The acoustic telescope[END_REF].

Ω [φ n (M )∆p(M ) -p(M )∆φ n (M )] dΩ = Σ [φ n (Q) ∂p(Q) ∂n -p(Q) ∂φ n (Q) ∂n ] dΣ + Σ [φ n (Q ) ∂p(Q ) ∂n -p(Q ) ∂φ n (Q ) ∂n ] dΣ + Σ [φ n (Q ) ∂p(Q ) ∂n -p(Q ) ∂φ n (Q ) ∂n ] dΣ (2)
∆p(M ) + k * 2 p(M ) = 0 ∀M ∈ Ω (3)
where, k * = ω/c * = k(1 + iη) -1 is the complex wave number, η the damping loss factor, ω the angular frequency and c the speed of sound in acoustic medium; and

∆φ n (M ) + k 2 n φ n (M ) = 0 ∀M ∈ Ω (4)
where, k n is the wave number at the eigen frequency ω n .

Then, to simplify the right-hand side of the Green's identity (2), the Euler's condition is applied :

∂p(Q) ∂n = -iωρ 0 V N (Q) ∀Q ∈ Σ (5) ∂p(Q ) ∂n = -iωρ 0 V N (Q ) ∀Q ∈ Σ (6) ∂p(Q ) ∂n = 0 ∀Q ∈ Σ (7)
where, V N is the outgoing normal velocity of the concerned surface.

Finally, boundary conditions of the virtual cavity have to be defined. Again, this choice has no physical meaning, the boundary conditions are defined by convenience as the Green's identity allows it. The classical H-iPTF method relies on a rigid wall cavity satisfying the homogeneous Neumann's condition. For the presented method, the source surface Σ and rigid surface Σ fulfill the homogeneous Neumann's condition [START_REF] Hadamard | Sur les problèmes aux dérivées partielles et leur signification physique[END_REF], while the virtual surface Σ satisfies the Dirichlet's condition [START_REF] Nikolaevich | Solutions of ill-posed problems[END_REF].

∂φ n (Q) ∂n = 0 ∀Q ∈ Σ, Σ (8)
φ n (Q ) = 0 ∀Q ∈ Σ (9)
Then, by decomposing the pressure p(M ) on a basis composed of the cavity pressure eigenmodes φ n (M ) and by using their orthogonality property, the pressure at any point M in the cavity can be expressed by :

p(M ) = ∞ n=1 iρ 0 φ n (M ) [k * 2 -k 2 n ] Λ n Σ ω φ n (Q) V N (Q) dΣ - Σ ω n χ n (Q ) p(Q ) dΣ (10)
where, Λ n represents the norm of the n-th mode, χ n are the mode shapes of the cavity expressed in velocity given by the Euler's equation ( 11)

∂φ n (Q ) ∂n = -iω n ρ 0 χ n (Q ) (11)
Finally, to be computed, Equation [START_REF] Per | Rank-deficient and discrete ill-posed problems: numerical aspects of linear inversion[END_REF] is written in a discrete form. All surfaces are divided into elementary surfaces named patches and the modal basis of the cavity is truncated to a mode N for which the convergence have to be ensured

p(M ) = N n=1 iρ 0 φ n (M ) [k * 2 -k 2 n ] Λ n   N j j=1 ω A j φ n j V N j - N k k=1 ω n A k χ n k p k   (12)
where, < X > S represents the space average of the variable X on surface S, and A S is the area of the surface S. Subscripts j and k denote respectively a patch on the source surface Σ and a patch on the virtual surface Σ .

In Equation [START_REF] Leclère | Acoustic imaging using under-determined inverse approaches: Frequency limitations and optimal regularization[END_REF], the physically rigid surface no longer appears in the calculation. It has been already taken into account through the rigid boundary condition of the virtual cavity.

Considering the pressure computation of several points M in the cavity, Equation ( 12) can be written in a matrix form :

p i = Z ij V j -Y ik p k ( 13
)
where

Z ij = N n=1 iωρ 0 (φ n ) i φ n j A j [k * 2 -k 2 n ] Λ n and Y ik = N n=1 iω n ρ 0 (φ n ) i χ n k A k [k * 2 -k 2 n ] Λ n (14)
subscript i denotes a measured point inside the cavity.

Equation (13) expresses the direct problem formulation with the computation of the radiation of the source S and its environment. Conversely, to retrieve the field coming from the source, the problem requires to be inverted :

V j = Z -1 ij (p i + Y ik p k ) (15)
Thus finally, by using Equation (15), it is possible to reconstruct the source velocity field from only pressure measurements performed on the virtual surface and in the cavity. These measurements can be done using a microphones array, allowing a cost-effective industrial measurement process. The use of the Green's identity has allowed to separate the contribution of the source from what coming from its environment, making the M-iPTF method usable in a non-anechoic environment. At last, the eigen-modes of the cavity expressed in pressure and in velocity can be computed using a finite elements solver, enabling the cavity and thus also the source to be of complex shapes .

Once the velocity field is known, the source is fully determined because all others acoustic quantities can be retrieved: the boundary pressure field can be computed using the direct problem formulation (16) and the normal intensity field can be obtained by multiplying the pressure and velocity fields (17). In parallel, other indicators like the radiated acoustic power or the radiation efficiency can be deduced. The source is then fully described by the proposed inverse method.

p j = Z jj V j -Y jk p k (16)
I j = 1 2 (p j • V * j ) (17)
where, (•) and • * represent the real part and the complex conjugate of a complex number. • is the Hadamard's product.

Tikhonov's regularization

The direct problem formulation of the M-iPTF method can be written on the general form of a linear system Ax = b :

Z ij V j = pi with pi = p i + Y ik p k (18)
In acoustics and more generally in all domains of physics, inverse problems are often ill-posed problems in the Hadamard's sens [START_REF] Hadamard | Sur les problèmes aux dérivées partielles et leur signification physique[END_REF][START_REF] Nikolaevich | Solutions of ill-posed problems[END_REF]. According to him, an inverse problem is well-posed when these three conditions are satisfied together : the solution exists, is unique and stable.

In general, the first condition will be satisfied, the M-iPTF method tries to identify an input from an observed output.

The second condition is relating to the problem determination, i.e the ratio between the number of measured points in the cavity and the number of source patches, and the rank of the system governed by the number of modes taken into account. In general, in order to limit the number of measurements, the system is under-determined. Thus, under these conditions, the solution is not unique and the Z ij matrix is not square [START_REF] Per | Rank-deficient and discrete ill-posed problems: numerical aspects of linear inversion[END_REF].

Finally, the solution have to be stable. Many sources of instabilities exist, in particular coming from the measurement process, which can disturb the identification. Indeed the pi vector is obtained by measurements and thus can contain errors and measurement noise. Also, the size of patches can limit the maximal frequency of identification and the modal truncation can affect the convergence of the solution.

Thus, due to the discretization and the measurement process, the inverse problem of the M-iPTF method is ill-posed and requires a reguralization. Many regularization methods exist in the literature [START_REF] Kim | Optimal regularisation for acoustic source reconstruction by inverse methods[END_REF][START_REF] Leclère | Acoustic imaging using under-determined inverse approaches: Frequency limitations and optimal regularization[END_REF], but only the so-called Tikhonov regularization has been used in this paper. The Tikhonov-regularized solution of Equation ( 18) is given by [START_REF] Totaro | Sound fields separation and reconstruction of irregularly shaped sources[END_REF] :

V β j = rank(Z ij ) m=1 v m u H m σ m σ 2 m σ 2 m + β 2 pi ( 19
)
where σ m , u m and v m are respectively the singular values in the decreasing order, the left and right singular vector of the Z ij matrix coming from its Singular Values Decomposition. • H denotes the Hermitian transpose.

β is a real positive parameter that prevents small singular values, associated with small data values and thus dominated by noise, to amplify the norm of the regularised solution during the inversion process. It is a tradeoff between measurement error and error of regularization. The optimal β parameter can be determined by the two well-known methods, the so-called L-Curv (LCV) and Generalized Cross Validation (GCV) methods [START_REF] Kim | Optimal regularisation for acoustic source reconstruction by inverse methods[END_REF][START_REF] Leclère | Acoustic imaging using under-determined inverse approaches: Frequency limitations and optimal regularization[END_REF]. These regularization methods were both applied to the M-iPTF method but only the LCV regularisation is presented thereafter due to similar results in between.

EXPERIMENTAL VALIDATION

System under study

A validation of the M-iPTF method was carried out on a real baffled oil pan (Figure 2(a)) excited by an electro-dynamic shaker (Figure 2(b)).

The measurement process was performed using a three axis robot and a 6x2 regular microphones array. The virtual surface was taken rectangular to simplify measurements (Figure 3(b)) and the FE modelisation of the cavity for its modal extraction (Figure 3(c)). As it can be seen in Equation (18), one needs the measurement of space averaged pressure p k on patches of the virtual surface and pressure p i at points inside the volume. In experiments, the space averaged pressure is approximated by the value at centers of patches. Also, points in the cavity were taken randomly in the cavity. Thus, combining pressure measurements and a numerical model, Equation (19) permits to reconstruct the velocity, pressure and intensity fields on the surface of the structure (Figure 3(a)). In parallel, to validate the present approach, some reference measurements have been done, by using a p-U probe, on the nearest rectangular box surrounding the oil pan (the measurement of pressure and particle velocity on the real surface of the oil pan was not possible). As it was mentioned in section 2.2, several factors can influence the identification results. A parametric study was carried out to assess their influence but it will be presented in further works. For the present paper, these parameters (Table 1) are taken optimal to ensure a proper identification and to illustrate the proposed method.

The ability of the M-iPTF method to correctly reconstruct the oil pan fields was tested with and without the presence of a secondary correlated source. Results are presented in sections 3.2 and 3.3. Table 1. Information concerning the identification process.

Identification results without a secondary correlated source

Figure 4 shows the comparison between the maps of the reference velocity and pressure fields measured with a p-U probe and those obtained from the computation of the M-iPTF method with the LCV -Tikhonov regularization at the frequencies of 554, 1004 and 1951 Hz. One has to remind that the "reference" measurements have been done on a rectangular box close to the surface of the oil pan. Considering this remark, Figure 4 demonstrates a very accurate reconstruction of the velocity field as well as the pressure field. The slight shift comes from the fact that the reference surface is slightly larger than the oil pan surface on which the identification is done. In Figure 5, the mean square velocity, pressure and the radiated acoustic power are plotted. p-U probe measurements are compared to the fields identified using the LCV regularization. Below 300 Hz, due to a low signal-to-noise ratio, the LCV strategy of regularization seems to be inappropriate in that case. Above 300 Hz, results are very good even if a slight shift for the amplitude of the velocity can be observed. This quantity is highly sensitive to the difference between the "reference" surface and the real surface of the oil pan used for identification. On the contrary, as expected, the acoustic radiated power is conservative and not affected by this difference between surfaces as it can be seen in Figure 5(c).

for industrial applications. The method relies on the application of the Green's identity on a virtual cavity defined around the source and the computation by FEM of the eigen-modes of this cavity. That permits to the M-iPTF method to be a separation method, i.e to be performed in a non anechoic environment, and the source to be of complex shape. Likewise, the mixed chosen boundary conditions for the cavity allow to use only pressure measurements. However, as inverse problems, the method suffers from ill-conditioning and have to be regularized. The M-iPTF method has shown its ability to correctly reconstruct all acoustic source fields in a large frequency range in an experimental application with a complex 3D source and a noisy environment.

Figure 1 :

 1 Figure 1: Definition of the virtual closed volume Ω and the boundary surfaces Σ (surface of the vibrating source), Σ (virtual surface surrounding the source) and Σ (physically rigid wall).

  (a) Source : baffled oil pan (b) Position of the shaker

Figure 2 .

 2 Figure 2. System under study : (a) a baffled oil pan excited by (b) an electro-dynamic shaker.

Figure 3 :

 3 Figure 3: (a) Patches model for the source, (b) for the virtual surface, (c) FE model of the cavity and (d) surface of the reference measurements.

Figure 4 :

 4 Figure 4: Comparison between p-U probe measurements (velocities ((a)-(e)-(i)) and pressures ((c)-(g)-(k)) and identification results with LCV regularization (velocities ((b)-(f)-(j)) and pressures ((d)-(h)-(l))) respectively at 554 Hz, 1004 Hz and 1951 Hz.

ACKNOWLEDGMENTS

This work was performed within the framework of Labex CeLyA of Université de Lyon, operated by the French National Research Agency (ANR-10-LABX-0060/ ANR-11-IDEX-0007).

Identification results in a noisy environment

The oil pan environment is now considered as noisy by the presence of a secondary correlated source placed outside the cavity domain (Figure 6). In Figure 7, the mean square velocity and pressure measured using the p-U probe are plotted in three different configurations : the oil pan radiates only, the secondary source is alone and both are activated. Here, depending on the frequency range, one source can prevail on the second. Thus, to validate the method in a noisy environment, only the configuration for which both sources generate an equivalent sound pressure level, i.e at the frequencies of 1713 Hz and 2080 Hz in Figure 7 In Figure 8, the maps of the velocity and pressure fields measured by using the p-U probe are compared to those obtained from the computation of the M-iPTF method with the LCV -Tikhonov regularization at the frequencies of 1713 Hz and 2080 Hz. At these frequencies, for which the sound generated by the oil pan and the secondary source is equivalent, the reconstruction of the velocity and pressure fields is quite satisfactory. The slight difference of amplitude, lower for the identified fields, might come from an over regularization due to the important noise coming from the correlated source.

CONCLUDING REMARKS

In acoustics, many identification methods exist, due to the variety of applications. The proposed method in this paper belongs to the inverse methods domain and has been developed in particular