Julien Nussbaum

Philippe Helluy

Jean-Marc Hérard

Alain Carrière

Jean-Marc H Érard

Alain Carri Ère

Numerical simulations of gas-particle flows with combustion

Keywords: finite volume, non-conservative scheme, two-phase flow, solid propellant, ignition, internal ballistic

come

Internal ballistic problem

The two-phase flows in guns are very difficult to model. Many inter phase interaction and complex phenomena occur while the powder burns and the bullet moves in the gun tube.

In order to simplify the physical geometry of a gun, we consider a cylindrical combustion chamber linked to a tube of the same constant cross section. At the initial time, the mixture of gas-powder grains is contained in the combustion chamber, limited by the breech at one end, and by the shot base at the other end. The initial geometry is illustrated on Figure 1.

Propellant is ignited by a hot gas stream from an igniter. The igniter is a cylinder filled with black powder, perforated by several holes. Hot combustion gases escape through the holes. The addition of energy increases the propellant grain surface temperature. The combustion occurs when the ignition temperature is reached. After ignition, the solid propellant burns and gases are produced. We assume that the combustion products are similar to the initial gas species. The pressure increases in the combustion chamber, while the front flame propagates in the powder bed according to the following mechanism: ignition of some grains produces a hot gas stream (as igniter does) that locally increases the temperature and pressure. The gases propagate in the domain, and by local heat transfer from the gas phase to the solid phase, the other grains are ignited. The bullet begins to move when the pressure at the shot base is greater than the start pressure, and goes into the tube until the muzzle. Internal ballistic studies stop when the bullet exits from the tube.

Experimental studies exist but are limited: pressure measurements are only possible at the breech or at the shot base and the only available velocity is the bullet's velocity at the muzzle. The temperature and pressure gradients, the high velocities of the flow and the particles volume fractions are still difficult to measure.

The mathematical model

One can find many two-phase flow models in the scientific literature, focusing on different approaches, with different advantages and drawbacks. We refer for example to the book of Gidaspow [START_REF] Ère | Gidaspow Multiphase flow and fluidization[END_REF] for a survey of such models. We will concentrate here on the Gough's model [START_REF] Gough | Modeling of two-phase flows in guns[END_REF] that is very popular in the ballistic community. Our work is devoted to find a numerical method in order to replace a Mac-Cormack integrator nuss_hell.tex; 10/02/2006; 11:02; p.2 that is not able to simulate non-classical weapons because of numerical instabilities. For details on the evolution of ballistic models, see [START_REF] Carrière | Sur l'évolution des modèles de Balistique Intérieure de Paul Vieille à nos jours, et leur validation sur diverses armes, des pistolets au canon de Jules Verne[END_REF].

The system is made of the mass and momentum conservation laws for both phases, energy conservation for the gas phase and two convection equations for the enthalpy of the grains and the thickness of burnt powder. The flow is supposed to be one-dimensional in the direction of the symmetry axis of the gun tube with constant section.

In a first stage after ignition, the shot base does not move. In this case, the governing set of equations takes the form

∂ W ∂ t + ∂ F (W ) ∂ x + C (W ) ∂ W ∂ x = S (W ) , (1)
The vector W = W (x, t) ∈ Ω ⊂ R 7 is the unknown vector. The time variable is noted t ∈]0, T [ and the space variable x ∈]0, L[, where L is the length of the tube. The flux vector F and the sources vector S are functions from ]0,

