B Emek Abali

Wolfgang H Miiller

Francesco Dell'isola

B E Abali
email: b.e.abali@berkeley.edu

In continuum mechanics, there exists a unique theory for elasticity, which includes the first gradient of displacement. Tue corresponding generalization of elasticity is referred to as strain gradient elasticity or higher gradient theories, where the second and higher gradients of displacement are involved. Unfortunately, there is a lack of consensus among scientists how to achieve the generalization. Various suggestions were made, in order to compare or even verify these, we need a generic computational tool. In this paper, we follow an unusual but quite convenient way of formulation based on action principles. First, in order to present its benefits, we start with the action principle leading to the well-known form of elasticity theory and present a variational formulation in order to obtain a weak form. Second, we generalize elasticity and point out, in which term the suggested formalism differs. By using the same approach, we obtain a weak form for strain gradient elasticity. The weak forms for elasticity and for strain gradient elasticity are solved numerically by using open-source packages-by using the finite element method in space and finite difference method in time. We present some applications from elasticity as well as strain gradient elasticity and simulate the so-called size effect.

Introduction

Starting with works Joseph-Louis de Lagrange, Pierre Louis Moreau de Maupertius, Leonhard Euler-for a short history see Romano et al. [38)-the theory of mechanics is formulated by an energy expression called LAGRANGEan that depends on primitive variables and their first derivatives in coordinates. Such a theory is of first order since the first derivative in coordinates describes the immediate neighborhood. This is an infnitesimally small region given by a sphere of radius s. In continuum mechanics, we declare the size of s as approaching zero. Therefore, from a conceptual point of view, no problerns arise and continuum mechanics is sufficient to characterize the immediate neighborhood with the first derivative. This approach is also called the locality argumentation and is valid for analytical mechanics, where we can set s infnitesimally small. However, if we want to exploit numerica} mechanics, then we have to set a finite value for s. It has to be smaller than the length-scale of the underlying structure.

In order to explain the length-scale, consider a measurement where the machine precision is limited to 1 μm. Therefore, we simply miss any deformation occurring on a smaller scale than a micrometer length-scale. We have to choose ε smaller than 1 μm, and then we can accurately model the measurement. Moreover, we set ε large enough such that sufficiently many atoms are included, so we may ignore any fluctuation effects. Suppose we can estimate 1 nm as the minimum size of ε, i.e., for a monatomic lattice we assume that first neighbor is in a distance less than 1 nm. In a smaller length-scale, the atomistic structure and interactions between lattice points become important. Between 1 nm and 1 μm, the size of ε can be chosen arbitrarily, and the solution of the problem remains the same; for the idea of self-similarity and intermediate distances, see Barenblatt [10,Chap. 2]. The first-order theory holds accurately for a material point in the given length-scale, i.e., ε is smaller than the upper limit of 1 μm and greater than the lower limit of 1 nm. We can measure with a measurement device capturing 1 μm and model the measurement with the first-order theory. Now consider a material with 1 μm grain size. In other words, the underlying material has a substructure within the length-scale. Now we have two different options to resolve this issue. Either we redeclare the lengthscale by decreasing the upper limit, simply stated, we buy a new measurement device with a higher precision; or we generalize the first-order theory in order to capture the effects of the substructure. This generalization introduces new coefficients responsible for the response effected by the substructure. In this work, we explain the method of least action and exploit it for generalizing the continuum theory from the first order to the second order.

Although the idea of an extended or generalized mechanics is older than a century, see dell'Isola et al. [14], the modern theory started to be developed by the early 1960s with Mindlin and Tiersten [START_REF] Mindlin | Effects of couple-stresses in linear elasticity[END_REF], Mindlin [START_REF] Mindlin | Micro-structure in linear elasticity[END_REF], Toupin [START_REF] Toupin | Theories of elasticity with couple-stress[END_REF], Eringen [18]. Today there are various theories describing an extension of the first-order theory for elasticity. Interestingly, these different theories are mathematically quite similar and can be viewed as specialized versions of a unified theory, see Neff et al. [START_REF] Neff | A unifying perspective: the relaxed linear micromorphic continuum[END_REF]. By starting with experimental evidence as in Morrison [START_REF] Morrison | The yield of mild steel with particular reference to the effect of size of specimen[END_REF] and then later on in Yang and Lakes [START_REF] Yang | Experimental study of micropolar and couple stress elasticity in compact bone in bending[END_REF], Brezny and Green [12], Chen and Fleck [13], Kesler and Gibson [START_REF] Kesler | Size effects in metallic foam core sandwich beams[END_REF], we know that a second-order theory-called the strain gradient theory-is necessary for an accurate prediction of experiments. For example, the prominent experiment of beam bending can be investigated. Consider a beam of 1 μm length and another beam of 1 mm, both out of the same material. The first-order theory is accurate for the beam in the meter length-scale; however, it is inaccurate for the beam in the micrometer length-scale. Apparently, for the micrometer beam, the substructure of the material becomes dominant. An internal substructure can be a crystallographic structure of aluminum or polymer chains of epoxy. If one wants to acquire accurate results with the first-order theory, then a precise and detailed model of this substructure needs to be involved. However, this method is not efficient at all. In other words, although we might try to model everything by molecular dynamics, we never do that for a structure bigger than a couple of micrometers due to the computational time required for such a detailed model. We will present a second-order theory, capable of modeling the deformation of beam bending even on smaller length-scales, with an adequate amount of computational cost.

Computational issues of generalized mechanics are heavily discussed in the literature. There are various propositions and numerical strategies. For example, Bilotta et al. [11] use a promising method by constructing special elements. The newly introduced material parameters can be motivated with a rigorous micro-macro identification procedure, see Pideri and Seppecher [START_REF] Pideri | A homogenization result for elastic material reinforced periodically with high rigidity elastic fibres[END_REF]. Based on this idea, Giorgio [19] provides a persuasive series of numerical results, starting from the first-order theory. In various applications, effects of substructure are unavoidable, for simulations of interesting examples, see Scerrato et al. [START_REF] Scerrato | On the effect of shear stiffness on the plane deformation of linear second gradient pantographic sheets[END_REF], Turco et al. [START_REF] Turco | Large deformations induced in planar pantographic sheets by loads applied on fibers: experimental validation of a discrete lagrangian model[END_REF], and Placidi et al. [START_REF] Placidi | Gedanken experiments for the determination of two-dimensional linear second gradient elasticity coefficients[END_REF]. Finite element implementation and comparable results to the underlying work have appeared recently in Reiher et al. [START_REF] Reiher | Finite-element analysis of polyhedra under point and line forces in second-strain gradient elasticity[END_REF]. We present our formulation and its implementation in a quite general setting; hence, an extension to many other branches than elasticity is possible. For example, the numerical solution procedure presented herein seems to be potentially very useful also for studying the equilibrium problems of second gradient fluids as discussed in Eremeyev [16], Eremeyev and Altenbach [17].