L[×]0, T [ to R 7 and C is a function from ]0, L[×]0, T [ to R 7×7 . We define W =            α 2 α 1 ρ 1 α 1 ρ 1 u 1 α 2 ρ 2 u 2 α 1 E 1 H ts d            , F (W ) =            α 2 u 2 α 1 ρ 1 u 1 α 1 ρ 1 u 2 1 + p 1 α 2 ρ 2 u 2 2 + p 2 α 1 u 1 (E 1 + p 1 ) 0 0            , C (W ) ∂W ∂x =            0 0 -p 1 ∂ x α 1 p 1 ∂ x α 1 p 1 ∂ x (α 2 u 2 ) u 2 ∂ x H ts u 2 ∂ x d            (2)
where α k is the volume fraction, ρ k the density, u k the velocity, p k the pressure and E k the total energy of phase k. The index k = 1 corresponds to the gas phase and k = 2 corresponds to the solid phase. The solid phase is supposed to be incompressible. The specific enthalpy of the grains is noted H ts and the thickness of the burnt powder of each grain is noted d.

REMARK 1.

In the sequel, we will give a relation between p 2 and p 1 (see formula [START_REF] Krier | Flame spreading and combustion in packed beds of propellant grains[END_REF]). It will then be possible to eliminate the pressure p 2 in the system (1)-(2). Thus, the Gough's system falls into the category of the so-called two-velocity one-pressure models for two-phase flows.

REMARK 2. In our application, the pressure p 1 in the solutions of ( 1)-(2) appears to be smooth: the apparition of shock waves would indicate a bad functioning of the gun. Thus, it is not necessary to give a more precise definition of the non-conservative products in [START_REF] Gough | Modeling of two-phase flows in guns[END_REF].

Let us define the source terms by

S (W ) =             -Γc ρ 2 Γ c + Γ ign Γ c u 2 -D + Γ ign u ign -Γ c u 2 + D Γ c Q ex + p 1 ρ 2 + u 2 2 2 -u 2 D -A s q t + Q ign Γ ign κ q t ṙ            
.

In this source term:

-Γ c is the mass transfer rate, due to combustion, from the solid phase to the gas phase;

-Γ ign is the mass addition rate from the igniter;

-u ign is the gas velocity from the igniter;

-D is the interphase drag force;

-Q ex is the exothermic energy released by the solid phase during the combustion;

-Q ign is the energy released by the igniter;

-ṙ is the combustion rate;

-q t is the heat flux per specific surface unit between the two phases;

-A s = α 2 S p
Vp is the specific surface of the solid phase; -S p and V p are respectively the instantaneous surface and volume of a propellant grain;

-κ corresponds to the thermal diffusivity of the solid phase.

REMARK 3. The thermal diffusion and radiation are taken into account only at the local level: in one grain, and between one grain and the surrounding gas. It means that our model will not be able to predict accurately combustion fronts that are driven by thermal diffusion or radiation. The combustion front will propagate because of the convection of hot gases. Pressure waves can also trigger the combustion by a local increase of the temperature.

In a second stage, when the pressure in the combustion chamber is greater than the resistive pressure of the bullet, the shot base begins to move.

With a moving boundary, the mesh must be adapted to the expansion of the computing domain. We use a rezoning technique that consists in a change of variables in order to obtain a virtual static computation domain. We introduce the change of variable

ξ = x x p ,
where x is the real position and x p is the position of the shot base. Thus, 0 ≤ ξ ≤ 1, and we obtain a new expression of the system (1)-( 2)

∂ ∂ t (x p W ) + ∂ ∂ ξ (F (W ) -v p ξ W ) + C (W ) ∂ W ∂ ξ = x p S (W ) (3)
In order to compute the position of the shot base x p , we apply the fundamental principle of dynamics

m p dv p dt = A (p m -p r ) , (4)
where v p = dxp dt is the bullet speed, m p its mass, A the tube section, p m the pressure at the shot base and p r the resistive pressure (induced by the bullet/tube friction, and depending on the geometry of the system). The resistive pressure is supposed to be constant.

We consider the following wall boundary conditions:

-at the breech u 1 = u 2 = 0;

-at the shot base

u 1 = u 2 = v p .
At the initial time, the two phases are supposed to occupy homogeneously the initial volume (the combustion chamber).

Constitutive laws

The definition of α k as volume fraction (also called porosity) for each phase gives

α 1 + α 2 = 1 . (5)
The perfect gas equation of state is inappropriate in internal ballistic problems. High temperatures and pressures require using real gas equation of state, such as Noble and Abel law that reads

p 1 (ρ 1 , e 1 ) = (γ -1) ρ 1 e 1 1 -η ρ 1 , ( 6
)
where γ is the specific heat ratio and η is the gas covolume. The internal specific energy e 1 of the gas is given by

E 1 = ρ 1 e 1 + u 2 1 2 . ( 7
)
p 2 is defined by

p 2 (α 1 , ρ 1 , e 1 ) = p 1 (ρ 1 , e 1 ) + R p (α 1 , ρ 2 ) , (8)
R p is the intra granular stress, given in [START_REF] Porterie | Modélisation de la phase d'allumage d'une charge propulsive en balistique intérieure[END_REF],

R p = 0 if α 1 > α c ρ 2 c 2 p αc (αc-α 1 ) α 1 (1-α 1 ) if α 1 ≤ α c (9)
The critical porosity α c is usually taken equal to the initial porosity. The velocity c p is a measured sound speed in a particle bed far from compaction. For our configuration, it is sufficient to suppose that it is constant.

From the intra granular stress expression [START_REF] Hirschefelder | Molecular theory of gases and liquids[END_REF], it is possible to compute the speed of propagation of infinitesimal granular disturbances. It is noted a and in the sequel we shall call it the sound speed of the powder. Of course, when the porosity α 1 is far from 0, it is much smaller than the real sound speed of the incompressible compacted powder. We have

a 2 = R p ρ 2 - α 2 R p,α 1 ρ 2 . (10)
In our configuration, the velocity a is given by

a = 0 if α 1 > α c αc α 1 c p if α 1 ≤ α c . ( 11
)
We observe that the intra granular stress is zero when the porosity (the volume fraction of gas) is sufficiently high. On the other hand, it tends to infinity when the porosity tends to zero. This latter case corresponds to the compaction of the powder. It is clear that the compaction will produce a very high sound speed in the powder. In numerical simulations this can lead to a very constraining CFL condition. Fortunately, in our application, it appears that the combustion of the powder tends to increase the porosity during the computation and that the dominant sound speed is the sound speed of the gas. Hyperbolicity of the system (1)-( 2) depends directly on the intra granular stress. A sufficient strength implies hyperbolicity of the system. The eigenvalues of the Gough's model are computed in Appendix B. For more details, we refer to [5], [START_REF] Combe | Finite volume algorithm to compute dense compressible gas-solid flows[END_REF].

Theoretical expressions have been developed to model the drag force of a single particle in a gas flow, but the more complex case of a particle bed has been modelled only from correlations. We use a limit of Ergun's correlation given in [START_REF] Ergun | Fluid flow through packed columns[END_REF]

D = f r ϕ (α 2 ) 6 ρ 2 (1 -α 2 ) S p V p (u 1 -u 2 ) |u 1 -u 2 | , (12)
with f r the resistive factor. The function ϕ depends on the shape of the grains. Here we use

ϕ (α 2 ) =      0.3 α 2 > 0.9 , 1.75 1-α 2 α 2 αc 1-αc 0.45 α c < α 2 < 0.9 , 1.75 α 2 < α c . (13)
At last, we use the combustion law of Vieille [START_REF] Vieille | Étude sur le mode de combustion des substances explosives[END_REF] ṙ = a r P n + b ,

where a r , b and n are experimentally determined constants.

We describe the mass transfer rate by

Γ c = A s ρ 2 ṙ = (1 -α 1 ) Sp Vp ρ 2 ṙ , (15)
where S p and V p are respectively the instantaneous surface and volume of a grain powder, computed by geometric formulas depending on the shape of the grains. Possible expressions for S p and V p are given below.

The heat flux q t per unit of surface depends on the gas temperature T and particle surface temperature T ps and is defined as follows

q t = h t (T -T ps ) , (16)
with h t the total thermal transfer coefficient, sum of convective and radiative coefficients

h t = h c + h r . (17)
The radiative coefficient is computed by

h r = ε p σ (T + T ps ) T 2 + T 2 ps , (18)
with the hypothesis that the gas emissivity is equal to unity. ε p corresponds to the particles emissivity and σ to the Stephan-Boltzmann constant.

The convective coefficient is obtained from the Nusselt number deduced from correlations (see [START_REF] Krier | Flame spreading and combustion in packed beds of propellant grains[END_REF] for example). It reads

Nu = 6 h c k V p S p , (19)
where the Eucken approximation [START_REF] Hirschefelder | Molecular theory of gases and liquids[END_REF] for polyatomic gas gives the thermal conductivity of the gas k using the viscosity coefficient µ, the universal gas constant R and the specific heat at constant volume

c v k = 15 4 R µ 4 15
c v R + 3 5 . ( 20
)
We follow Porterie [START_REF] Porterie | Modélisation de la phase d'allumage d'une charge propulsive en balistique intérieure[END_REF] in order to find another expression of the Nusselt number Nu = 2 + 0.4 Re

2/3 p Pr 1/3 , ( 21
)
where Pr is the Prandlt number defined by Pr = µ c p /k with c p the specific heat of the fluid at constant pressure. From (20), we deduce the expression of the Prandlt number for polyatomic gas

Pr = 4 γ 9 γ -5 , (22)
γ corresponding to the specific heat ratio.

The expression of the surface temperature T ps is given after the integration of Fourier's law and by supposing a parabolic temperature profile in a spherical grain. See the Appendix in [START_REF] Porterie | Modélisation de la phase d'allumage d'une charge propulsive en balistique intérieure[END_REF] for more details. We compute T ps from

T ps = T ps 0 -3 2 h t Hts k 2 p + T ps 0 -3 2 ht Hts k 2 p 2 + 3 ht Hts T k 2 p -T 2 ps 0 1/2 (23)
where k p is the thermal conductivity of particles and T ps 0 is the initial surface temperature.

Numerical method

We use a version of the Rusanov scheme, a Godunov scheme based on an approximate Riemann solver. In general, Finite Volume Schemes are used for conservation laws. In our case, non-conservative terms induce an adaptation of the scheme. We follow the idea presented in [START_REF] Gallouët | Numerical modelling of two-phase flows using the two-fluid two-pressure approach[END_REF] for a two-fluid two-pressure model.

In order to approximate the solution, consider a space step h, a time step τ , the points ξ i = x i /x p = ih/x p and the instants t n = nτ . The computations cells are

C i =]ξ i-1/2 , ξ i+1/2 [.
The solution W = x p W of ( 3) is approximated in each cell C i and at each time t n by a constant vector

W n i = x n p W n i ≃ W (ξ, t n ) , ξ ∈ C i . (24)
The non-conservative finite volume scheme reads

h W n+1 i -W n i + τ F n i+1/2 -F n i-1/2 + τ G n i+1/2,--G n i-1/2,+ = τ S n i , (25)
where S n i = x n p S n i . We define the numerical conservative flux by the classical Rusanov flux

F n i+1/2 = 1 2 F n i+1 + F n i - s n i+1/2 2 W n i+1 -W n i , (26)
with

F n i = F n i -v n p ξ i W n i , (27)
and the numerical non-conservative fluxes by

G n i+1/2,-= C(W n i ) W n i+1 + W n i 2 , G n i-1/2,+ = C(W n i ) W n i + W n i-1 2 , (28)
The velocity s i+1/2 is the maximal wave speed at the interface i + 1/2. It is defined by

s n i+1/2 = max s n i , s n i+1 , s n i = max (|(u 1 ) n i | + c n i , |(u 2 ) n i | + a n i ) . (29)
In practice, we only have to consider the wave speed in the gas phase: in our configuration c n i is always greater than the granular wave speed a n i and |u 2 -u 1 | ≪ c. In order to satisfy the CFL stability condition we take the time step as τ = δ h max i,n

s n i (30)
The CFL number δ has to be < 1. In practice, we observe a stability of the scheme when this condition is satisfied. It indicates that the source terms are not stiff compared with the convective terms.

Let us notice that the chosen non-conservative numerical flux (28) is arbitrary. Other choices are possible. They will give the same result when the solution pressure is smooth. On the other hand, to compute discontinuous pressure solutions, it is important to define precisely the non-conservative products and the associated numerical approximations. For more details on these topics we refer for example to the work of Sainsaulieu [START_REF] Sainsaulieu | Traveling-wave solutions of convection-diffusion systems in nonconservative form[END_REF] and the cited references. Recently, relaxation two-velocity two-pressure models have been designed where the nonconservative products are associated to linearly degenerated waves and are thus naturally defined [START_REF] Gallouët | Numerical modelling of two-phase flows using the two-fluid two-pressure approach[END_REF]. It may be an elegant way to circumvent this difficulty.

The Rusanov scheme is known to be very robust but also very dissipative. This robustness is interesting here because a part of the computation occurs in a non-hyperbolic regime. The hyperbolicity of the model is discussed in Appendix B. For more precise computations it is necessary to use a higher order scheme and to modify the model in such a way that it is always hyperbolic. REMARK 4. Recently, relaxation models for two-velocity two-pressure flows have been designed that are always hyperbolic [START_REF] Gallouët | Numerical modelling of two-phase flows using the two-fluid two-pressure approach[END_REF], [START_REF] Saurel | A multiphase Godunov method for compressible multifluid and multiphase flows[END_REF]. It is possible, for some one-pressure models, to provide an approximation by a two-pressure model. In this approach it is necessary to fix a time scale τ r at which the two pressures equilibrate. When τ r → 0 we recover the one-pressure model. It would be interesting to construct such a twopressure model for approximating the Gough's system (but it is not an immediate application of [START_REF] Gallouët | Numerical modelling of two-phase flows using the two-fluid two-pressure approach[END_REF] or [START_REF] Saurel | A multiphase Godunov method for compressible multifluid and multiphase flows[END_REF]).

However, numerical experiments have shown that this approach, with τ r = 0, does not permit to eliminate the non-hyperbolic instability on fine meshes [START_REF] Herard | A relaxation method to compute two-fluid models[END_REF]. To our knowledge, the construction of a general onepressure hyperbolic model for gas-particle flows is still not achieved.

Numerical results

Academic validation

In order to validate our code, we first simulate a flow with constant pressure and constant velocity. It is a numerical test (no physical meaning) where we consider that no interaction occurs between the two phases and we set the intra granular stress R p = 0.

The initial condition is made of two different constant states (Riemann problem) in a 1 meter length computation domain.

The solution can be computed explicitly (with a source S(W ) = 0). The velocity and pressure remain constant, and the other quantities are simply convected at the constant velocity. The proposed numerical scheme ( 26), (28) has the important property that it preserves exactly the constant pressure and velocity fields. Despite its simplicity, this property is not satisfied by all the classical finite volume schemes.

The choice of the pressure law ( 6) is also important (see [START_REF] Gallouët | A hybrid scheme to compute contact discontinuities in one-dimensional Euler systems[END_REF]). This property also appears to be important to compute more complicated configurations.

Figure 2 displays the results obtained with a 1000 cells mesh at a final time t = 3 ms (the constant velocity and pressure are not represented).

We notice the high diffusion of the density and porosity by the Rusanov scheme on Figure 2. The method is only first order in time and space. It remains interesting for 1D computations, because it is always possible to refine the mesh.

Virtual 132 mm gun

Secondary, the simulation of a ballistic cycle in a virtual 132 mm gun (see [START_REF]Fluid dynamic aspect of interior ballistics[END_REF] and Appendix A) is computed. Only the expected bullet velocity at the muzzle, the maximal pressure at the breech and shot base and shot exit time are available. We use a 100 cells mesh, with the empirical values γ = 1.27, η = 1.0838 • 10 -3 and α c = 0.4225. The powder used is a 7-holes type and the corresponding geometrical functions are

S p = π (L 0 -2d) [D 0 -2d + 7 (d 0 + 2d)]+ π 2 (D 0 -2d) 2 -7 (d 0 + 2d) 2 , V p = π 4 (L 0 -2d) (D 0 -2d) 2 -7 (d 0 + 2d) 2 ,
where L 0 , D 0 and d 0 are respectively the length, the external diameter and the perforation diameter before the combustion. The evolutions of the breech and shot base pressures on time are displayed on Figure 3. Different simulations with other internal ballistic codes have given realistic ranges for the computed values. Table II summarizes the results and shows that our results are in agreement with the expected ones.

Real 60 mm gun

At last, simulations on a 60 mm gun are compared to experimental measurements given in [START_REF] Kay | Étude de l'influence des charges propulsives à haute densité de chargement[END_REF]. We compare the maximal breech pressure and the muzzle velocity coming from the experiments and the simulations, by using three different initial powder masses: 1.17 kg, 1.365 kg and 1.4625 kg. In Table III and IV the data are compared. We observe a good agreement between the simulations and the experiments even if the precision of the model decreases with the mass of powder.

Conclusion

In this paper, we have adapted the Gough's model to describe the two-phase flow with combustion in a gun. This is a non-conservative, two-velocity, one-pressure model that involves complicated source terms and lacks of hyperbolicity. We have proposed a numerical method, based on the Rusanov scheme, to simulate this kind of flow. We have been able to obtain satisfactorily on academic and real configurations.

We have now to improve the model in several directions:

-It is important, for practical configurations, to take into account combustion chambers with variable section;

-It would be interesting to find a relaxation two-velocity two pressure model to approximate the Gough's model. In this way we could avoid the lack of hyperbolicity and then envisage higher order numerical schemes. With an appropriate two-velocity two pressure model we also hope to be able to give a simple definition of the non-conservative products (even if it is not crucial for interior ballistics);

-Finally, a huge work has still to be done on the modelling of the combustion process: chemical model, combustion of grains with special shapes, dependence of combustion front velocity with the convection and the diffusion, etc.

B. Hyperbolicity domain of the model

We observe first that the hyperbolicity of the model is not affected by the rezoning: it is sufficient to study the eigenvalues of ( 1)-( 2). The last two equations in (1)-( 2) are convection equations associated to the eigenvalue λ = u 2 . Thus, we concentrate on the first five equations. We rewrite the system in a new set of variables Y = (α 2 , ρ 1 , u 1 , u 2 , e 1 ). In this way, without the source terms, the system becomes

Y t + B(Y )Y x = 0 , (31)
with

B(Y ) =         u 2 0 0 α 2 0 ρ 1 α 1 (u 2 -u 1 ) u 1 ρ 1 ρ 1 α 1 α 2 0 0 p 1,ρ 1 ρ 1 u 1 0 p 1,e 1 ρ 1 Rp α 2 ρ 2 - Rp,α 1 ρ 2 p 1,ρ 1 ρ 2 0 u 2 p 1,e 1 ρ 2 p 1 (u 2 -u 1 ) α 1 ρ 1 0 p 1 ρ 1 p 1 ρ 1 α α 1 u 1         . (32)
In order to simplify it further, we consider the gas entropy s(ρ 1 , e 1 ) satisfying

T 1 ds = de 1 + p 1 d 1 ρ 1 , (33)
where T 1 is the gas temperature. The entropy also satisfies

ρ 1 s ρ 1 = - p 1 ρ 1 s e 1 . (34)
Multiplying the last equation in (31) by s e 1 and the second by s ρ 1 and adding the two we find a new system

U = (α 2 , ρ 1 , u 1 , u 2 , s), U t + C(U )U x = 0, C(U ) =         u 2 0 0 α 2 0 ρ 1 α 1 (u 2 -u 1 ) u 1 ρ 1 ρ 1 α 1 α 2 0 0 p 1,ρ 1 ρ 1 u 1 0 p 1,s ρ 1 Rp α 2 ρ 2 - Rp,α 1 ρ 2 p 1,ρ 1 ρ 2 0 u 2 p 1,s ρ 2 0 0 0 0 u 1         . (35)
The sound speed in the gas is noted c. It is given by

c 2 = ∂ ∂ρ 1 p 1 (ρ 1 , s). (36)
The sound speed a in the powder bed is given by [START_REF] Gallouët | Numerical modelling of two-phase flows using the two-fluid two-pressure approach[END_REF].

In this way, we have

C(U ) =         u 2 0 0 α 2 0 ρ 1 α 1 (u 2 -u 1 ) u 1 ρ 1 ρ 1 α 1 α 2 0 0 c 2 ρ 1 u 1 0 p 1,s ρ 1 a 2 α 2 c 2 ρ 2 0 u 2 p 1,s ρ 2 0 0 0 0 u 1        
The characteristic polynomial is then

(u 1 -λ) (u 2 -λ) 2 -a 2 (u 1 -λ) 2 -c 2 -c 2 α 2 α 1 ρ 1 ρ 2 (u 1 -λ) 2
We define y = λ-u 1 c that has to be a root of

Q 4 (y) = Q 2 (y) Q 4 (y) := y 2 -1 (y -u 2 -u 1 c ) 2 -a 2 c 2 Q 2 (y) := 1-α 1 α 1 ρ 1 ρ 2 y 2
It is easy to check that Q 2 and Q 4 have at least two real intersections. For the two other roots, general analytic conditions cannot be obtained but we can indicate sufficient conditions for hyperbolicity (or ellipticity)

1. a = 0, u 1 = u 2 : the system is hyperbolic.

2. a = 0, u 1 = u 2 , -1 ≤ u 2 -u 1 c ≤ 1: the system is not hyperbolic. Let us observe that in our computations we are generally in the cases 2 or 3. During the ballistic process, the combustion of the powder tends to increase the porosity. Thus, case 3 is observed only at the beginning of the computation, when α 1 ≃ α c (see [START_REF] Gallouët | Numerical modelling of two-phase flows using the two-fluid two-pressure approach[END_REF]). The intra granular stress plays here an important role to stabilize the computation.

Although the rest of the computation occurs in a non-hyperbolic regime we did not observe instabilities. This is certainly due to several dissipative mechanisms:

-the numerical viscosity of the Rusanov scheme is known to be high, compared to other schemes. Typically, if the drag force is set to 0, instabilities begin to appear for 10 000 cells in 1D computations (see [START_REF] Herard | A relaxation method to compute two-fluid models[END_REF]);

-the drag force plays an important role to damp the oscillations as demonstrated by Hérard in [START_REF] Herard | A relaxation method to compute two-fluid models[END_REF]. The scheme can remain stable up to 500 000 cells computations;

-our case corresponds to small velocity differences u 2 -u 1 and the size of the imaginary part of the complex eigenvalues tends to 0 with u 2 -u 1 .

nuss_hell.tex; 10/02/2006; 11:02; p.18 With finer meshes, or higher order schemes, it will be necessary to improve the stability of the initial model. A simple way to do that is to modify the intra granular stress expression [START_REF] Hirschefelder | Molecular theory of gases and liquids[END_REF] in such a way that a is always > 0, as in [START_REF] Combe | Finite volume algorithm to compute dense compressible gas-solid flows[END_REF] or [5].

Figure 1 .

 1 Figure 1. Initial geometry in a gun

Figure 2 .Figure 3 .

 23 Figure 2. Porosity and density

3 .

 3 a < c, |u 2 -u 1 | ≤ max(a, c -a): the system is hyperbolic. This sufficient condition is simply obtained by requiring that Q 4 (0) > Q 2 (0) when all the roots of Q 4 are in [-1, 1].

Table I .

 I Initial conditions, constant velocity, constant pressure. α1 ρ1 (kg.m -3 ) ρ2 (kg.m -3 ) u1 (m.s -1 ) u2 (m.s -1 ) p1 (Pa) p2 (Pa)

	Left	0.5	0.870	1587	100	100	10 5	10 5
	Right	0.8	0.512	1587	100	100	10 5	10 5

Table II .

 II Result of simulations

	Computed value	Acceptable range Algorithm result
	Maximal shot base pressure (M P a)	325 -360	344
	Maximal breech pressure (M P a)	355 -400	377
	Muzzle velocity (m.s -1 )	660 -705	694
	Shot exit time (ms)	14.66 -16.58	15.75

Table III .

 III Muzzle velocity datas

	Initial powder mass	Experimental result Simulation result
	1.17 kg	1000 m.s -1	1001 m.s -1
	1.365 kg	1119 m.s -1	1116 m.s -1
	1.4625 kg	1194 m.s -1	1172 m.s -1
	Table IV. Breech pressure datas
	Initial powder mass	Experimental result Simulation result
	1.17 kg	367 M P a	363 M P a
	1.365 kg	534 M P a	521 M P a
	1.4625 kg	665 M P a	626 M P a

nuss_hell.tex; 10/02/2006; 11:02; p.13

nuss_hell.tex; 10/02/2006; 11:02; p.3

nuss_hell.tex; 10/02/2006; 11:02; p.4Gas-particle flows with combustion

nuss_hell.tex; 10/02/2006; 11:02; p.5

nuss_hell.tex; 10/02/2006; 11:02; p.6

nuss_hell.tex; 10/02/2006; 11:02; p.7

nuss_hell.tex; 10/02/2006; 11:02; p.8

nuss_hell.tex; 10/02/2006; 11:02; p.9

nuss_hell.tex; 10/02/2006; 11:02; p.10

nuss_hell.tex; 10/02/2006; 11:02; p.11

nuss_hell.tex; 10/02/2006; 11:02; p.15

nuss_hell.tex; 10/02/2006; 11:02; p.16

nuss_hell.tex; 10/02/2006; 11:02; p.17

nuss_hell.tex; 10/02/2006; 11:02; p.19

Appendix

A. Data set for AGARD 132 mm

We used a 7-holes powder. Igniter's holes are placed at 5 points x = 0 mm, x = 31.75 mm, x = 63.5 mm, x = 95.24 mm, x = 127 mm.

The computations data are given in Table V.


M Ladieu

Sanquer

Michael Specht

Laurent P Levy

Francois Ladieu

Mare Sanquer

Magnetopolarizability at the Metal-Insulator Transition

Keywords: numbers: 71, 30, +h

In 3D Anderson insulators, each electronic state is localized over a characteristic distance s = I/tc, the localization length, as a result of interferences induced by disorder [1]. Numerical [2] and theoretical [START_REF] Pichard | Quantum Coherence in Mesoscopic Systems[END_REF] studies of the Anderson transition have demonstrated that the values of the inverse of the localization length ~have a Gaussian distribution centered around (~) = I/(g), deep in the insulating regime. As the metal-insulator (MI) transition is approached, the width of the distribution var[t~j becomes comparable to its mean (t~). In this Letter, we demonstrate the experimental relevance of the distribution of localization lengths at the MI transition.

At zero temperature, ac conduction in an Anderson insulator can take place through the polarization current associated with the displacement of localized states. As the transition is approched, this "intrinsic polarizability" e;", diverges [4,5]. At finite temperature, there is an additional contribution to the ac response, because charges can hop from one state to another by phonon assisted tunneling (Mott hopping). An occupied state, localized on a length scale s = I/tc, can make a few "optimal" hops to neighboring empty sites. These hops can be represented by bonds with a conductance work: deviations are seen only at the lowest temperature on the sample closest to the transition. At finite frequency, this network must be modified, since the electric field mod- ulates the occupation n; (and hence the charge en;) of states (i) close to the Fermi level. This introduces a capac- itance [8] C; = e tin;/ciV; = e /4ekttT cosh (E,/2kttT) between each node and ground as shown in the inset of Fig. 1 (e is the dielectric constant averaged over the locali- zation length and F, is the energy of the state i measured from the Fermi level). The typical value of C for states with inverse-localization length ~, 10 -10 0 1 2 0.5 2.5

1.5 = gp exp( -2~8ivt, «), (1) dependent on ~through kttTp(~) = 18~/n(eF) These.

bonds form a Miller-Abrahams resistor network [6]. When the bond conductances are distributed over a sufficiently broad scale, the overall conductance is dominated by the few bonds with inverse localization length K"which estab- lish the sample-spanning cluster [7] over the network. The overall conductance is of the order of g(tc, ), yielding the well-known Mott hopping law [Eq. (1)]. This amorphous material is known to be a good An- derson insulator with strong spin-orbit scattering due to the yttrium [9,10]. In particular, the large electronic density at the Fermi level [n(eF) = 7 X 104s J 'm at the critical concentration x, = 0.18] partially screens electron-electron interactions on distances short compared to s. Finally, the typical values of the localization length, 3 ( ( ( 40 nm, are much larger than all microscopic length scales and can easily exceed the magnetic length accessible in our laboratory. The frequency dependent capacitance and conductance are measured using a three terminal capacitance bridge as in earlier studies [11 -13] in a dilution refrigerator. The three curves taken at 0.59, 0.4, and 0.33 K differ only by the characteristic frequency cu, at which the capacitance drops, and can be overlapped by a simple translation on the logarithmic frequency axis, implying a universal fre quency dependence of the ac response. This is illustrated in Fig. 2(b), where nine data sets for the capacitance and conductance taken at different temperatures and magnetic fields have been overlayed after rescaling the frequency axis by the characteristic frequencies.

This frequency scaling involves a single parameter (iu, ) instead of two as observed in glasses [START_REF] Dixon | [END_REF].

Such behavior can be accounted for [8] by applying an effective medium approximation to the impedance network shown in Fig. 2(b). The total admittance of all the bonds leaving site 0 is replaced by a local effective admittance Y(cu). If the network is locally homogeneous, its admittance seen from the neighboring sites i is also equal to Y(ru). Y(cu) can therefore be determined self- consistently.

The equivalent admittance of all the bonds connected to site 0 leaves its admittance Y(m) invariant:

Y(Kp, M) g(Kp) Y(Kp, Cd) + Cr(Kp) dC (3)
where z is the effective number of bonds and the brackets represent an average over local disorder. This approach can be justified in a percolation picture, where all the bonds of conductance g ( g(~, ) are removed and the bonds of conductance g ) g(~, ) are replaced by a conductance g(i~, ), effectively making the network homogeneous. In this case, there is only one characteristic A characteristic frequency can be identified as the inHection point. At 30 mK, the capacitance is frequency independent: Mott's polarizability is quenched leaving only the "intrinsic" one. (b) Scaling plot of the capacitance and conductance for the same sample as a function of scaled frequency for nine different temperatures and magnetic fields.

(

) 4
where g and C are defined in Eqs. (1) and (2). At finite frequency, the current flows from site 0 through the capacitor Co and in parallel through the network represented by the admittance Y(ir, cu). At frequencies small compared to ~" the network current can only fiow through the critical conductances [Y(x."cu) = g, ].

When cu exceeds ~" the capacitors allow the current to be shunted to ground, effectively breaking the sample in subcritical clusters which are divided further as the frequency increases. The overall frequency dependence of the local admittance is conveniently represented by the scaled admittance y(ru) = [Y(tp) + iCcu] jY(0), which depends only on the scaled frequency 0 = ru/tu, . The self-consistent solution of (3) has been shown [8] to be equivalent to ylny =iA as long as Intel «(To/T)'t . For insulating samples far from the transition, the percolating bonds dominate the impedance network and the sample conductance is proportional to y(rr"or). At low frequency (static limit), y = exp(iA), the capacitance and the conductance of the network are constant. At high frequencies (II » 1), y = iA/In(iA), yielding a sample conductance Re(g) Re(y) = (~/2)A/[InA]2 increasing sublinearly with fre- quency and a capacitance C ~Im(y)/0 = I/InfI de- creasing logarithmically with ~, in excellent agreement with the data shown in Fig. 2(b). At the lowest tem- perature [30 mK, Fig. 2(a)], the capacitance reaches a frequency independent limit, indicating that the intrinsic polarizability takes over Mott s contribution.

In these samples, a magnetic field can induce a large change in eM, «, which saturates at sufficiently high field (cf. Fig. 3). Since data taken in a field can be overlaid on the same universal plot [Fig. 2(b)], this dependence is purely dynamical and can be described by a magnetic field dependence of cu, . This dependence may be attributed to a magnetic field dependence of rr, (see below). Since ~, is here proportional to the static conductance, mag- netocapacitance and magnetoconductance are related to one another. In particnlar, the negative magnetoconductance of these samples [9] is entirely consistent with the decrease of ~, with magnetic field. At high tempera- ture or frequencies small compared to co" the response is static, and changes in characteristic frequencies have little effect on the capacitance [cf. Fig. 2(b)]. When or is comparable to ~"a change in cu, of approximatively one decade is sufficient to go from the low to the high fre- quency regime. This produces a purely dynamical capacitance drop of an order of magnitude. At low temperature or high frequency, changes in cu, have only logarithmic effects [cf. Fig. 2(b)] on the capacitance.