First, we discuss the variational formulation in the first-order theory by means of an action principle. Instead of the least action principle, another similar theory known as principle of virtual work (or power) can be used, too. We start by postulating an action and omit to give its detailed form. If a principle of virtual work needs to be used, then we immediately have to describe the form. There are no big discrepancies between these methods. Second, we generalize the variational formulation for the second-order theory by following Abali and Müller [3]. The outcome is a weak form necessary for the computation. Third, we present a numerical tool solving the weak form. For the case of elasticity and strain gradient elasticity, we perform several simulations in order to justify the benefits of the approach. We use open-source packages and make our codes publicly available in order to encourage further studies.

Variational formulation for the first-order theory

The general formulation relies on a mathematical construct called action. The first postulate in this methodology is an existing Lagrangean density, L, describing the underlying system

L = L x μ , φ A , φ A,μ , φ A,μ = ∂φ A ∂ x μ , (1)
which depends on Cartesian coordinates, x μ ∈ , μ = 1, 2, . . . m, in m-th-dimensional space, ⊂ R m ; on primitive variables, φ A = φ A (x μ ); as well as on their derivatives, φ A,μ = φ A,μ (x μ ). Only the first derivatives are included in the list of arguments. The Lagrangean density describes the state of the underlying system. Thus, we call {φ A , φ A,μ } the state space. Consider a continuum body defined in the space x μ = (t, x, y, z) is being deformed under a given mechanical loading. Then the state is the deformation of the body defined by the displacement field, φ A = u i = u i (t, x, y, z), and its first gradient in coordinates. The space spanned by the displacement field and its first gradient in coordinates creates the state space. This concept is often used in terms of rigid bodies, where the displacement and velocity (derivative in t) create the state space.

The second postulate is the so-called action functional:

A = L d + ∂ W s d , d = dx 1 dx 2 . . . dx m , d = dx 1 dx 2 . . . dx m-1 , ( 2
)
where we have used the Cartesian coordinates such that the metric determinant of the space is equal to 1.

If we call d an infinitesimal volume element, then one-dimension less d would be a surface element. We have a first-order theory such that only one-dimension lower space is included. Technically, we know W s on the surface and want to determine the appropriate L within the space. Thus, W s is given and we restrict it be depending on φ A . The action is a scalar, i.e., it is invariant with respect to the transformation of primitive variables. First, we introduce variations of the primitive variables,

φ A = φ A + εδφ A , (3)
and of its rates,

φ A,μ = φ A,μ + εδφ A,μ , δφ A,μ = (δφ A ) ,μ = ∂δφ A ∂ x μ . ( 4
)
The test functions, δφ A , are arbitrary; however, they are chosen to vanish on the Dirichlet boundaries, ∂ D , where the primitive variables are known. The constant factor ε is such a small number (with respect to the upper limit in the length-scale) that the above transformation is accurate in the first order in ε. Now, with the same accuracy, the variation of the action reads

δA = L d - L d + ∂ W s d - ∂ W s d , L = L x μ , φ A , φ A,μ , W s = W s x μ , φ A .
(5) By using a Taylor expansion linear in ε, we obtain

L = L x μ , φ A , ∂φ A ∂ x μ = L x μ , φ A , φ A,μ + ∂L x μ , φ A , φ A,μ ∂φ A (φ A -φ A ) + ∂L x μ , φ A , φ A,μ ∂φ A,μ (φ A,μ -φ A,μ ) = L + ∂L ∂φ A εδφ A + ∂L ∂φ A,μ εδφ A,μ , (6)
and analogously,

W s = W s + ∂ W s ∂φ A εδφ A , (7)
where and henceforth over repeated indices we apply the Einstein summation convention. The variation of action becomes

δA = ∂L ∂φ A εδφ A + ∂L ∂φ A,μ εδφ A,μ d + ∂ ∂ W s ∂φ A εδφ A d . ( 8
)
According to the principle of least action, the action functional is a scalar such that its variation vanishes

δA = 0 , (9)
for non-dissipative (reversible) systems. Since ε is a constant nonzero number, we can divide by it such that we obtain

∂L ∂φ A δφ A + ∂L ∂φ A,μ δφ A,μ d + ∂ ∂ W s ∂φ A δφ A d = 0 . ( 10
)
This integral form is general for a system having a Lagrangean density depending on the primitive variables and their first derivatives. Often, an integration by parts is established on the second term within such that the integrand contains only δφ A . This step introduces a surface integral derived from L set in with ∂ W s /∂φ A . Within the so-called Euler-Lagrange equations are obtained; over the surface, a boundary condition is defined. For enabling this relation, we have restricted W s depending only on φ A and not on φ A,μ . Although we omit a direct use of this relation for generating the weak form, it serves an understanding of a boundary term and its specifications in an experimental setting.

Generating the weak form

Consider a material system composed of a three-dimensional continuum body in a reference frame, B 0 ⊂ R 3 , with its closure, ∂B 0 , such that B = B 0 ∪ ∂B 0 . Material particles of the continuum body are denoted by their positions, X, in the reference frame. The reference frame is an inertial frame, and we choose it as fixed. In the beginning of the simulation, where the continuum body is undeformed, we set the initial frame as the reference frame and indicate this choice with the use of an index "0." Material particles' positions are expressed in a Cartesian coordinate system, X i ∈ B with i = 1, 2, 3, in the initial frame, and they form an infinitesimal volume element dV = dX 1 dX 2 dX 3 as well as an infinitesimal surface element d A. For this configuration, the space is four-dimensional, x μ = (t, X i ), with time t ∈ τ = [t 0 , t 1 ] and physical space X i such that the integral form in Eq. ( 10) reads

τ B 0 ∂L ∂φ A δφ A + ∂L ∂φ • A δφ • A + ∂L ∂φ A,i δφ A,i dt dV + τ ∂B 0 ∂ W s ∂φ A δφ A dt d A = 0 , (11)
where (•) • denotes a partial time derivative since the reference frame is fixed. The integral form is continuous.

In a computational environment, we will solve it for discrete time steps in discrete elements (subsets of the physical space). First, we discretize in time. For the sake of a simplified notation (and programming), we choose constant time steps such that the time becomes the following list:

t = {0, t, 2 t, 3 t, . . . , t end } (12)
As a consequence, the discrete (in time) representation of φ A is a list,

φ A = {φ A (X i , t = 0), φ A (X i , t = t), φ A (X i , t = 2 t), . . . , φ A (X i , t = t end )} . ( 13
)
The time rate of the primitive variables read

φ • A = ∂φ A ∂t = φ A -φ 0 A t , ( 14
)
where φ A denotes the unknown value to be computed and φ A indicates the already computed value from the last time step. This formulation converges to the analytic time rate for t → 0. Since the test functions δφ A are arbitrary, we may set them as constants in the beginning of a simulation. For doing so, we have to eliminate δφ • A from the integral form by integrating by parts

τ B 0 ∂L ∂φ A δφ A - ∂L ∂φ • A • δφ A + ∂L ∂φ A,i δφ A,i dV dt + B 0 ∂L ∂φ • A δφ A dV τ + τ ∂B 0 ∂ W s ∂φ A δφ A dt d A = 0 . ( 15
)
By given values of φ A at two subsequent time instants like t 0 and t 1 , we set the test functions zero on τ = [t 0 , t 1 ] and obtain for every other time instant

B 0 ∂L ∂φ A δφ A - ∂L ∂φ • A • δφ A + ∂L ∂φ A,i δφ A,i dV + ∂B 0 ∂ W s ∂φ A δφ A d A = 0 . ( 16
)
Second, we discretize the physical space. The primitive variables as well as the test functions are defined on a Hilbert space with continuity k in a discrete element, E, with the Dirichlet boundary condition on ∂ E D ,

V = φ A , δφ A ∈ [H k (E)] 3 : φ A = given, δφ A = 0 on ∂ E D . ( 17
)
It is technically a Sobolev space where also the differentiability within the domain is ensured. The continuity k shall be as high as necessary. It depends on the chosen Lagrangean density, which has a derivation in φ A,i .

In the Euler-Lagrange equations, one more integration by parts is used to "switch" the derivative in the test function to that term. Hence, the regularity condition is stronger than the latter form. Therefore, the latter is called a weak form to be fulfilled in each element, E, finite in size. The totality of all finite elements reads the following weak form:

Form = nr of ele. n=1 E n ∂L ∂φ A δφ A - ∂L ∂φ • A • δφ A + ∂L ∂φ A,i δφ A,i dV + ∂B 0 ∂ W s ∂φ A δφ A d A . ( 18
)
The external work done on the surface is given, W s = W s (t, X i , φ A ), and the surface element d A denotes an integration over the surface belonging to boundary ∂B 0 excluding Dirichlet boundaries. On the Dirichlet boundaries, the numerical value of φ A is given, since we know them δφ A on the Dirichlet boundaries vanishes, called Galerkin approach. The weak form is a general integral form applicable to all reversible systems defined by a Lagrangean depending on primitive variables and their derivatives.

Elasticity theory

In elasticity the primitive variable is the displacement:

φ A = u i , u i = x i -X i , ( 19
)
where the displacement is the deviation from the reference (initial) position, X i , to the current position, x i . The Lagrangean density for the elasticity reads

L = 1 2 ρ 0 u • i u • i -w + ρ 0 f i u i , ( 20
)
where the first term denotes the kinetic energy, the second term indicates the stored or deformation energy, and the third term is called the potential energy. In the case of elasticity, the stored energy density (energy per volume) w depends on displacement and its first derivative in space. Its exact form relies on the underlying material. The conservative force f i is simply the gravitational specific force (force per mass). We present two different material models in the following sections, and they only differ in the definition of the stored energy density, w. On Neumann boundaries, the energy expression

W s = ti u i ( 21
)
implements a traction vector ti in the unit of stress, modeling a mechanical loading on a part of the surface. By inserting the latter in Eq. ( 18), we obtain the weak form discrete in time and space, Form = nr of ele.

n=1 E n f i δu i -ρ 0 u i -2u 0 i + u 00 i t t δu i - ∂w ∂u i, j δu i, j dV + ∂B 0 ti δu i d A . ( 22
)
This weak form can be solved by using k = 1 since we have only first-order derivative in u i . In other words, linear continuous standard elements will be used. Indeed, the latter weak form can be solved for various definitions of the stored energy density. It is of importance to note that the term:

∂w ∂u i, j , (23)
is called Piola stress:

P ji = ∂w ∂ F i j , ( 24
)
since the deformation gradient, F ji , leads to

F ji = ∂ x j ∂ X i = δ ji + u j,i , ∂w ∂u i, j = ∂w ∂ F kl ∂ F kl ∂u i, j = P lk ∂u k,l ∂u i, j = P lk δ ki δ l j = P ji . ( 25
)
The polynomial degree of Eq. ( 23) states the degree of the material model. The stored energy density is one degree higher, for example, for a linear material model w is quadratic. By defining the Lagrangean density in Eq. ( 20), we have introduced several assumptions. First, we have assumed that it be independent of time. In other words, f i is a conservative force and does not depend on time explicitly. Second, we assume that the kinetic energy depends on u • i and the stored energy depends on u i and u i, j . In other words, we exclude viscoelasticity by neglecting a dependence on u • i, j . Third, we use the deformation gradient or the Green-Lagrange strain measure

E i j = 1 2 u k,i u k, j + u (i, j) , ( 26
)
instead of displacement dependence in the stored energy. So the stored energy is

w = w(X i , F i j ) = w(X i , E i j ) . ( 27
)
As a special case, for homogeneous materials the dependency on X i vanishes. Due to the principle of objectivity, the scalar function depends on the invariants of the deformation gradient or strain.

Variational formulation for the second-order theory

In order to generalize the proposed method, we redefine the Lagrangean density as including the second derivatives of the primitive variables

L = L(x μ , φ A , φ A,μ , φ A,μν ) , φ A,μ = ∂φ A ∂ x μ , φ A,μν = ∂ 2 φ A ∂ x μ ∂ x ν . ( 28
)
In this configuration, the state space reads

{φ A , φ A,μ , φ A,μν } (29)
with the following variations:

φ A = φ A + εδφ A , φ A,μ = φ A,μ + εδφ A,μ , φ A,μν = φ A,μν + εδφ A,μν . ( 30
)
The same action definition and the variation of it lead to

L = L + ∂L ∂φ A εδφ A + ∂L ∂φ A,μ εδφ A,μ + ∂L ∂φ A,μν εδφ A,μν , (31)
and provide the following integral form:

∂L ∂φ A δφ A + ∂L ∂φ A,μ δφ A,μ + ∂L ∂φ A,μν δφ A,μν d + ∂ ∂ W s ∂φ A δφ A d + ∂ ∂ W s ∂φ A,μ δφ A,μ d + ∂∂ ∂ W e ∂φ A δφ A d + ∂∂ ∂ W e ∂φ A,μ δφ A,μ d = 0 , ( 32
)
after using the principle of least action and dividing by ε. By using the same justification as in the last section, we allow W s on surface d and W e on edge d depending only on φ A and φ A,μ but not on φ A,μν .

Generating the weak form

As in the last section, again for x μ = (t, X i ) with t ∈ τ = [t 0 , t 1 ] and X i ∈ B 0 , we obtain