An important field scale, however, can be identified from Fig. 3. For all temperatures, one can define the characteristic field B, for which the capacitance change C(0) -C(B) is half of its total variation C(0) -C(~). The resulting plot of B, as a function of T, shown in the inset of Fig. 3 is a power law with an exponent close to 1/2. This means that the product B,SM."=B, --= 1.3 X 10 "Tm' = 0.33eo This relates directly to the observed magnetocapacitance to quantum-interference effects in an area limited by ZM, « [START_REF] Shklovskii | Hopping Transport in Solids[END_REF]. In insulators, this length plays the role of the phase-coherence length since the phonon absorbed in a Mott hop is an inelastic process [START_REF] Milliken | [END_REF].

Samples close to the MI transition have a different be- havior. Figure 4 shows the frequency dependence of the geneous, with an inverse localization length sc, the overall conductance is the average of I'(tc, co) over a relevant dis- tribution of inverse localization length 2 (tc). Hence, the frequency dependence of the normalized complex conduc- tance g(co)/g(0) is approximatively g(co)/g(0) = 2 (tr)y(tc, co) dtc. the temperature-dependent frequency stretching of the crossover region. After differentiating Eq. ( 4), the vari- ance in characteristic frequencies is var[in co, ] = i6[To((tc))/T]' var[tc]/(t~), (8) which spreads on a logarithmic scale as (To/T)'l . This is illustrated in Fig. 4, where two line segments stretched by the temperature ratios (0.806/0. 103)'t4 = 1.7 of two different conductance curves reproduce the stretching in the crossover region. Equations ( 5) and ( 7) define the solid line fits shown in Fig. 4 using the following fitting parameters: (a) the average localization length (g) = I/(tc) = 30 nm and the conductance prefactor go are determined by the temperature dependence of the conductance (cf. Fig. 1), (b) the dielectric constant e = 500 sets the overall frequency scale, and (c) the width var[tc] = 0.3(tc). No other parameters are necessary to reproduce the magnitude of the capacitance. The overall behavior is well reproduced, except for a slower drop of the capacitance in the high frequency regime. This is not surprising, considering (a) the limited validity of Eq. ( 5) close to the MI transition and (b) the role of the intrinsic polarizability which was left out of the analysis.

The overall frequency dependence in large magnetic fields (not shown) has the same shape as the dependence in zero field but at a lower temperature. For example, in a field of 7 T, the effective temperature drops roughly by a factor of 2. This is not only apparent in the smaller dc value of the conductance (negative magnetoconductance), but also in the additional stretching of frequency scales observed in the frequency dependent capacitance. Several interpretations may be considered. This "effective cooling" in the ac response can be mimicked either by a reduction of the average localization length [increasing (tc) = I/($), see Eq. ( 4)] or by an increase in the width of the distribution var[tc] [cf. Eq. ( 8)]. The observed "cool- ing" in a field of 7 T can be reproduced by a 20% to 30% change of (tc), or by a 15% to 20% increase of var[K].

Such significant changes (considering that they come in the exponent for most experimentally accessible quantities) are consistent with interference effects in Anderson insulators with strong spin-orbit scattering [10]. If the MI is driven only by disorder, a one-parameter scaling theory [1]is applicable and changes in var [tc] and (tc) are related and occur sirnultanously.

In conclusion, the ac response of Anderson insulators is sensitive not only to the zero-temperature displacement current, but also to an activated or "Mott" polarizability. The frequency dependence of the real and imaginary part of the conductivity probes the relevant length scales in the problem. Deep in the insulating regime, the localization length is the only relevant length scale observed. A temperature-dependent characteristic magnetic field can be identified, consistent with interference effects within the Mott hopping length. We relate the observed stretch- ing of the frequency scales close to the MI transition to the distribution of inverse localization lengths at the tran- sition. A magnetic field modifies this distribution, but it is at present not possible to distinguish between a shift or a broadening of the distribution.

We wish to thank the DRET for financial support under Grant No. 92/181. One of us (M. Sp.) acknowl- edges financial support from the DAAD foundation (HSPII/AUFE program).

  g(tc) = gp exp( -[Tp(K)/T]' InC(t~) = (InC;(t~))

Figure 2 (

 2 Figure 2(a) shows the typical frequency dependence of the capacitance at different temperatures.

10 FIG. 2 .

 102 FIG. 2. (a) Frequency dependence of the capacitance of samples deep in the insulating regime. A characteristic

  , which, from the experimentally determined value of To and s (cf. Fig.1), is of the order of the flux quantum.

FIG. 3 .FIG. 4 .

 34 FIG. 3. Magnetic field dependence of the capacitance (mea- sured at fixed frequency) of samples deep in the insulating regime at four different temperatures.Inset: Temperature de- pendence of the characteristic field 8,.

  and its width var[tc], 2 (v) ~tc'exp(z(tc -(tc)) /var[t~]) is used to fit the experimental data. The width vat[tc] reproduces


P Pari

P Forgep

Ph Hernandez

F Ladieu

M Sanquer

I M Broto B'

H Rakoto

G Coffe

S Askenazy

Design of a dilution refrigerator and sample holder operating in pulsed field

We have constructed a compact dilution refrigerator which can lower the temperature from 300 K to 50 mK in less than 3 h, and operate in a pulsed magnetic field. To avoid eddy currents due to contacts with the sample we design an original sample holder.

We describe how to measure magnetoresistance in pulsed magnetic field up to 36 T at T = 60 mK without significant heating of the sample.

To avoid large eddy currents, the metallic parts of the refrigerator (4 K flange, Joule-Thomson exchanger, metallic parts of the composite boiler and heat exchangers) are mainly situated outside the magnet. But, even with the use of a complete plastic lower part (mixing chamber and its connections), metallic wires to the sample are possible sources of heating in large varying flux. It is important to avoid the diffusion of this heat to the sample during the measurement time of approximately 1 s. For this purpose, we design a composite sample holder which includes a quartz radiator (dissipating the heat coming from the upper wires in short times) in series with highly resistive electrical contacts down to the sample (in order to reduce the diffusion of residual heat). To test the whole system, we measure the magnetoresistance of a sample exhibiting a large temperature dependence of its resistance. We compare the results in pulsed field with that in a continuous 8 T field and obtain that * Corresponding author. the sample temperature is increased from 60 to 70 mK during the pulse.

The characteristics of the pulsed field are detailed in Ref. [1]: the field is obtained by discharge of a capacitor bank storing 1.25 MJ at 10 kV. The coil produces 42 T in an inner diameter of 15 mm at T = 4.2 K. The field first increases to its maximum value in 0.1 s, and then decreases in 1 s. During this decrease, we measure the voltage across the sample either with DC or AC fixed injected current and 4-probe technique.

We use a dilution refrigerator without separated condensation stage, because of its compactness and its easyto-use qualities. The cooling time from room temperature to T = 50 mK is only 3 h. The cooling power at T = 100 mK is 20 laW.

Due to the small available diameter, we cannot build a plastic heat exchanger near the mixing chamber, but we construct composite fiberglass-epoxy resin tubes (diameter 1 mm) to connect the mixing chamber to the metallic heat exchanger located out of the coil where the residual field is only 1/100 the nominal value (see Fig. 1).

The mixing chamber itself is from epoxy resin and sample and thermometer are placed inside the helium mixture (Fig. 2 Due to the low heat capacity of the quarl~z, the heat coming to the pads diffuses very rapidly in! the quartz (large diffusivity), and then is dissipated in a short time (small time constant of thermal coupling tq the bath). Moreover, because the quartz itself has a large contact area with the helium bath, a large volume Of the liquid is involved in dissipation during a small time; this is necessary to maintain the quartz at low temperature. We have calculated an increase of 5 mK at 60 mK for 1 ~J of dissipated energy by eddy current in the metallic part of the sample holder into the mixture inside the mixing chamber. In addition, there is energy dissipated in the metallic parts of diffusion stage outside the coil, and the associated heating of the helium bath diffuses only after 10 s into the mixing chamber due to the flow rate of the concentrated phase.

To test the dilution refrigerator in the coil, we have measured both a RuO2 calibrated sensor (1.58 kf~ Dale resistor) and our sample. The sample itself is a very good thermometer due to the large temperature dependence of its resistance. The fit at H = 0 is given by R = 27.71 exp(3.145/T) 1/4. It consists of an amorphous Yo.2Sio.a alloy which is in the variable Range Hopping Mott regime below T = 1 K [2]. In addition it presents a very large positive magnetoresistance, such that one discriminate between heating, which produces a negative variation of the resistance, and magnetoeffect. First we measure the magnetoresistance in constant magnetic field of 0 and 7 T at 60 inK: with the field the resistance increases from 400 to 947 l'~. This magnetoresistance is reported on Fig. 3 (O).

Then we measure the whole magnetoresistance curve in a pulsed magnetic field at the same temperature (i.e. we choose the same H = 0 sample resistance). We observe coincidence to well within 10 mK of the measured values of AR (H) at H = 7 T. In addition, we make another pulsed field experiment with the variation of flux increased by a factor 3, and confirm that the measured magnetoresistance depends only smoothly on this factor. To conclude, the value of R(7 T ) obtained with the faster field sweep is the same as obtained in DC field at 7 T and 70 mK (830 ~ measured).

Fig. 1 .Fig. 2 .

 12 Fig. 1. Lower part of dilution refrigerator in the coil.

Fig. 3 .

 3 Fig. 3. AR/R = [R(H) -R(O)]/R (0) versus H l]z a! T = 60 mK in constant field of 7 T (O), and for 2 pulsed field experiments.


Thierry Huillet

Servet Mart Ínez

Martin M Öhle

ON POLYMORPHISM FOR DISCRETE EVOLUTIONARY DYNAMICS DRIVEN EITHER BY SELECTION OR SEGREGATION DISTORTION

Keywords: species frequencies dynamics, selection, segregation, polymorphism, stability, potential and stochastic matrices

We revisit some problems arising in the context of multiallelic discrete-time and deterministic evolutionary dynamics driven first by fitness differences and then by segregation distortion. In the model with fitness, we describe classes of fitness matrices exhibiting polymorphism. In the segregation case, still in search for conditions of polymorphism, we focus on a class of skew-symmetric matrices with a unique strictly positive kernel vector. Our main results reduce the study of these cases to the analysis of stochastic matrices.

Introduction

We will first briefly revisit the basics of the deterministic dynamics arising in discrete-time asexual multiallelic evolutionary genetics driven only by fitness differences (or selection) and we will mainly consider the diploid case with K alleles. In the diploid case, there is a deterministic updating dynamics of the full array of the genotype frequencies that involves the nonnegative symmetric fitness matrix W 0 (with W = W , W the transpose of W ) attached to the genotypes. When the Hardy-Weinberg law applies, it is sufficient to treat the induced marginal allelic frequencies dynamics. The updating dynamics on the simplex involves the mean fitness as a quadratic form in the current frequencies whereas marginal fitnesses are affine functions in these frequencies. The induced dynamics turns out to be gradient-like. We will also consider an alternative updating mechanism of allelic frequencies on the simplex, namely segregating distortion: here the fitness matrix W is based on skew-symmetric matrices and the fitness landscape is flat. The induced relative frequencies dynamics is divergence-free like.

In the first case driven by viability selection and in the diploid context, the Fisher theorem, stating that the mean fitness increases as time passes by, holds true but, as a result of the fitness landscape being quadratic, a polymorphic equilibrium state will possibly emerge. Due to its major evolutionary interest, our subsequent concern is to identify examples of diploid dynamics leading to a unique polymorphic state x eq on the simplex, either unstable or stable. We start with the unstable case and draw the attention on a class of fitness matrices leading to a unique unstable polymorphic equilibrium state: the class of potential matrices which includes the class of ultrametric matrices. For such potential fitness matrices, the mean fitness quadratic form is positive-definite. We will then consider a related class of fitness matrices with a mean fitness quadratic form which is negative-definite on a hyperplane of codimension one. For this class of matrices, there will also be a unique polymorphic equilibrium state for the diploid dynamics and it will be stable.

In the second case driven by segregation distortion, fitness is flat in the sense that its mean fitness quadratic form is a constant (the segregation ratio or fraction of A k gametes contributed by A k A l individuals is W k,l ≥ 0 with W k,l + W l,k =Const). Conditions under which its skew-symmetric part (W -W ) /2 has a strictly positive kernel will be addressed. If this is indeed the case, such models will also display a polymorphic equilibrium state and the flat fitness dynamics will spiral away from it to eventually reach the boundary of the state-space (as a simplex) in infinite time. In [START_REF] Weissing | Selection and segregation distortion in a sex-differentiated population[END_REF], a population without sex-differentiation driven simultaneously by viability selection and segregation distortion is considered in Section 2, together with the adjunction of the effect of sex in Section 3. In this note, we shall study separately systems driven first by selection (in the absence of segregation distortion when segregation is Mendelian W k,l = W l,k =Const/2) and then by segregation distortion (in the absence of selection) and we shall not consider the effect of sex.

Our main results will therefore characterize different types of nonnegative or skewsymmetric matrices (relevant for polymorphism) in terms of stochastic matrices. Theorem 1 of Section 2 characterizes the symmetric nonnegative fitness matrices for which x eq is an isolated polymorphic stable state, in terms of stochastic matrices. In Theorem 2 of Section 4, we characterize skew-symmetric matrices having a strictly positive probability kernel x eq also in terms of stochastic matrices. The associated dynamics with a flat fitness landscape will then display an isolated polymorphic state x eq , but it will be unstable.

2. Single locus: diploid population with K alleles 2.1. Joint and marginal allelic dynamics [START_REF] Bürger | The mathematical theory of selection, recombination, and mutation[END_REF][START_REF] Ewens | Mathematical population genetics. I.Theoretical introduction[END_REF]). Consider K alleles A k , k ∈ I K := {1, ..., K} attached at a single locus. Let W = (W k,l : k, l ∈ I K ) be some nonnegative fitness matrix. The coefficient W k,l ≥ 0 stands for the absolute fitness of the genotypes A k A l . Since W k,l is proportional to the probability of an A k A l surviving to maturity, it is natural to assume that W is symmetric.

Let X = (x k,l : k, l ∈ I K ) be the current frequency distribution at (integer) time t of the genotypes A k A l , so x k,l ≥ 0 and k,l x k,l = 1. The joint evolutionary fitness dynamics in the diploid case is given by X(t + 1) = P(X(t)), with (1)

P(X) k,l = x k,l W k,l ω(X) and ω(X) = k,l x k,l W k,l .

The relative fitness of

A k A l is W k,l /ω(X)
and ω(X) is the mean fitness. The genotypic variance in absolute fitness is

σ 2 (X) = K k,l=1 x k,l (W k,l -ω(X)) 2 and the diploid variance in relative fitness is σ 2 (X) = σ 2 (X)/ω(X) 2 . Note that (2) ∆ω(X) = k,l ∆x k,l W k,l = k,l x k,l W 2 k,l ω(X) -W k,l = ω(X)σ 2 (X) ≥ 0
with a relative rate of increase: ∆ω(X)/ω(X) = σ 2 (X). This expression vanishes only when W k,l is constant. This is the full diploid version of the Fisher theorem for asexuals; see Okasha 2008 for a deeper insight on the meaning of this theorem.

Assuming Hardy-Weinberg proportions, the frequency distribution at time t of the genotypes A k A l is given by: x k,l = x k x l where x k = l x k,l is the marginal frequency of allele A k in the whole genotypic population. The frequency information is now contained in x = X1 1 , where 1 := (1, ..., 1) is the 1-row vector of dimension

K. x := (x k : k ∈ I K ) belongs to the K-simplex S K = {x := (x k : k ∈ I K ) ∈ R K : x 0, |x| = 1}.
Here |x| := K k=1 x k and x 0 means that all x k ≥ 0, k = 1, ..., K. The elements of S K are called states. The mean fitness is given by the quadratic form:

ω(x) := k,l x k x l W k,l = x W x. We now have σ 2 (x) = K k,l=1 x k x l (W k,l -ω(x))
2 and

σ 2 (x) = σ 2 (x)/ω(x) 2 .
We will consider the update of the allelic marginal frequencies x themselves. Define the frequency-dependent marginal fitness of

A k by w k (x) = (W x) k := l W k,l x l .
For the vector x denote by D x :=diag(x 1 , ..., x K ), the associated diagonal matrix. The marginal mapping p : S K → S K of the dynamics is given by:

(3)

x(t + 1) = p(x(t)), where p(x) = 1 ω(x) D x W x = 1 ω(x) D W x x .
This dynamics involves a multiplicative interaction between x k and (W x) k , the kth entry of the image W x of x by W and a normalization by the quadratic form ω(x) = x W x. Iterating, the time-t frequency distribution is:

x(t) = p t (x(0))
where x(0) ∈ S K is some initial condition and p t = p • ... • p, t-times.

Remark. If W k,l = w k w l , then (3) boils down to

x(t + 1) = p(x(t)), where p(x) = 1 w(x) D w x = 1 w(x) D x w.
Here w = (w 1 , ..., w K ) and w(x) = l w l x l is linear. This yields the haploid asexual selection dynamics.

For an alternative representation of the dynamics (3) take ∆x = p(x) -x and define the symmetric positive-definite matrix G(x) = D x (I -1x ) with entries:

G(x) k,l = x k (δ k,l -x l ) . Let V W (x) = 1 2 log ω(x).
Then, (3) may be recast as the gradient-like dynamics [START_REF] Shashahani | A new mathematical framework for the study of linkage and selection[END_REF]):