τ B 0 ∂L ∂φ A δφ A + ∂L ∂φ • A δφ • A + ∂L ∂φ A,i δφ A,i + ∂L ∂φ •• A δφ •• A + ∂L ∂φ • A,i δφ • A,i + ∂L ∂φ A,i j δφ A,i j dV dt + τ ∂B 0 ∂ W s ∂φ A δφ A dt d A + τ ∂B 0 ∂ W s ∂φ • A δφ • A + ∂ W s ∂φ A,i δφ A,i dt d A + τ ∂∂B 0 ∂ W e ∂φ A δφ A dt d + τ ∂∂B 0 ∂ W e ∂φ • A δφ • A + ∂ W e ∂φ A,i δφ A,i dt d = 0 . ( 33
)
In the time discrete fashion, we eliminate terms with rates of test functions by integrating by parts. Then we use the fact that the values φ 0 A and φ 00 A are known in the beginning such that δφ A vanishes. After discretization in space, we obtain the following weak form Form = nr of ele.

n=1 E n ∂L ∂φ A δφ A - ∂L ∂φ • A • δφ A + ∂L ∂φ A,i δφ A,i + ∂L ∂φ •• A •• δφ A - ∂L ∂φ • A,i • δφ A,i + ∂L ∂φ A,i j δφ A,i j dV + + ∂B 0 ∂ W s ∂φ A - ∂ W s ∂φ • A • δφ A + ∂ W s ∂φ A,i δφ A,i d A + ∂∂B 0 ∂ W e ∂φ A - ∂ W e ∂φ • A • δφ A + ∂ W e ∂φ A,i δφ A,i d . ( 34
)
The extension to the second-order theory introduced many different boundary conditions. All of them have to be given, and they are discussed by many researchers. Mainly we have boundary terms acting on surface and additionally terms acting on edges. Since now not only φ A but also its derivative is involved, we need to approximate the discrete representation of φ A from a Sobolev space,

V = φ A , δφ A ∈ [H k (E)] 3 : φ A = given, δφ A = 0 on ∂ E D , (35)
with k ≥ 2.

Strain gradient elasticity theory

In strain gradient elasticity, the primitive variable is the displacement:

φ A = u i . ( 36
)
For being able to compare the elasticity and strain gradient elasticity, we implement exactly the same surface energy and neglect the edge energy,

W s = ti u i , W e = 0 , (37)
where the traction vector ti denotes again the mechanical load on the Neumann boundary. These boundary conditions are heavily discussed in the literature, and their application to experimental mechanics is very challenging. For some useful interpretations of these surface and edge energies and how they affect the solution, among others, we refer to Auffray et al.

[8], Javili et al. [START_REF] Javili | Geometrically nonlinear higher-gradient elasticity with energetic boundaries[END_REF], Steigmann and dell'Isola [START_REF] Steigmann | Mechanical response of fabric sheets to three-dimensional bending, twisting, and stretching[END_REF]. Although we obtained these additional terms on boundaries, they are challenging to justify. For the specific case of a quadratic energy, by using variational formulation in Mindlin [START_REF] Mindlin | Second gradient of strain and surface-tension in linear elasticity[END_REF], by using principle of virtual power in Polizzotto [START_REF] Polizzotto | A second strain gradient elasticity theory with second velocity gradient inertia-part i: constitutive equations and quasi-static behavior[END_REF], it is shown that W s and W e may include relations leading to a dependency in u i and u i, j . Moreover, the aforementioned relations can be rewritten, see [START_REF] Polizzotto | A second strain gradient elasticity theory with second velocity gradient inertia-part i: constitutive equations and quasi-static behavior[END_REF]Sect. 5], in order to relate them to measurable quantities like surface tension (energy). Such an analysis results in W s depending on u i and surface normal component of gradient of displacement. Herein, we have even proposed to involve dependencies on u • i . Unfortunately, we lack a clear justification of experimental feasibility for such dependencies. For simplicity, we choose Eq. ( 37) to implement such that a comparison with the well-known linear elasticity is possible.

For nonpolar materials, we propose the following Lagrangean density:

L = 1 2 ρ 0 u • i u • i -w + ρ 0 ( f i u i + l i j u j,i ) . ( 38
)
The first term is an inertial term. Herein, we simplify the formulation by assuming that only the velocity and ordinary momentum (mass times velocity) take part in this term, see Polizzotto [START_REF] Polizzotto | A second strain gradient elasticity theory with second velocity gradient inertia-part i: constitutive equations and quasi-static behavior[END_REF], Polizzotto [START_REF] Polizzotto | A second strain gradient elasticity theory with second velocity gradient inertia-part ii: dynamic behavior[END_REF] for an elaborate discussion. So the inertial term consists of the kinetic energy due to the macroscopic velocity for nonpolar materials. For a polar material, this term would have a term with spin, too. The second term is again the stored energy; however, the list of argumentation is amended such that the stored energy depends not only on the first but also on the second gradient of the primitive variable, i.e., displacement,

w = w(X i , E i j , E i j,k ) . ( 39
)
We will restrict the formulation to homogeneous materials such that the dependency on material particles' positions, X i , vanishes. The third term denotes the energy due to the volumetric effects. The first part, f i u i , is the specific energy because of the gravitational forces creating a translational displacement (implied by the linear momentum). The second part, l i j u j,i , is the specific volumetric energy leading to a rotational displacement (implied by the angular momentum). For a nonpolar medium, this term vanishes. By inserting the Lagrangean density, we obtain the weak form for a nonpolar material with traction-type boundaries:

Form = nr of ele. n=1 E n ρ 0 f i δu i -ρ 0 u i -2u 0 i + u 00 i t t δu i - ∂w ∂u i, j δu i, j + ρ 0 l ji δu i, j - 1 t t ∂w ∂u i, j -2 ∂w 0 ∂u 0 i, j + ∂w 00 ∂u 00 i, j δu i, j - ∂w ∂u i, jk δu i, jk dV + ∂B 0 ti δu i d A , ( 40
)
where the following notation is introduced

w 0 = w(E 0 i j , E 0 i j,k ) , w 00 = w(E 00 i j , E 00 i j,k ) , E 0 i j = 1 2 u 0 k,i u 0 k, j + u 0 (i, j) , E 00 i j = 1 2 u 00 k,i u 00 k, j + u 00 (i, j) , (41)
with the Green-Lagrange strain tensor, E i j , and its values in the last and second last time steps, E 0 i j , E 00 i j , respectively. The latter form can be solved by using k = 2 since we have second-order derivative in u i . In other words, we will use quadratic continuous standard elements.

Applications

The weak forms in Eq. ( 22) for elasticity and in Eq. [START_REF] Steigmann | Mechanical response of fabric sheets to three-dimensional bending, twisting, and stretching[END_REF] for strain gradient elasticity are general forms applicable to linear as well as nonlinear materials. Moreover, since we use a nonlinear strain measure, namely the Green-Lagrange strain tensor, the formulation captures geometric nonlinearities as well. Other than elastic response, we have not made any assumptions so far. In this section, we will present a computational tool for solving the presented nonlinear weak forms numerically. By using the Python programming language, see Oliphant [START_REF] Oliphant | Python for scientific computing[END_REF], we develop a code for obtaining solution of weak forms by exploiting the collection of open-source packages developed under the FEniCS project, see Hoffman et al. [START_REF] Hoffman | Fenics[END_REF], Logg et al. [START_REF] Logg | Automated Solution of Differential Equations by the Finite Element Method the FEniCS Book[END_REF], Jones et al. [START_REF] Jones | SciPy: Open source scientific tools for Python[END_REF]. We use finite element method in space and finite difference method in time. For the visualization of results, we use ParaView, see Ahrens et al.

[5], Ayachit [9].
Three examples are shown for linear elastic, hyperelastic, and linear strain gradient elastic materials. All examples have nonlinear weak forms such that the weak form needs to be linearized before solving. This linearization is gained by a symbolic derivation of the weak form, before assembly. Therefore, the linearization is at the level of the partial differential equations. The derivation is fully automatized by using symbolic derivative functionality, see Alnaes and Mardal [6], Alnaes and Mardal [7]. Moreover, we have implemented the code by defining the stored energy function. Its derivation, precisely the solution of Eq. ( 23), is computed by using the aforementioned symbolic derivative functionality. This abstraction is of paramount importance, and the following three examples differ only in the definition of the stored energy function.

For learning the basics of programming in Python for FEniCS ecosystem, see Langanten and Logg [START_REF] Langanten | of Simula SpringerBriefs on Computing[END_REF]. For several engineering applications by using FEniCS, see Abali [1]. In order to encourage further research, we publish all used codes in Abali [2] under the GNU Public license as stated in Gnu Public [START_REF]Gnu general public license[END_REF].

Example for linear elasticity

A structure made of an engineering steel can be modeled with a linear elastic material if the mechanical loading generates stresses less than the yield stress such that no plastic deformation occurs. Consider a thin plate modeled as a three-dimensional structure under a mechanical loading. Under a concentrated loading, plate underlies a large deformation (with respect to the geometric dimensions). For the linear model of a homogeneous steel, the generic form of the quadratic scalar function reads w = E i j C i jkl E kl , [START_REF] Turco | Large deformations induced in planar pantographic sheets by loads applied on fibers: experimental validation of a discrete lagrangian model[END_REF] where C i jkl is called the stiffness or elasticity tensor with 21 independent coefficients since

C i jkl = C jikl , C i jkl = C i jlk , C i jkl = C kli j , ( 43
)
effected by the symmetric strain, E i j = E ji . This case is applicable for any material class. Existing crystal symmetries decrease the number of independent coefficients. For example, in the case of an isotropic material, the representation of this rank four tensor simplifies to

C i
w = λ 1 2 E ii E j j + μE i j E ji , (45)
where μ = μ 1 = μ 2 since strain is symmetric and the so-called Lame parameters, λ, μ, are given by Young's modulus, E, and Poisson's ratio, ν as follows

λ = Eν (1 + ν)(1 -2ν) , μ = E 2(1 + ν) . ( 46
)
The strain measure is nonlinear; thus, the geometric nonlinearities are computed accurately. The stored energy density is quadratic in strain; therefore, the material model in Eq. ( 23) is linear. We simulate a plate of an AISI steel with the following parameters:

ρ 0 = 8.1 × 10 -9 Mg/mm 3 , E = 200 × 10 3 MPa , ν = 0.3 . ( 47
)
The plate is of 1000 × 600 × 20 mm, and its surface lies in X 1 X 2 -plane. The concentrated force is modeled by a Gaussian distribution, increasing linearly in time

t = at exp b(X 1 -X1 ) 2 + b(X 2 -X2 ) 2 , ( 48
)
such that ti = -t N i applies on a small area on top of the plate against the plane normal, N i . The force increases linearly in time. The gravitational forces are implemented, f i = (0, 0, -9850) N/Mg, although the deformation due to them is negligible (small with respect to the deformation caused by the loading). The transient solution of the weak form in Eq. ( 22) provides the displacement field, u(t, X i ), as expected, for the convergence see Appendix A. The response is immediate, and the deformation goes back to the initial value after unloading. For example, the deformation is shown in Fig. 1 at the time instant of maximum force applied on the top of the plate.

  jkl = λδ i j δ kl + μ 1 δ ik δ jl + μ 2 δ il δ jk , (44) see, for example, Abali et al. [4, Appendix] for a derivation. By inserting the stiffness tensor into the energy density, we obtain the so-called St. Venant-Kirchhoff model for the stored energy:


Nicolas Champagnat
email: nicolas.champagnat@sophia.inria.fr

Amaury Lambert
email: amaury.lambert@ens.fr

Adaptive dynamics in logistic branching populations

Keywords: secondary 60J25, 60J70, 60J75, 60J85, 60K35, 92D10, 92D15, 92D25, 92D40

Introduction

The recent biological theory of adaptive dynamics [START_REF] Hofbauer | Adaptive dynamics and evolutionary stability[END_REF][START_REF] Metz | How should we define 'fitness' for general ecological scenarios ?[END_REF] proposes a description of the long term evolution of an asexual population by putting emphasis on the ecological interactions between individuals, in contrast with classical population genetics models which focus on the genetic structure of the population. The basic models are individual-based models in which the population dynamics is precisely described and includes birth, death, competition and mutation. The basic idea of the theory of adaptive dynamics is to try to get insights about the interplay between ecology and evolution by studying the invasion of a single mutant type appearing in a simplified (monotype stable) resident population. The evolution of the population can then be described as a sequence of mutant invasions in the population. If the resident type goes extinct when the mutant type invades (we say that the mutant type fixates), the evolution is described by the so-called 'trait substitution sequence' (TSS) [START_REF] Metz | Adaptive Dynamics, a geometrical study of the consequences of nearly faithful reproduction[END_REF]. This appraoch has revealed powerful to predict the qualitative behaviour of complicated evolutionary dynamics. In particular, it allows to determine the (local) direction of evolution in the space of phenotypic traits (or simply traits) from the individual ecological parameters, and to predict and explain the phenomenon of evolutionary branching [START_REF] Metz | Adaptive Dynamics, a geometrical study of the consequences of nearly faithful reproduction[END_REF], where a population, initially composed of individuals with traits concentrated around a single trait value, divides into two (or more) subpopulations concentrated around distinct trait values that stably coexist because of their ecological interactions. The description of this phenomenon is an important achievement of this theory, as well as the 'canonical equation of adaptive dynamics' [START_REF] Dieckmann | The dynamical theory of coevolution: A derivation from stochastic ecological processes[END_REF], which describes the evolution of the dominant trait of the population as a deterministic "hill-climbing" process on a fitness landscape which depends on the current state of the population (see (1) below).