(4) ∆x = 1 ω(x) G(x)W x = G(x)∇V W (x), with |∆x| = 1 ∆x = 0 as a result of 1 G(x) = 0 . Note that ∇V W (x) ∆x = ∇V W (x) G(x)∇V W (x) ≥ 0.
1 Throughout, a boldface variable, say x, will represent a column-vector and its transpose, say

x , will be a row-vector. And x 0 means all entries of x are positive.

One can easily check that (3), as a replicator dynamics, can also be recast under the (non-linear) more conventional Fokker-Planck-like form x (t + 1) = x (t)P (x(t)), where P (x) is the x-dependent matrix with (k, l)-entry P (x) k,l = W k,l x l /x W x ∈ (0, 1). This matrix is not stricto sensu stochastic, but we note however that x P (x)1 = 1 when x ∈ S K .

The mean fitness ω(x), as a Lyapunov function, ∆ω(x) ≥ 0, see [START_REF] Atkinson | A matrix inequality[END_REF]Kingman 1961b). Its partial rate of increase due to frequency shifts only is δω(x)

:= k ∆x k w k (x). It satisfies (5) δω(x) ω(x) = k x k w k (x) ω(x) -1 2 = k (∆x k ) 2 x k = σ 2 A (x) 2 , σ 2
A (x) being the allelic variance in relative fitness. See [START_REF] Price | Fisher's fundamental theorem made clear[END_REF][START_REF] Ewens | Mathematical population genetics. I.Theoretical introduction[END_REF][START_REF] Castilloux | The fundamental theorem of natural selection in Ewens' sense (case of many loci)[END_REF] and the discussion of these results in (Bürger 2000).

Equilibria (diploid case).

The mean fitness increase phenomenon occurs till the evolutionary dynamics reaches an equilibrium state. We wish to briefly discuss the questions relative to equilibria in the diploid case.

In contrast with the haploid case, in the diploid situation, the dynamics (3) can have more complex equilibrium points, satisfying w k (x eq ) = w 1 (x eq ), k = 2, ..., K and |x eq | = 1. To avoid linear manifolds of equilibria, we first assume that all principal minors of W are nonsingular. In particular W k,k > 0 for all k. In this case, from the Bézout theorem, the number of isolated equilibria is bounded above by the number 2 K -1 of faces of S K , see [START_REF] Tallis | Equilibria under selection for k alleles[END_REF].

When W = I is the identity matrix, there are 2 K -1 equilibrium points (the barycenters of the K l+1 l-dimensional faces, l = 0, .., K -1), but only one polymorphic equilibrium which is the barycenter x B of S K . The barycenter has minimal mean fitness and it is unstable. The 0-faces are stable fixed points whereas the barycenters of the l faces with l ∈ I K-2 are unstable saddle-points. The simplex S K can be partitioned into the attraction basins of the stable 0-face states.

Similar conclusions can be drawn when we start with W = (I -D λ ) -1 , where λ := (λ k : k ∈ I K ) satisfies 0 λ ≺ 1 (meaning 0 ≤ λ k < 1, ∀k). In this case again, there is only one unstable polymorphic equilibrium which is x eq = (1-λ)/(K -|λ|) ∈ S K . Again, all fixed points but the 0-faces are unstable and the simplex S K can be partitioned into the attraction basins of these stable 0-face monomorphic states: however, in contrast with the haploid case, the type of the survivor is not necessarily the one of the fittest (with largest λ k ); it will depend on the initial condition.

Due to its evolutionary interest, we would like now to discuss the stability of a polymorphic state. Under the above assumptions on the principal minors of W , a unique polymorphic equilibrium state exists if and only if (i) there is a unique z 0 for which W z = 1.

Under this condition, the unique equilibrium polymorphic state is

x eq = z |z| .
Note that its mean fitness is given by ω(x eq ) = 1/|z|.

Let us analyze the stability of x eq . Under (i), by Cramer's rule, with W (k) a K ×K matrix obtained from W while substituting its k-th column by 1 (satisfying