More formally, as soon as eternal coexistence of two or more types is not permitted by the model, the evolution can be described by the sequences (T n ) n and (V n ) n , where T n is the n-th time where the population becomes monomorphic (i.e. composed of only one type) and V n is the surviving type at time T n . The sequence (V n ) is the above-mentionned TSS. It is possible to prove the convergence of an individual-based model to the TSS under two biologically motivated assumptions [START_REF] Metz | Adaptive Dynamics, a geometrical study of the consequences of nearly faithful reproduction[END_REF][START_REF] Champagnat | A microscopic interpretation for adaptive dynamics trait substitution sequence models[END_REF]. First, the assumption of rare mutations guarantees that, in the timescale of mutations, the widths of time intervals during which the population is polymorphic vanish, so that there is only one type surviving at any time t. To prevent the population from rapidly becoming extinct in the new timescale, one also has to rescale population sizes, thereby making the assumption of large populations.

Subsequently, the TSS is a Markov jump process on the trait space whose semigroup is shown [START_REF] Champagnat | A microscopic interpretation for adaptive dynamics trait substitution sequence models[END_REF] to depend on the invasion fitnesses (as defined in [START_REF] Metz | How should we define 'fitness' for general ecological scenarios ?[END_REF]) f (x, y), x, y ∈ X , where f (x, y) is defined as the expected growth rate of a single individual of type y -the mutant -entering a monomorphic population of type x 'at equilibrium' -the residents. Note that this fitness is not given a priori, but derived from the microscopic model of individual interactions. Because of the assumption of large population, the sign of this fitness determines the possibility of invasion of a mutant type: if f (x, y) < 0, the mutant type y cannot invade a resident population of type x. Thus, evolution proceeds by successive invasions of (only) advantageous mutant types replacing the resident one.

The 'canonical equation of adaptive dynamics' [START_REF] Dieckmann | The dynamical theory of coevolution: A derivation from stochastic ecological processes[END_REF], which describes the evolution of a one-dimensional trait x as the solution of the following ODE, is obtained from the TSS in the limit of small mutations:

dx dt = 1 2 σ(x) 2 µ(x)n(x) ∂ ∂y f (x, x), (1)
where σ(x) 2 stands for the (rescaled) variance of the mutation step law, n(x) for the equilibrium size of a pure x-type population, and f (x, y) for the invasion fitness mentioned above. Note how only advantageous types get fixed (the trait follows the fitness gradient) and how the fitness landscape y → f (x, y) depends on the current state x of the population. However, it is well-known that slightly deleterious types can be fixed by chance in finite populations. This phenomenon is known under the name of genetic drift. Depending on the strength of genetic drift, selection is said to be strong (genetic drift has negligible effects) or weak. In the large population asymptotic from which the TSS of adaptive dynamics is derived, genetic drift has negligible impact compared to the action of selection. Therefore, the fixation of slightly deleterious types cannot be observed. Our goal here is to include genetic drift in the adaptive dynamics models by considering finite populations under weak selection. We continue using the bottom-up approach of adaptive dynamics; that is, model (macroscopic) evolution from (microscopic) populations. In particular, we allow the population sizes to fluctuate randomly through time and we aim to reconstruct a fitness function from the microscopic parameters.

After the description of the model (Section 2), we derive a new TSS in the limit of rare mutations (Section 3), from which a limit of small mutations gives what we call the 'canonical diffusion of adaptive dynamics' (Section 4). The coefficients of this diffusion involve the first-order derivatives of the fixation probabilities, which are computed in Section 5 as a linear combination of four fundamental components associated to fertility, defence, aggressiveness and isolation. New numerical results on the robustness of the population with respect to these fundamental components are also given, as well as some consequences on the canonical diffusion of adaptive dynamics in large populations.

The microscopic model

We will restrict here to logistic interaction. More general models are considered in [START_REF] Champagnat | Evolution of discrete populations and the canonical diffusion of adaptive dynamics[END_REF].

A monotype (binary) logistic branching process (LBP, see [START_REF] Lambert | The branching process with logistic growth[END_REF]) with dynamical parameters (b, c, d) is a Markov chain in continuous time (X t ; t ≥ 0) with nonnegative integer values and transition rates

q ij =        bi if j = i + 1 ci(i -1) if j = i -1 -i(b + c(i -1)) if j = i 0 otherwise.
The non-linear term ci(i-1) describes competition mortality due to random encounters between individuals. Other terms correspond to independent birth events with constant individual rates. This Markov chain is positiverecurrent and converges in distribution to a r.v. ξ, where ξ is a Poisson variable of parameter θ := b/c conditioned on being nonzero

P(ξ = i) = e -θ 1 -e -θ θ i i! i ≥ 1. ( 2
)
We consider in the sequel, a multitype asexual birth and death process with mutation, generalizing this LBP. At any time t, the population is composed of a finite number N (t) of individuals characterized by their phenotypic traits x 1 (t), . . . , x N (t) (t) belonging to a given trait space X , assumed to be a closed subset of R k . The population state at time t is represented by the counting measure on X

ν t = N (t) i=1 δ x i (t) .
The population dynamics is governed by the following parameters.

• b(x) is the rate of birth from an individual of type x. The function b is assumed to be C 2 b . • c(x, y) is the rate of death of an individual of type x due to the competition with another individual of type y. Therefore, the total death rate of an individual of type x in a population ν may be written as c(x, y)(ν(dy) -δ x (dy)). In this expression, the Dirac mass at x substracted to the measure ν means that the individual does not compete with himself. The function c is assumed to be C 2 b and bounded away from 0 on X 2 .

• γµ(x) is the probability that a birth from an individual with trait x produces a mutant individual, where µ(x) ∈ [0, 1] and where γ ∈ (0, 1] is a parameter scaling the frequence of mutations. When there is no mutation, the new individual inherits the trait of its progenitor. In Section 3, we will be interested in the limit of rare mutations (γ → 0).

• M (x, dh) is the law of the trait difference h = y -x between a mutant individual with trait y born from an individual with trait x. We assume that M (x, dh) has 0 expectation (no mutation bias), i.e. hM (x, dh) = 0, and has a density on R k which is uniformly bounded in x ∈ X by some function M (h) with finite third-order moment.

We will denote the dependence of ν t on the parameter γ with the notation ν γ t . Observe that such a population cannot go extinct because the death rate is 0 when there is only one individual in the population. Since we want to apply a limit of rare mutations while keeping the population size finite, this is necessary to prevent the population to become extinct before any mutation occur.

Let ξ(x) be a random variable whose law is the stationary distribution of a pure x-type population with no mutation (µ ≡ 0). This law is given by ( 2) where θ is replaced by θ(x) := b(x)/c(x, x).