W (k) | = (adj(W )1 k
, where adj(W ) is the adjugate matrix of W , (see Nagylaki 1992, page 56), we have |W | := det(W ) = 0 and ( 6)

z k = |W (k) | |W | > 0, k ∈ I K .
In (Kingman 1961a and1961b), the condition to have stability of the polymorphic state x eq = z/ |z| was given. This occurs only when the following additional condition holds, [START_REF] Lessard | A criterion for stability-instability at fixation states involving an eigenvalue one with applications in population genetics[END_REF] (ii) W has exactly one strictly positive dominant eigenvalue and at least one strictly negative eigenvalue or else the sequence of the leading principal minors of W alternates in sign, starting positive [START_REF] Mandel | The stability of a multiple allelic system[END_REF].

Under (i) and (ii), x eq is a stable equilibrium polymorphic state, and the stability is in the sense that it is a local maximum of the mean fitness ω(x), so a perturbation of it will result in a mean fitness decrease. The condition for x eq to be an isolated global maximum of the fitness landscape ω(x) is that for all x ∈ S K with x = x eq we have x W x < x eq W x eq . A direct computation shows (Kingman 1961a)

(7) ∀y = 0, y ⊥ 1 : y W y < 0.
In this case, starting from any initial condition in the interior of S K , all trajectories will be attracted by x eq . In Theorem 1 below, we characterize the symmetric nonnegative fitness matrices for which x eq is an isolated polymorphic stable state, in terms of stochastic matrices.

When there is no such unique globally stable polymorphic equilibrium, all trajectories will still converge but perhaps to a local equilibrium state where some alleles have gone extinct [START_REF] Lyubich | Selection-Induced Convergence to Equilibrium in a Single-Locus Autosomal Population[END_REF]. Which allele and how many alleles are concerned seems to be an unsolved problem in its full generality.

We note that when the leading principal minors of W are all positive, x eq is unstable and the quadratic form x W x is positive-definite.

In the sequel we will denote by J := 11 the flat K × K matrix whose entries are all 1.

Theorem 1. Let W be a nonnegative symmetric matrix and suppose there exists z 0 such that W z = 1. Then, W admits the representation W = P D -1 z , where P is a stochastic irreducible matrix obeying detailed balance with respect to z , and conversely.

Assume W satisfies the above conditions. Then, (i) W is positive-definite if and only if P has a positive spectrum;

(ii) The matrix W satisfies sup{y W y : y ⊥ 1, y = 0} < 0, if and only if P has 1 as unique positive eigenvalue, all the other ones being negative. Therefore, this is a sufficient condition in order that x eq be a polymorphic stable state for W ;

(iii) If P is a stochastic matrix that satisfies detailed balance with respect to some z 0 and λ > λ 0 , where λ 0 = min(spectrum(P )) then the stochastic matrix P λ = (1 + λ 0 ) -1 (P + λ 0 I) satisfies detailed balance with respect to z and it has a positive spectrum. Furthermore, every stochastic matrix that satisfies detailed balance and has positive spectrum can be constructed in this way;

(iv) Assume W is a positive-definite nonnegative symmetric matrix for which there exists z 0 with W z = 1. For ξ > max i,j W i,j , define

(8) W = ξJ -W.
Then W is a nonnegative symmetric matrix for which z/|z| is a polymorphic stable state.

Proof: Assume that W ≥ 0 is a symmetric matrix and W z = 1 for some z 0.

Then W D z 1 = 1 and so P = W D z is a stochastic matrix; clearly, W = W is equivalent to P = D z P D -1 z . Since P z = z, z is proportional to the invariant probability measure of P . The reciprocal follows immediately, because detailed balance for the stochastic matrix P means P = D z P D -1 z and so W = P D -1 z is a nonnegative symmetric matrix for which W z = 1.

Let us now prove condition (i), namely that W has a positive spectrum if and only if P also does. Since P obeys detailed balance we have that

Q = D 1/2 z P D -1/2 z
is symmetric with the same spectrum as P . We have

W = P D -1 z = D -1/2 z QD -1/2 z
. Hence, W is positive-definite if and only if Q is positive-definite, and this happens if and only if Q has a positive spectrum by an argument given in (Kingman 1961b). And so if and only if P has a positive spectrum.

Let us show (ii) which is a rephrasing of (Kingman 1961a, section 3). From P 1 = 1, we get that the matrix Q has eigenvalue 1 with eigenvector D 1/2 z 1. Then, all other eigenvalues of P are negative if and only if all other eigenvalues of Q are negative and since Q is symmetric, this condition is equivalent to sup{y Qy : y ⊥ D

1/2 z 1, y = 0} < 0. Substituting y for D 1/2
z y, and because z > 0, this is equivalently

{y D -1/2 z QD -1/2 z y : y ⊥ 1, y = 0} < 0. Recalling Q = D 1/2 z P D -1/2 z , this is also y P D -1 z y : y ⊥ 1, y = 0 < 0. Which is sup{y W y : y ⊥ 1, y = 0} < 0, because W = P D -1
z . This is the condition [START_REF] Ganikhodzhaev | Quadratic stochastic operators and processes: results and open problems[END_REF] in order that x eq be a polymorphic stable state for W .

Part (iii) is a direct computation.

Let us prove part (iv). Let W be a positive-definite nonnegative symmetric matrix and z 0 such that W z = 1.

Then W = ξJ -W with ξ > max ij W i,j satisfies W 0 and W z = (ξ|z| -1)1. Note that ξ > max i,j W i,j and W z = 1 imply ξ|z| > 1.
Moreover, since W is positive-definite, we get that for all vector y ⊥ 1 we have y W y = y (ξ11 -W ) y = -y W y < 0, so the stability condition ( 7) is satisfied. 2

We note that in (Kingman 1961a, section 3), the symmetric matrix Q was used for the study of stability.

In the case when W is a symmetric bi-stochastic matrix we have W 1 = 1, and since D 1 = I, the identity matrix, we get P = W in the representation given in the theorem. Let S be a symmetric permutation matrix so S = S = S -1 and it is bistochastic. Let y 0 and consider W = D -1 y SD -1 y . Putting z = W -1 1 = D y Sy, we have z 0 and W z = 1. In this case D z = D y SD y and so P = W D z is the identity matrix, which is stochastic and every positive vector satisfies the detailed balance with respect to it, in particular z.

Potential matrices.

Let us introduce potential matrices associated to strictly substochastic matrices. Let Λ be a strictly substochastic matrix, so Λ 0 and 0 ≺ Λ1 := q ≺ 1. We will also assume that Λ is irreducible. Then, W = (I -Λ) -1 exists and it is the potential matrix associated to Λ (else the resolvent matrix R Λ (z) = (zI -Λ) -1 , evaluated at z = 1). We have W 0 and W is obviously nonsingular. For any λ > 0, the matrix W = λ -1 (I -Λ) -1 will also be called a potential matrix, parameterized by λ.

The equilibrium potential (vector) of a potential matrix W is by definition

z = W -1 1. Hence, for all λ > 0, λ1 z = λ(I -Λ)1 = λ(1 -q) 0 is satisfied. We have: |z| = λ(K -|q|) and (9) z = λ(I -Λ)1 = λ(1 -q) 0.
When Λ is symmetric it has real eigenvalues and all them are strictly smaller than 1. Hence, the symmetric potential matrix W is positive-definite because when solving W u = ξu we get Λu = (1 -1/(ξλ))u and we deduce ξ > 0. The potential matrices W are row diagonal dominant in the sense that

W k,k = max l W k,l .
We can make the construction done in Theorem 1 part (iv), that is matrices of type [START_REF] Hofbauer | Evolutionary games and population dynamics[END_REF]. Let W = λ -1 (I -Λ) -1 be a symmetric potential matrix. Then

z = λ(1-q) 0 satisfies W z = 1. For ξ > max k,l W k,l = max k W k,k
, define the matrix W as in [START_REF] Hofbauer | Evolutionary games and population dynamics[END_REF]:

W = ξJ -W . Then, W 0 and W z = (ξ -|z|)1. We note that ξ > max k,l
W k,l and W z = 1 implies that ξ -|z| > 0, and so z = z(ξ -|z|)z is an equilibrium potential of W . From Theorem 1 parts (ii) and (iv), we get that x eq is a stable polymorphic state for W and ω(x eq ) = (ξ -|z|)/|z| is the mean fitness at equilibrium.

Ultrametric matrices.

Let us introduce ultrametric matrices, which in [START_REF] Martínez | Inverse of strictly ultrametric matrices are of Stieltjes type[END_REF]) were proven to be potential matrices. The symmetric nonnegative ma-

trix W = (W i,j : i, j ∈ I K ) is called an ultrametric matrix if it satisfies ∀j, k, l ∈ I K : W k,l ≥ min {W k,j , W j,l } . This condition implies W k,k ≥ max l =k {W k,l }, ∀k ∈ I K . If W satisfies this last condition but in a strict way, that is W k,k > max l =k {W k,l } ∀k ∈ I k , then
W is called a strictly ultrametric matrix. An ultrametric matrix is nonsingular if no row and no column is equal to 0, and no two rows (or columns) are equal. So, a strictly ultrametric matrix is nonsingular. If W is a nonsingular ultrametric matrix, it is a potential matrix having a nonnegative equilibrium potential z = W -1 1 and if W is a strictly ultrametric matrix then z 0.

When K > 1, by a permutation of indices in I K , an ultrametric matrix can always be written in a nested block form (NBF), meaning that we can write it in the form,

(10) W = W 1 α W 1 I 1 1 I 2 α W 1 I 2 1 I 1 W 2 ,
where

I K = I 1 ∪ I 2 , α W = min i,j
W ij and W 1 and W 2 are ultrametric matrices in NBF. Here, the vector 1 I l denotes the 1-constant vector of dimension |I l |, l = 1, 2.

Then,

W = α W 11 + W 1 0 0 W 2 ,
where

W i := W i -α W 11 , i = 1,
W 1 = α W 1 1 I 1 1 I 1 + W 3 0 0 W 4 ,
and so, with 1 I 1 the vector of dimension |I| whose coordinates in I 1 are equal to 1 and those in I \ I 1 are equal to 0, we get

W = α W 11 + α W 1 1 I 1 1 I 1 +    W 3 0 0 0 W 4 0 0 0 W 2    .
If (extending what we have already used for l = 1 to all l) we define 1 I l to be the vector of dimension |I| whose coordinates in I l are equal to 1 and those in I \ I l are equal to 0, iterating the above construction with K > 1, we can partition W while generating a sequence of partition matrices W 1 , ..., W r with r ≥ K -2, till we end in a diagonal matrix D. But only K -2 of these matrices are not 1 × 1 matrices and so can be subsequently partitioned. Let us denote this subsequence by W l1 , ..., W l K-2 . Starting from W , and ending with the diagonal matrix, we obtain

W = α W 11 + K-2 r=1 α W lr 1 I lr 1 I lr + D.
Note that all the diagonal elements of the terminal matrix D are strictly positive if W is strictly ultrametric.

Stated differently, the above construction gives a real sequence s l ≥ 0 and a sequence of vectors u l for l = 1, ..., 2K -1 with entries in {0, 1} such that (see Nabben and Varga 1994)

(11) W = 2K-1 l=1 s l u l u l .
For strictly ultrametric matrix all the 1×1 matrices corresponding to the K diagonal terms should be strictly positive, so there should be at least K terms s l > 0.

If W is an ultrametric matrix, then the associated matrix W defined in [START_REF] Hofbauer | Evolutionary games and population dynamics[END_REF], satisfies the following ultrametric property,

(12) ∀j, k, l : W k,l ≤ max W k,j , W j,l ,
and if W is strictly ultrametric then we have a strict fitness domination of the heterozygotes A k A l over the homozygotes, ( 13)

∀k : W k,k < min l =k W k,l .
Then, W will display a stable polymorphic equilibrium state x eq .

Remark. The strictly ultrametric conditions ( 12) and ( 13) should be compared with the so-called 'triangle inequality' conditions, which are necessary for stable polymorphism, as first pointed out in [START_REF] Lewontin | Heterosis as an explanation for large amounts of genic polymorphism[END_REF]). They read:

(i) ∀k = l ∃j, k = j = l with: W k,l < W k,j + W l,j ; (ii) ∀k = l : W k,l > 1 2 W k,k + W l,l .
It is clear that the matrices W of type ( 8) derived from ultrametric matrices W will satisfy these conditions. 2

Assume K ≥ 2 and W is a strictly ultrametric matrix, then

max{W ij : i, j ∈ I K } < 2K-1 l=1 s l .
We can construct W = ξJ -W in (8) with ξ = 2K-1 l=1

s l and so

W = s 1 J + V with V := 2K-1 l=2 s l (J -u l u l ) .
So, V is of type (8) as well, then there is a z

V 0 such that V z V = 1. Thus W z V = (s 1 J + V )z V = (s 1 |z V | + 1)1 showing that, with z = z V / (s 1 |z V | + 1), W z = (s 1 J + V )z = 1.
In the case s l = 1, ∀l, any positive matrix of the form W = γJ + V with γ > 0 will admit a stable polymorphic equilibrium state and the equilibrium mean fitness for

x eq = z V /|z V | is x eq W x eq = (γ|z V | + 1)/ |z V |. 2.5. Examples. (i) Let K = 2, s > -1, h > 0 and sh > -1. The positive matrix W = 1 + s 1 + sh 1 + sh 1
defines the fitness matrix with selection parameter s and dominance h. Here x eq = (h/(2h -1); (h -1)/(2h -1)). The matrix W is strictly ultrametric if and only if s < 0 and h > 1, and in this case the equilibrium state is unstable. The matrix W is of type [START_REF] Hofbauer | Evolutionary games and population dynamics[END_REF] if and only if s > 0 and h > 1, and in this case the equilibrium state is stable.

(ii) Let K = 4, that is four alleles, with W = 7 l=1 s l u l u l , being s l = 1 for l = 1, .., 7, u 1 = (1, 1, 1, 1), u 2 = (1, 0, 0, 0), u 3 = (0, 1, 1, 1), u 4 = (0, 1, 1, 0), u 5 = (0, 0, 0, 1), u 6 = (0, 1, 0, 0), u 7 = (0, 0, 1, 0) and ξ = 7. We find W = ξJ -W = J +V with:

V =     2 1 1 1 1 4 3 2 1 3 4 2 1 2 2 3     ,
which is a symmetric matrix of type [START_REF] Hofbauer | Evolutionary games and population dynamics[END_REF]. For this example, z V = 1 150 (27, 1, 1, 3) and x eq = 1 32 (27, 1, 1, 3) and the equilibrium mean fitness is x eq W x eq = (|z V | + 1)/|z V | = 91/66. 2

Flat fitness and segregation distortion (Weissing and van Boven 2001)

We now address the flat fitness model. Let A be a real skew-symmetric matrix, so obeying A = -A.

Let σ, λ > 0. We wish here to consider evolutionary dynamics of the form (3) but now when W is of the form W = λ -1 (J + σA) 0 when A = -A and such that |A k,l | ≤ 1/σ. The mean fitness function ω(x) appearing in ( 3) is a constant ω(x) = x W x = λ -1 , and in this sense the fitness matrix W is called flat. Because W k,l + W l,k = 2λ -1 , these models correspond to constant-sum games in which each pair of two players has opposed interest or to evolution under the effect of segregation distortion in population genetics; See [START_REF] Karlin | Mathematical models, problems, and controversies of evolutionary theory[END_REF], Hofbauer and Sigmund 1998, Weissing and van Boven 2001). An interesting sub-family of such models is when σ ∈ (0, 1] and A k,l ∈ {-1, 0, 1}. Some authors recently and interestingly considered the interaction between selection differences and segregation distortion effects but we shall not run into this problem (Weissing and van Boven 2001; Sanyal and Sarkar 2013).

The dynamics (3) for this particular form of W boils down to [START_REF] Kingman | Mathematics of genetic diversity[END_REF] x(t + 1) = p(x(t)), where p

(x) = 1 ω(x) D x W x = x+σD x Ax, independent of λ. With ∆x k (t) = x k (t + 1) -x k (t), with H (x) := K k=1 x 2
k /2 the simple energy quadratic form, we have [START_REF] Lessard | A criterion for stability-instability at fixation states involving an eigenvalue one with applications in population genetics[END_REF] ∆x k (t)

x k (t) = (W x) k -1 ω(x) = σ(A∇H(x)) k , k = 1, ..., K,
with σ∇ • (A∇H (x)) = 0: the relative allelic frequencies obey a divergence-free dynamics.

Such p(•) are called quadratic stochastic or Volterra operators, [START_REF] Ganikhodzhaev | Quadratic stochastic operators and processes: results and open problems[END_REF]. When σ = 1/ (max A k,l + c) with c ≥ 0 some constant, it can be checked that the K extremal points (0-faces) of S K are fixed points of [START_REF] Kingman | Mathematics of genetic diversity[END_REF].

A key problem we shall face in the next Section is to determine those skewsymmetric matrices A for which Az = 0 for some z 0. Indeed, when such a z exists, with λ = |z|, the new W = λ -1 (J + σA) obeys W z = 1, with z an 'equilibrium potential' of W ; x eq := z/|z| will be a polymorphic equilibrium state of the dynamics [START_REF] Kingman | Mathematics of genetic diversity[END_REF] governed by the flat fitness matrix W .

Note that since ω(x) is constant, it is not a Lyapunov function. Define the partial rate of increase of the mean fitness due to frequency shifts as

δω(x) := k ∆x k w k (x), where w k (x) = (W x) k . Observing x Ax = k x k (Ax) k = 0, it is satisfied δω(x) ω(x) = δω(x) = |z 0 | k ∆x k (W x) k = |z 0 | k ∆x k (1 + (Ax) k ) = |z 0 | k x k (Ax) k 1 |z 0 | + (Ax) k = |z 0 | k x k (Ax) 2 k = σ 2 A (x) 2 > 0,
where σ 2 A (x) is the allelic variance in relative fitness [START_REF] Price | Fisher's fundamental theorem made clear[END_REF]). Consider the dynamics p

(x) = 1 ω(x) D x W x = 1 ω(x) D W x x, where W = J + A 0, with A skew-symmetric, |A k,l | ≤ 1 and such that there is a z 0 for which Az = 0. So p(x) = x + D x Ax. Consider V (x) = K k=1 x z k k , where x ∈ S K . Then V is maximum at x eq = z/|z| with V (x eq ) = |z| |z| K k=1 z z k k . Clearly, in view of (W x) k ≥ 0 and x W x = 1, V (p(x)) = V (x) K k=1 (W x) z k k = V (x) K k=1 (1 + (Ax) k ) z k < V (x),
the latter inequality following from the information inequality (Bapat and Raghavan 1997, p. 80), the concavity of the function log and the Jensen's inequality.

This shows [START_REF] Liberman | On the relation between the instability of ESS in discrete dynamics and segregation distortion[END_REF][START_REF] Ganikhodzhaev | Quadratic stochastic operators and processes: results and open problems[END_REF], that x eq in the interior of S K is a repeller and that the trajectory of x spirals away x eq towards the boundary of S K , without ever reaching it in finite time (a quasi-fixation type event as defined in Kimura 1954).

Examples. (i) Let

A be a K × K skew-symmetric matrix. Let us show that for K odd or even, the conditions under which there is a positive state x for which Ax = 0, is amenable to the sign structure of some determinants by Cramer's rule.

-The odd case. Let K be odd. Suppose A has rank K -1. So 0 is an algebraically and geometrically simple eigenvalue of A. We want to check under what conditions Ax = 0 for some unique x 0 with |x| = 1. We can decompose A as follows, where a is a size K -1 column vector and rank(A) = K -1,: -The even case. Let K be even. Suppose A has rank K -2 so that 0 is a double eigenvalue of A. Then, up to a permutation of the rows and columns, A can be written as:

A = A a -a 0 . If Az =
A =   A a b -a 0 -a K -b a K 0   ,
with (i) rank(A) = K -2 and (ii) there exist c and d of dimension K -2 such that:

  A -a -b   c =   a 0 a K   and   A -a -b   d =   b -a K 0   .
Thus, Ac = a defines uniquely c), a c = 0 because A is skew-symmetric and

-b c = a K fixing a K ). Similarly, A d = b defines uniquely d, b d = 0 because A is skew-symmetric and a d = -b c = a K .
Under these conditions, the vectors c = [c, -1, 0] and d = d, 0, -1 obey Ac = 0 and Ad = 0; so they are the two linearly independent vectors of the kernel of A. If both c ≺ 0 and d ≺ 0, then [-c, 1, 0] and -d, 0, 1 are two nonnegative linearly independent vectors of the kernel of A. Then A can be written as:

A =    A b 1 b 2 -b 1 0 0 -b 2 0 0    , The sub-matrix A := A b 1 -b 1 0
is amenable to the previous study with K replaced by K -1; so it has rank K -2.

Asymmetric permutations.

Let us note by • the Hadamard product of matrices. Let S be a permutation matrix defined by a permutation σ that has no fixed points nor 2-cycles. Then S • S = 0 and S is called a strictly asymmetric permutation matrix. Let Q = (S -S )/2 and z 0. Define A = D -1 z QD -1 z . We have A = -A and, in view of Q1 = 0 and 1 Q = 0, we find Az = D -1 z Q1 = 0 and z A = 1 QD -1 z = 0.

Let K be odd. Let us define the following sequence (S l : l = 1, .., (K -1)/2) of strictly asymmetric permutation matrices. Let S 1 be one of these matrices, then the matrix

Q 1 := S 1 -S 1 has 2K entries ±1 and A 1 := D -1 z Q 1 D -1 z
has the desired properties. Now, for each k = 2, .., (K -1)/2 we take a strictly asymmetric permutation matrix

S k such that S k • (S k-1 + S k-1 ) = 0. Define S k := S k-1 + S k , so S k • S k = 0. Then Q k := S k -S k has 2kK entries ±1 and A k := D -1 z Q k D -1
z has the same properties as A. The matrix

Q (K-1)/2 := S (K-1)/2 -S (K-1)/2 has K(K -1) entries ±1 and A (K-1)/2 := D -1 z Q (K-1)/2 D -1
z has the desired properties but it has no zero entries outside the diagonal. This shows how to generate skewsymmetric matrices Q with entries Q k,l ∈ {-1, 1} for all k = l, iteratively nesting 'mutually orthogonal' strictly asymmetric permutation matrices.

Nonsymmetric Nonnegative Matrices sharing same Perron eigenvector

Let B be a nonnegative nonsymmetric irreducible matrix such that B and B share the same Perron eigenvector z 0, so with Bz = B z = ρz, where ρ > 0 is the real Perron eigenvalue (the spectral radius of B, see [START_REF] Horn | Matrix analysis[END_REF]. Then, for the skew-symmetric part

A := (B -B )/2 of B, Az = 0. For λ > 0 define W = J/ |z| + λ -1 A.
There is a value λ c such that W 0 for λ ≥ λ c and we have

W z = 1 + λ -1 Az = 1.
Let K be odd. We will construct a class of matrices B satisfying the above properties by using a perturbation analysis of the symmetric case. Let C 0 be a symmetric irreducible nonnegative matrix and z =: z 0 0 be the right and left Perron eigenvector associated to the Perron eigenvalue ρ > 0. We can assume z 0 z 0 = 1.

Let z i , i = 1, ...K -1 be such that z i , i = 0, ...K -1 constitute an orthonormal system, so z i z j = δ i,j . Let D = K-1 i=0 γ i z i z i , where P i = z i z i are projectors with P 2 i = P i and choose the real numbers γ i in such a way that γ 0 = 1 and D 0 (this is always possible because D = z 0 z 0 + K-1 i=1 γ i z i z i and z 0 z 0 is a matrix with strictly positive entries). Clearly D is symmetric with

D 2 = DD = K-1 i=0 γ 2 i z i z i . Then B = CD is

4.1.

Characterization in terms of stochastic matrices. Let K > 2. It turns out that the class of K × K skew-symmetric matrices A for which there is a unique z 0 such that Az = 0 admits a representation in terms of some stochastic irreducible matrices. We note that the matrix A has rank K -

1 if K is odd or K -2 if K is even.
In the next result we shall consider P a stochastic irreducible matrix. Then, it has a unique invariant probability vector π. We recall that P obeys detailed balance if π k P k,l = π l P l,k for all k, l ∈ I K . This is equivalent to the fact that P is self-adjoint with respect to the inner product defined by π. Theorem 2. Let K > 2.

(i) Let P be a stochastic irreducible matrix not obeying detailed balance. Let π be its invariant row vector, so π = π P . Let z 0 and define B = D -1 z D π P D -1 z 0. Then D π P is an irreducible nonnegative matrix whose row and column sum vectors coincide. Furthermore, Bz =B z and A = (B -B )/2 is skew-symmetric such that Az = 0.

(ii) Any matrix A with the latter property can be represented as (B -B ) /2, where

B = D -1 z D π P D -1 z .
Proof: (i) Bz = B z = D -1 z π and so Az = 0. We have (D π P ) 1 = (D π P ) 1 = π. When P obey detailed balance, π k P k,l = π l P l,k and so B k,l = B l,k leading to a degenerate A ≡ 0. The one-dimensional simplex for which case K = 2 is ruled out because in this case P always obeys detailed balance.

(ii) Let A be a skew-symmetric matrix such that Az = 0 for some z 0. Assume rank(A) = K -1, which entails that K is odd. Let w 0 and S = (ww )/w z a rank 1 symmetric matrix.

= ρ -1 D -1 z BD z . Then P 1 = ρ -1 D -1 z Bz = 1 so P is stochastic. We have π = π P , where π = D 2 z 1. 2
Let K be even, A = (B -B )/2 and S = (B + B )/2. Assume rank(A) = K -2. Ker(A) consists of two complex conjugate linearly independent and orthogonal eigenvectors z 1 and z 2 = z 1 (with

z 1 • z 2 := k z 1,k z 2,k = z 2 • z 1
) and z = ξz 1 + ξz 2 = 2 Re(ξz 1 ) 0 for some complex number ξ with complex conjugate ξ. Then B = S + A with rank(B) = K and Bz = B z = Sz = w. Note that when we take a realization of a random stochastic matrix P obtained while considering that its rows are drawn independently from Dirichlet(1, .., 1) distributions. With probability 1, P 0 and it is irreducible with invariant measure π. Let B = D -1 z D π P D -1 z ; then 0 is a double eigenvalue of A = (B -B )/2 with eigenvectors z 1 and z 2 = z 1 generating Ker(A) in C and z can be chosen as z = ξz 1 + ξz 2 = 2 Re(ξz 1 ) 0 for some ξ. 2 4.2. Stability of polymorphic states. Let Λ be a nonsymmetric bi-substochastic matrix, that is Λ = Λ , Λ 0 with Λ1 := q < 1 and Λ 1 := r < 1. Let λ > 0. Then, the matrices B = λ -1 (I -Λ) -1 and B = λ -1 (I -Λ ) -1 are potential matrices satisfying B = B . Denote by z and v the solutions of Bz = 1 and B v = 1, which are nonnegative. Since v Bz = v 1 = 1 z, we get |z| = |v|. Assume q = r. Then, z = v and z = B -1 1 = λ(1 -q) . Let A := (B -B )/2 and z = λ(1 -q). Clearly Az = (Bz -B z) /2 = 0 but also Az 0 = 0, where z 0 := 1 -q.

Let W := |z 0 | -1 J + A. It has nonnegative entries if and only if λ ≥ λ c = |z 0 | max k,l (I -Λ) -1
k,l /2), it satisfies W i,j +W j,i = 2/ |z 0 | and W has a positive right equilibrium potential obeying W z 0 = 1.