The last notation needed concerns a population with initially only two types x and y and with no mutation. Then ν t = X t δ x + Y t δ y , where (X t , Y t : t ≥ 0) is a bivariate Markov chain. For this Markov chain, P(T < ∞) = 1, where T is the first time where either X t or Y t reach 0. We call fixation (of the mutant y) the event {X T = 0}. The probability of fixation will be denoted by u n,m (x, y)

u n,m (x, y) := P(X T = 0 | X 0 = n, Y 0 = m).
3 The trait substitution sequence in finite populations

In this section, we apply the limit of rare mutations (γ → 0) to the process ν γ , in order to describe the evolution of the population as a TSS in finite population. This limit requires to rescale time properly, as t/γ, to describe the evolution on the mutation timescale.

Theorem 3.1 Fix x ∈ X . Assume that ν γ 0 = N γ 0 δ x where sup γ∈(0,1) E((N γ 0 ) p ) < ∞ for some p > 1. Then, for any 0 < t 1 < . . . < t n , the n-tuple (ν γ t 1 /γ , . . . , ν γ tn/γ ) converges in law for the weak topology to (N t 1 δ St 1 , . . . , N tn δ St n ) where (1) (S t ; t ≥ 0) is a Markov jump process on X with initial value S 0 = x and whose jumping rates q(x, dh) from x to x + h are given by

q(x, dh) = β(x)χ(x, x + h)M (x, dh),
where

β(x) = µ(x)b(x)E(ξ(x)) = µ(x)b(x)θ(x)/(1 -e -θ(x) ) and χ(x, y) = n≥1 nP(ξ(x) = n) E(ξ(x)) u n,1 (x, y) = n≥1 e -θ(x) θ(x) n-1 (n -1)! u n,1 (x, y). (3)
(2) Conditional on (S t 1 , . . . , S tn ) = (x 1 , . . . , x n ), the N t i are independent and respectively distributed as ξ(x i ).

Therefore, in the limit of rare mutations, on the mutation timescale, the population is always monomorphic and the dominant trait of the population evolves as a jump process over the trait space, where a jump corresponds to the appearance and fixation of a mutant type. Moreover, at any time, the population size stationary (i.e. has the stationary distribution corresponding to the dominant trait of the population). The fixation rate of a mutant is governed by the function χ(x, y), which is therefore the random analogue of the traditional invasion fitness [START_REF] Metz | How should we define 'fitness' for general ecological scenarios ?[END_REF], defined as the probability of invasion of a mutant type y is a resident population of type x at equilibrium. Observe that, as usual in adaptive dynamics, the fitness landscape depends on the current state of the population. Moreover, in contrast with the classical TSS [START_REF] Metz | Adaptive Dynamics, a geometrical study of the consequences of nearly faithful reproduction[END_REF], from a given monomorphic resident population, any mutant trait has a positive probability to invade (by genetic drift). Therefore, evolution is possible in any direction of the trait space. However, a directional selection still exists, as will appear in the next section.

We refer to [START_REF] Champagnat | Evolution of discrete populations and the canonical diffusion of adaptive dynamics[END_REF] for the proof of this result. However, this convergence is natural in view of the following interpretation of each parameter. β(x) can be seen as the mean mutant production rate of a stationary x-type population (i.e. with size ξ(x)), and χ(x, y) is the probability of fixation of a single ytype mutant entering a pure x-type population with size-biased stationary size. The size bias comes from the fact that the mutant appears at a birth time in the stationary population (since the birth rate is proportional to the population size, the population size after a birth event in the stationary population is given by the size-biaised stationary population size).

The canonical diffusion of adaptive dynamics

Let us assume for simplicity that X = R k . Let σ(x) be the square root matrix of the covariance matrix of M (x, •). We also need to assume that the matrix σ(x) is a Lipschitz function of x.

In order to obtain the equivalent of the canonical equation of adaptive dynamics in a finite population, we want to apply a limit of small mutation steps (weak selection) to the TSS S. To this aim, we introduce a parameter ǫ > 0 and replace the mutation kernels M (x, •) by their image by the application h → ǫh. Time also has to be rescaled in order to obtain a non-degenerate limit. The correct time scaling is 1/ǫ 2 , which leads to the following generator for the rescaled TSS (S ǫ t ; t ≥ 0)

A ǫ ϕ(x) = 1 ǫ 2 R k (ϕ(x + ǫh) -ϕ(x))β(x)χ(x, x + ǫh)M (x, dh).
Using the assumption that the mutation kernels M (x, •) have 0 expectation, it is elementary to compute the limit of this expression as ǫ → 0 (for sufficiently regular ϕ). This limit, which takes the form of a diffusion generator, explains the following result (its full proof can be found in [START_REF] Champagnat | Evolution of discrete populations and the canonical diffusion of adaptive dynamics[END_REF]).

Theorem 4.1 If the family (S ǫ 0 ) ǫ>0 has bounded first-order moments and converges in law as ǫ → 0 to a random variable Z 0 , then the process S ǫ with initial state S ǫ 0 converges in law for the Skorohod topology on D(R + , R k ) to the diffusion process (Z t ; t ≥ 0) with initial state Z 0 unique solution to the stochastic differential equation

dZ t = β(Z t )σ 2 (Z t ) • ∇ 2 χ(Z t , Z t )dt + β(Z t )χ(Z t , Z t )σ(Z t ) • dB t ( 4
)
where ∇ 2 χ denotes the gradient w.r.t. the second variable y of χ(x, y) and B is a standard k-dimensional Brownian motion.

This result gives the equivalent of the canonical equation of adaptive dynamics (1) when the population is finite. It is no longer a deterministic ODE, but a diffusion process, in which the genetic drift remains present (in the form of a stochastic diffusion term), as a consequence of the population finiteness and of the limit of weak selection. The deterministic drift part of ( 4) is very similar to the standard canonical equation of adaptive dynamics (1), and involves in particular the gradient of the fitness function χ. The process (4) provides a diffusion model describing the evolution of the dominant trait value in a population [START_REF] Lande | Natural selection and random genetic drift in phenotypic evolution[END_REF][START_REF] Hansen | Stabilizing selection and the comparative analysis of adaptation[END_REF], grounded on a precise microscopic density-dependent modelling of the population dynamics. It also gives the precise balance between directional selection and genetic drift as a function of the individual's dynamical parameters.

Fixation probability near neutrality

The SDE (4) involves the fixation probability χ(x, x) and the fitness gradient with respect to the second variable ∇ 2 χ(x, x). In this section, we explain how these quantities can be explicitly computed.

We need to compute the derivatives of the fixation probabilities u n,m (x, y) when y is close to x. Recall that the law of the two-types LBP without mutation (X, Y ) used to define u n,m in the end of Section 2 is characterized by the birth vector B and the competition matrix y,y) .