Consider the following dynamics over S

K : p(x) = D x W x/ω(x) with ω(x) := x W x. Since x Ax = 0 we get x W x = 1/|z 0 |. This is also p(x) = |z 0 |D x W x = |z 0 |D W x x = x+|z 0 |D x Ax or ∆x = |z 0 |D x Ax , Note that ω(x) = x W x is a constant of motion because W is a flat fitness matrix.
This dynamics has a polymorphic equilibrium state x eq = z 0 / |z 0 |. The stability condition of x eq results from the study of the eigenvalues of the Jacobian matrix J eq := ∂p(x) ∂x xeq at equilibrium, acting to the right on zero sum vectors δ = x-x eq , namely

J eq = I + 1 ω (x eq ) D xeq W = I + D z0 J |z 0 | + A = I + x eq 1 + D z0 A.
On the zero-sum eigenspace, J eq now reduces to J eq = I + D z0 A. The equilibrium state x eq is stable if the spectral radius of J eq is < 1 and unstable otherwise. Since the diagonal terms of A vanish, the centers of the Gerschgorin disks of J eq are located at point (1, 0) of the complex plane, so ρ(J eq ) > 1 and x eq is unstable (a result due to Liberman 1990).

4.3. GUM matrices. We will introduce GUM matrices and show that the class of GUM nonsymmetric matrices W for which the equilibrium potential vectors of W and W are strictly positive and coincide, are necessarily symmetric. Hence, this class of potential matrices do not admit nontrivial skew-symmetric matrices supplied by Theorem 2.

For the definition and main properties of GUM matrices see [START_REF] Mcdonald | Inverse M-Matrix Inequalities and Generalized Ultrametric Matrices[END_REF][START_REF] Nabben | Generalized ultrametric matrices -a class of inverse M-matrices[END_REF][START_REF] Dellacherie | Description of the sub-Markov kernel associated to generalized ultrametric matrices. An algorithmic approach[END_REF]. Here we give the GUM matrix in its nested block form (NBF) matrices rather than in terms of inequalities as in the case of ultrametric matrices, but such a characterization exists. The matrix W = (W i,j : i, j ∈ I K ) is a GUM matrix in its NBF if there exists a non-trivial partition of the index set I K = I 1 ∪ I and such that W 1 and W 2 admit also this representation in NBF. This algorithmic way of defining a nested block form matrix generates a family of nested block form submatrices, we simple call it a family of nested block form matrices defining W .

The GUM matrix W is nonsingular if no row and no column is equal to 0, and no two rows nor two columns are equal. We will assume that this is the case. On the other hand it was shown in [START_REF] Mcdonald | Inverse M-Matrix Inequalities and Generalized Ultrametric Matrices[END_REF]Nabben and[START_REF] Nabben | Generalized ultrametric matrices -a class of inverse M-matrices[END_REF] that a nonsingular GUM matrix W is a potential matrix, so it has right and left nonnegative equilibrium potential vectors denoted respectively by z and v.

Let us mention that a GUM matrix is symmetric if and only if it is ultrametric. So, for dimension 1 a nonsingular GUM matrix is ultrametric and z = v = α -1 W . If W is a GUM matrix has a constant row and a constant column (we can assume that it is the first row and the first column), both with constant value α W then z = v = α -1 W e 1 , where e 1 is the vector with a 1 in coordinate 1 and all the other coordinates are 0. Moreover in the decomposition ( 16), we have that W 1 is a symmetric matrix whose off-diagonal matrix is the constant α W ; we also have β W = α W and W 2 is a GUM matrix (necessarily with β W 2 ≥ β W and α W 2 ≥ α W ).

Hence, a necessary and sufficient condition in order that the right and left z and v equilibrium potential vectors are (strictly) positive is that every matrix belonging to a family of nested block form matrices defining W has no constant row and constant column equal to its minimum value. Theorem 4. Let W be a nonsingular GUM matrix. Then, if z 0 and v 0 we have z = v if and only if W is symmetric (so W is an ultrametric matrix). We have

W z = W 1 z 1 + α|z 2 |1 I 1 β|z 1 |1 I 1 + W 2 z 2 .
and

z W = z 1 W 1 + β|z 2 |1 I 1 , α|z 1 |1 I 2 + z 2 W 2 .
From W z = 1 and z W = 1 we find

W 1 z 1 = (1 -α|z 2 |)1 I 1 , W 2 z 2 = (1 -β|z 1 |)1 I 1 ; (17) z 1 W 1 = (1 -β|z 2 |)1 I 1 , z 2 W 2 = (1 -α|z 1 |)1 I 2 .
Therefore

z 1 = (1 -α|z 2 |)z W 1 , z 2 = (1 -β|z 1 |)z W 2 z 1 = (1 -β|z 2 |)z W 1 , z 2 = (1 -α|z 1 |)z W 2 .
Then, from α ≤ β we get α = β. We also deduce [START_REF] Liberman | On the relation between the instability of ESS in discrete dynamics and segregation distortion[END_REF] z

W 1 = v W 1 and z W 2 = v W 2 .
Moreover, z 1 0 and z 2 0 also implies z W 1 = v W 1 0 and z W 2 = v W 2 0. A recurrence argument gives α W k = β W k for every matrix W k belonging to a family of nested block form matrices defining W ; so W is symmetric.

We note that we can make a perturbation of a strictly ultrametric matrix to get a nonsymmetric potential matrix such that it shares its Perron eigenvector with its transpose. But the above result implies that this matrix is not a GUM matrix. 2

-a b= a 0 .

 0 0 with z 0, and x := z/|z|, we have Ax = 0 for some x 0 with |x| = 1. Since rank(A) = rank(A) = K -1, there exists a vector b of size K -1 such that A Then b is uniquely given by inversion of A b = a. We have a b = 0 because b a = b A b = 0, in view of A being skew-symmetric. The vector b := b -1 solves Ab = 0. If b < 0, then z := -b is a positive vector of the kernel of A. Whether b ≺ 0 or not can be checked while using Cramer's rules to compute b = A -1 a. If b ≺ 0, x = z/|z| defines a vector in the one-simplex, solution to Ax = 0.For instance the matrix A = has the required property if ac > 0 and bc < 0 (or ac > 0 and ab < 0. Let W = J/|z| + λ -1 A, where λ > 0. There is a value λ c of λ such that W 0 for λ ≥ λ c and we have W z = 1.

  nonnegative and asymmetric because [C, D] = B -B = 0. Note also that Dz 0 = z 0 = D z 0 and Dz i = γ i z i i = 1, ..., K -1, so that z i (and z i ) is the right (left) eigenvector of D with eigenvalue γ i . We have Bz = Bz 0 = CDz 0 = Cz 0 = ρz 0 and z B = z 0 B = z 0 CD = ρz 0 D = ρz 0 . So B is an asymmetric nonnegative matrix with the same left and right Perron eigenvector z =: z 0 0, associated to the Perron root ρ of C.

2 π 1 ,and

 21 Let B = S +A with rank(B) = K. Then, A = (B -B )/2, where B is an irreducible nonnegative matrix such that Bz = B z = w. Further, C := D z BD z 0 is an irreducible nonnegative matrix whose row and column sum vectors coincide with C1 = D z Bz = C 1. Let c = C1 = C 1. Then P = D -1 c C is an irreducible stochastic matrix with invariant measure π = c. Finally, one can choose w =D -1 z π, fully characterizing B. Because B is full rank, so is C and then also P . 2 As an example consider P a bi-stochastic nonsymmetric matrix. The invariant vector is π = 1 and so D π is the identity. Therefore the representation of the above Theorem gives A = D -1 z QD -1 z 0, where Q = (P -P )/2. Corollary 3. (i) If z = π, B = P D -1 π and Bz = B z = 1. Then A = (B -B ) /2 is skew-symmetric such that Az = 0. (ii) If z = D 1/Bz = B z = z so z is the common left and right eigenvector of B associated to the eigenvalue 1. Then A = (B -B ) /2 is skew-symmetric such that Az = 0. (iii) Let ρ > 0. Then B = ρD matrix with z as the common left and right Perron eigenvector associated to the eigenvalue ρ. Any matrix B with this latter property can be represented in this way. Proof: Only (iii) deserves a proof. Let B 0 be an irreducible nonnegative matrix sharing z 0 as a common left and right Perron eigenvector. Define P

Proof:= z 1 z 2 ,

 2 Assume z = v, we must show W is symmetric. Let us make the partition z with z 1 of dimension |I 1 | and z 2 of dimension |I 2 |. Denote α = α W and β = β W .

  2, are again ultrametric.matrices. If |I 1 | > 1, by permuting the indices in I 1 if necessary, we now have

  W , 1 I 1 1 I 2 β W 1 I 2 1 I 1 W 2

		2 , such that
	(16) α with W = W 1
	α W = min i,j∈I K	W ij and β W ≤ min{ min i,j∈I 1 ,i≥j	W ij , min i,j∈I 2 ,i≥j	W i,j };

Acknowledgments:

T. Huillet and S. Martínez acknowledge support from the Project Basal PFB 03 of the CONICYT, Chile. S. Martínez thanks the University of Cergy-Pontoise for a visiting Professor opportunity. T. Huillet also acknowledges support from the "Chaire Modélisation mathématique et biodiversité " and from the labex MME-DII Center of Excellence (Modèles mathématiques et économiques de la dynamique, de l'incertitude et des interactions, ANR-11-LABX-0023-01 project).