C B = b(x) b(y) , C = c(x, x) c(x, y) c(y, x) c(
We will say that the mutant is neutral if all individuals are exchangeable, i.e. when b(y) = b(x) and c(x, y) = c(y, x) = c(y, y) = c(x, x) (this holds in particular when y = x). As will appear below, using the notation (b, c) := (b(x), c(x, x)), it is natural to focus on deviations from the neutral case expressed as

B = b1 + 0 λ , C = c1 - 0 0 δ δ + 0 α 0 α - 0 ε ε 0
In words, deviations from the neutral case are a linear combination of four fundamental selection coefficients λ, δ, α, ε, that are chosen to be positive when they confer an advantage to the mutant. It is convenient to assess these deviations to the neutral case in terms of Under neutrality, an elementary martingale argument shows that the fixation probability equals the initial mutant frequency p := m/(m + n). This implies in particular that

χ(x, x) = e -θ(x) -1 + θ(x) θ(x) 2 . ( 5
)
The following theorem unveils the dependence of u upon λ, δ, α, ε, when they slightly deviate from 0, and explains why these four selection coefficients provide a natural basis to decompose the gradient of the fixation probability.

Theorem 5.1 As a function of the multidimensional selection coefficient s = (λ, δ, α, ε), the probability u is differentiable, and in a neighbourhood of s = 0 (selective neutrality),

u = p + v ′ .s + o(s), (6)
where the selection gradient v = (v λ , v δ , v α , v ε ) can be expressed as

v ι n,m = p (1 -p) g ι n+m ι = ε, v ε n,m = p (1 -p) (1 -2p) g ε n+m
where the g's depend solely on the resident's characteristics b, c, and on the total initial population size n+m. They are called the invasibility coefficients.

The invasibility coefficients of a pure resident population are interesting to study, as they provide insights about how the fixation probability deviates from p as the selection coefficients of the mutant deviate from 0. They also provide information about the robustness of the resident population, i.e. its resistance to mutant invasions. In particular, this allows to compare the sensitivity of the invasion probability in a given monomorphic resident population with respect to the four fundamental selection coefficients. In the simplest case where mutations in the parameter space are isotropic, the biggest invasibility coefficent gives the direction of the parameter space where a mutant is more likely to invade. More generally, when there are correlations between mutations in the parameter space (either because of the phenotypic structure in the functions b(•) and c(•, •), or because the covariance matrix of the mutation steps is non-diagonal), the likeliest direction of evolution in the trait space is given by the deterministic coefficient of the canonical diffusion (4), in which the fitness gradient is given by

∇ 2 χ(x, x) = a λ (x)∇b(x) -a δ (x)∇ 1 c(x, x) + a α (x)∇ 2 c(x, x), (7)
where, for ι = λ, δ, α,

a ι (x) = e -θ(x) ∞ n=1 ng ι n+1 (x)θ(x) n-1 (n + 1) 2 (n -1)! .
It is possible to obtain explicit expressions for the invasibility coefficients g ι as series. We refer to [START_REF] Champagnat | Evolution of discrete populations and the canonical diffusion of adaptive dynamics[END_REF] for the exact expressions. In particular, these expressions yield that a ι (x) = âι (θ(x))/c(x, x) for some function âι . Moreover, they allow one to compute numerically the invasibility coefficients, and therefore the quantities âι for ι = λ, δ, α as functions of the parameter θ(x). These numerical results can be used to make simulations of the canonical diffusion of adaptive dynamics in various ecological examples. In particular, in contrast with the classical canonical equation of adaptive dynamics, the presence of a genetic drift can induce the evolutionary dynamics to drift away from evolutionary stable strategies, where the fitness gradient is zero.

When the fitness gradient admits several zeros, it can visit several basins of attraction on various timescales. This numerical study is a work in progress, that is quite delicate because the series involved in the computation of g ι are slowly converging, with first terms that grow exponentially fast with θ. We shall give here our first results. Fig. 1 shows the functions âι for ι = λ, δ, α. Several comments can be made from this figure. First, for any θ > 0, âδ (θ) > âα (θ) > âλ (θ). This means that, for equal mutation steps in the parameter space, a mutation is always more advantageous in the direction δ than in the direction α, which is itself more advantageous than in the direction λ. In other words, in a given population, a better defence capacity is more beneficial than a better aggressiveness, which is more beneficial than a better fertility.

Moreover, as θ goes to infinity, these functions have different asymptotic behaviors. âδ seems to converges to 1/2, whereas âλ (θ) and âα (θ) are both equivalent to 1/2θ (see Fig. 2). However, this does not mean that mutations are much more likely to fixate in the δ direction, because, when θ is large, b = θc is larger than c, and, in [START_REF] Lambert | The branching process with logistic growth[END_REF], a δ and a α are multiplied by ∇ 1 c and ∇ 2 c respectively, whereas a λ is multiplied by ∇b. More formally, to compute the limit of the canonical diffusion when θ goes to infinity, one can divide the competition kernel c(•, •) by a constant K in the microscopic model and then let K go to infinity in the canonical diffusion. Denoting by χ K (•, •) the fitness function obtained this way and using the asymptotic behaviors given above, one gets that lim K→+∞ ∇ 2 χ K (x, x) = 1 2b(x) (∇b(x) -θ(x)∇ 1 c(x, x)).

Moreover, by (5), χ K (x, x) converges to 0 when K → ∞. Now, as proved in [START_REF] Champagnat | A microscopic interpretation for adaptive dynamics trait substitution sequence models[END_REF], with our notation, the fitness function of the canonical equation of adaptive dynamics (1) is given by f (x, y) = b(y) -c(y, x)θ(x). Therefore, as K → +∞, the canonical diffusion converges to a deterministic ODE which is precisely the canonical equation of adaptive dynamics. This gives a new justification of this equation, and this also allows one to study the fluctuations around the canonical diffusion when K is large. In particular, this diffusion equation with small diffusion term enters the framework of Freidlin-Wentzell's theory [START_REF] Freidlin | Random Perturbations of Dynamical Systems[END_REF], which can be applied to predict the long time behaviour of the diffusion, and its chain of visit of basins of attractions when K is large. This kind of information is biologically very relevant, since it allows one to predict in which order all the evolutionary stable strategies will be visited by the population and on which timescale.

1 .

 1 fertility (λ): positive λ means increased mutant birth rate 2. defence capacity (δ): positive δ means reduced competition sensitivity of mutant individuals w.r.t. the total population size 3. aggressiveness (α): positive α means raised competition pressure exerted from any mutant individual onto the rest of the population 4. isolation (ε): positive ε means lighter cross-competition between the two different types, that would lead, if harsher, to a greater probability of exclusion of the less abundant one

Figure 1 :Figure 2 :

 12 Figure 1: The functions âλ , âδ and âα as functions of θ.