Elie Bretin
email: bretin@polytechnique.fr

Valerie Perrier
email: valerie.perrier@imag.fr

Phase field method for mean curvature flow with boundary constraints

This paper is concerned with the numerical approximation of mean curvature flow t → Ω(t) satisfying an additional inclusion-exclusion constraint Ω 1 ⊂ Ω(t) ⊂ Ω 2 . Classical phase field model to approximate these evolving interfaces consists in solving the Allen-Cahn equation with Dirichlet boundary conditions. In this work, we introduce a new phase field model, which can be viewed as an Allen Cahn equation with a penalized double well potential. We first justify this method by a Γ-convergence result and then show some numerical comparisons of these two different models.

Introduction

In the last decades, a lot of work has been devoted to the motion of interfaces, and particularly to motion by mean curvature. Applications concern image processing (denoising, segmentation), material sciences (motion of grain boundaries in alloys, crystal growth), biology (modeling of vesicles and blood cells), image denoising, image segmentation and motion of grain boundaries.

Let us introduce the general setting of mean curvature flows. Let Ω(t) ⊂ R d , 0 ≤ t ≤ T , denote the evolution by mean curvature of a smooth bounded domain Ω 0 = Ω(0) : the outward normal velocity V n at a point x ∈ ∂Ω(t) is given by

V n = κ, ( 1
)
where κ denotes the mean curvature at x, with the convention that κ is negative if the set is convex. We will consider only smooth motions, which are well-defined if T is sufficiently small [START_REF] Ambrosio | Geometric evolution problems, distance function and viscosity solutions[END_REF]. Since singularities may develop in finite time, one may need to consider the evolution in the sense of viscosity solutions [START_REF] Barles | Solutions de viscosité des équations de Hamilton-Jacobi[END_REF][START_REF] Evans | Phase transitions and generalized motion by mean curvature[END_REF].

The evolution of Ω(t) is closely related to the minimization of the following energy:

J(Ω) = ∂Ω 1 dσ.
Indeed, (1) can be viewed as a L 2 -gradient flow of this energy.

Following [START_REF] Modica | Un esempio di Γ --convergenza[END_REF][START_REF] Modica | Il limite nella Γ-convergenza di una famiglia di funzionali ellittici[END_REF], the functional J can be approximated by a Ginzburg-Landau functional :

J ǫ (u) = R d ǫ 2 |∇u| 2 + 1 ǫ W (u) dx.
where ǫ > 0 is a small parameter, and W is a double well potential with wells located at 0 and 1 (for example W (s) = 1 2 s 2 (1 -s) 2 ).

Modica and Mortola [START_REF] Modica | Un esempio di Γ --convergenza[END_REF][START_REF] Modica | Il limite nella Γ-convergenza di una famiglia di funzionali ellittici[END_REF] have shown the Γ-convergence of J ǫ to c W J in L 1 (R d ) (see also [START_REF] Bellettini | Variational approximation of functionals with curvatures and related properties[END_REF]), where

c W = 1 0 2W (s)ds. ( 2
)
The corresponding Allen-Cahn equation [START_REF] Allen | A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening[END_REF], obtained as the L 2 -gradient flow of J ǫ , reads

∂ u ∂t = ∆u - 1 ǫ 2 W ′ (u). (3)
Existence, uniqueness, and a comparison principle have been established for this equation (see for example chapters 14 and 15 in [START_REF] Ambrosio | Geometric evolution problems, distance function and viscosity solutions[END_REF]). To this equation, one usually associates the profile

q = arg min R 1 2 γ ′ 2 + W (γ) ; γ ∈ H 1 loc (R), γ(-∞) = 1, γ(+∞) = 0, γ(0) = 1 2 (4)
Remark 1. The profile q (when W is continuous) can also be obtained [START_REF] Alberti | Variational models for phase transitions, an approach via γ-convergence[END_REF] as the global decreasing solution of the following Cauchy problem q ′ (s) = -W (s), s ∈ R q(0) = 1 2 , and satisfies

R 1 2 q ′ (s) 2 + W (q(s)) = 1 0 2W (s)ds.
Then, the motion Ω(t) can be approximated by

Ω ǫ (t) = x ∈ R d ; u ǫ (x, t) ≥ 1 2 ,
where u ǫ is the solution of the Allen Cahn equation ( 3) with the initial condition

u ǫ (x, 0) = q d(x, Ω(0)) ǫ .
Here d(x, Ω) denotes the signed distance of a point x to the set Ω.

The convergence of ∂Ω ǫ (t) to ∂Ω(t) has been proved for smooth motions [START_REF] Chen | Generation and propagation of interfaces for reaction-diffusion equations[END_REF][START_REF] Bellettini | Quasi-optimal error estimates for the mean curvature flow with a forcing term[END_REF] and in the general case without fattening [START_REF] Barles | Solutions de viscosité des équations de Hamilton-Jacobi[END_REF][START_REF] Evans | Phase transitions and generalized motion by mean curvature[END_REF]. The convergence rate has been proved to behaves as O(ǫ 2 |log ǫ| 2 ).

From the practical point of view, this equation is usually solved in a box Q, with periodic boundary conditions, which allow solutions to be compute via a semi-implicit Fourier-spectral method as in the paper [START_REF] Chen | Applications of semi-implicit Fourier-spectral method to phase field equations[END_REF].

In this article, we investigate the approximation of interfaces evolving in a restricted area, which usually occurs in several physical applications. More precisely, we will consider mean curvature flow t → Ω(t) evolving as the L 2 gradient flow of the following energy :

J Ω 1 ,Ω 2 (Ω) = ∂Ω 1 dσ if Ω 1 ⊂ Ω ⊂ Ω 2 , +∞ otherwise.
Here Ω 1 and Ω 2 are two given smooth subsets of R d such that dist(∂Ω 1 , ∂Ω 2 ) > 0. These evolving interfaces will clearly satisfy the following constraint Ω 1 ⊂ Ω(t) ⊂ Ω 2 .

Ω Ω Ω [START_REF] Owen | Minimizers and gradient flows for singularly perturbed bi-stable potentials with a dirichlet condition[END_REF]. Yet, some limitations appear in this model :

• The Dirichlet boundary conditions prevent interfaces to reach boundaries ∂Ω 1 and ∂Ω 2 . This can be seen as a consequence of thickness of the interface layer which is about O(ǫ ln(ǫ)). This highlights the fact that the convergence rate of this model can not be better than O(ǫ ln(ǫ)).

• From a numerical point of view, the resolution of the Allen Cahn equation with Dirichlet boundary conditions can be performed using a finite element method. This appears less efficient and more difficult to implement in dimensions greater than 2 than the semiimplicit Fourier-spectral method of [START_REF] Chen | Applications of semi-implicit Fourier-spectral method to phase field equations[END_REF].

To overcome these limitations, we introduce in this paper a new phase field model. The main idea will be to consider the Allen-Cahn equation in the whole domain with a penalization technique to take into account the boundary constraints.

The paper is organized as follows:

In section 2, we present in detail the two phase field model. In section 3, we justify our penalized approach by a Γ-convergence result. In section 4, we compare our method to the classical Finite Element model of [START_REF] Owen | Minimizers and gradient flows for singularly perturbed bi-stable potentials with a dirichlet condition[END_REF], through numerical illustrations. These simulations will clarify the numerical convergence rate of each model.

Phase field model with boundary constraints

In this section, we will introduce our Allen Cahn model for the approximation of mean curvature flow t → Ω(t) evolving as the L 2 gradient flow of the following energy

J Ω 1 ,Ω 2 (Ω) = ∂Ω 1 dσ if Ω 1 ⊂ Ω ⊂ Ω 2 +∞ otherwise
where Ω 1 and Ω 2 are two given smooth subsets of R d satisfying dist(∂Ω 1 , ∂Ω 2 ) > 0. We begin with the description of the classical approach.

Classical model with Dirichlet boundary conditions

The classical strategy, see for instance [START_REF] Owen | Minimizers and gradient flows for singularly perturbed bi-stable potentials with a dirichlet condition[END_REF][START_REF] Bourdin | Design-dependent loads in topology optimization[END_REF], consists in introducing the function space

X Ω 1 ,Ω 2 = u ∈ H 1 (Ω 2 \ Ω 1 ) ; u |∂Ω 1 = 1 , u |∂Ω 2 = 0 ,
and a penalized Ginzburg-Landau energy of the form

Jǫ,Ω 1 ,Ω 2 (u) =    Ω 2 \Ω 1 ǫ 2 |∇u| 2 + 1 ǫ W (u) dx if u ∈ X Ω 1 ,Ω 2 +∞ otherwise.
In such framework, Chambolle and Bourdin [START_REF] Bourdin | Design-dependent loads in topology optimization[END_REF] have shown the Γ-convergence of Jǫ,

Ω 1 ,Ω 2 to c W J Ω 1 ,Ω 2 in L 1 (R d ) (c W has been introduced in (2)
). This approximation conduces to the following Allen-Cahn equation

       u t = △u -1 ǫ 2 W ′ (u), on Ω 2 \ Ω 1 u |∂Ω 1 = 1, u |∂Ω 2 = 0 u(0, x) = u 0 ∈ X Ω 1 ,Ω 2 .
A more general Γ-convergence result for the Allen-Cahn equation with Dirichlet boundary conditions can be found in [START_REF] Owen | Minimizers and gradient flows for singularly perturbed bi-stable potentials with a dirichlet condition[END_REF].

Novel approach with a penalized double well potential

Now, we describe an alternative approach to force the boundary constraints, based on a penalized double well potential. From W considered in section 1.1, we define two continuous and positive potentials W 1 and W 2 satisfying the following assumption :

(H1) W 1 (s) = W (s) for s ≥ 1/2 W 1 (s) ≥ max(W (s), λ) for s ≤ 1/2 and W 2 (s) = W (s) for s ≤ 1/2 W 2 (s) ≥ max(W (s), λ) for s ≥ 1/2, .
where λ > 0.

For α > 1, ǫ > 0, and x ∈ R d , we introduce also a penalized double well potential W ǫ,Ω 1 ,Ω 2 ,α defined by where dist(x, Ω 1 ) and dist(x, Ω c 2 ) are respectively the signed distance functions to the sets Ω 1 and Ω c 2 , and q is the profile function associated to W defined in (4).

W ǫ,Ω 1 ,Ω 2 ,α (s, x) = W 1 (s) q dist(x, Ω 1 ) ǫ α + W 2 (s) q dist(x, Ω c 2 ) ǫ α +W (s) 1 -q dist(x, Ω 1 ) ǫ α -q dist(x, Ω c 2 ) ǫ α , ( 5
)
Our modified Ginzburg-Landau energy J ǫ,Ω 1 ,Ω 2 ,α reads

J ǫ,Ω 1 ,Ω 2 ,α (u) = R d ǫ 2 |∇u| 2 + 1 ǫ W ǫ,Ω 1 ,Ω 2 ,α (u, x) dx. ( 6
)
We will prove in the next section that this energy Γ-converges to c W J Ω 1 ,Ω 2 . The associated Allen-Cahn equation reads in this context:

∂ t u(x, t) = △u(x, t) -1 ǫ 2 ∂ s W ǫ,Ω 1 ,Ω 2 ,α (u(x, t), x), for all (x, t) ∈ R d × [0, ∞( u(0, x) = u 0 (x)

Approximation result of the penalized Ginzburg-Landau energy

In this section, we prove the convergence of the Ginzburg-Landau energy J ǫ,Ω 1 ,Ω 2 ,α introduced in (6), to the following penalized perimeter

J Ω 1 ,Ω 2 (u) = |Du|(R d ) if u = 1l Ω and Ω 1 ⊂ Ω ⊂ Ω 2 +∞ otherwise . Remark 2. Given u ∈ L 1 (R d ), |Du|(R d ) is defined by |Du|(R d ) = sup R d u div(g)dx ; g ∈ C 1 c (R d , R d ) ,
where

C 1 c (R d ; R d ) is the set of C 1 vector functions from R d to R d with compact support in R d . If u ∈ W 1,1 (R d ), |Du| coincides with the L 1 -norm of ∇u and if u = 1l Ω where Ω has a smooth boundary, |Du| coincides with the perimeter of Ω. Moreover, u → |Du|(R d ) is lower semi-continuous in L 1 (R d ) topology.
We assume in this section that Ω 1 and Ω 2 are two given smooth subsets of R d satisfying dist(∂Ω 1 , ∂Ω 2 ) > 0, and that ǫ is sufficiently small such that

1 -q dist(x, Ω 1 ) ǫ α -q dist(x, Ω c 2 ) ǫ α > 1/2, (7)
for all x in Ω 2 \ Ω 1 .

We now state the main result for the modified Ginzburg-Landau energy J ǫ,Ω 1 ,Ω 2 ,α :

Theorem 1. Assume that W is a positive double-well potential with wells located at 0 and 1, continuous on R and such that W (s) = 0 if and only if s ∈ {0, 1}. Assume also that W 1 and W 2 are two continuous potentials satisfying assumption (H1). Then, for any α > 1, it holds

Γ -lim ǫ→0 J ǫ,Ω 1 ,Ω 2 ,α = c W J Ω 1 ,Ω 2 in L 1 (R d ).
Proof. We first prove the liminf inequality. i) Liminf inequality :

Let (u ǫ ) converge to u in L 1 (R d ). As J ǫ,Ω 1 ,Ω 2 ,α ≥ 0, it is not restrictive to assume that the lim inf of J ǫ,Ω 1 ,Ω 2 (u ǫ ) is finite. So we can extract a subsequence u n = u ǫn such that lim n→+∞ J ǫ h ,Ω 1 ,Ω 2 ,α (u n ) = lim inf ǫ→0 J ǫ,Ω 1 ,Ω 2 ,α (u ǫ ) ∈ R + .
Note that from remark [START_REF] Alberti | Variational models for phase transitions, an approach via γ-convergence[END_REF] and assumption [START_REF] Bourdin | Design-dependent loads in topology optimization[END_REF], it holds

               q dist(x,Ω 1 ) ǫ α ≥ 1/2 for x ∈ Ω 1 , q dist(x,Ω c 2 ) ǫ α ≥ 1/2 for x ∈ Ω c 2 , 1 -q dist(x,Ω 1 ) ǫ α -q dist(x,Ω c 2 ) ǫ α ≥ 1/2 for x ∈ Ω 2 \ Ω 1 , 1 -q dist(x,Ω 1 ) ǫ α -q dist(x,Ω c 2 ) ǫ α ≥ 0 for x ∈ R d ,
for ǫ sufficiently small. This implies that

Ω 1 W 1 (u n )dx ≤ Ω 1 2q dist(x, Ω 1 ) ǫ α n W 1 (u n )dx ≤ 2 R d W ǫn,Ω 1 ,Ω 2 ,α (u n , x)dx ≤ 2ǫ h J ǫn,Ω 1 ,Ω 2 ,α (u n ).
In the same way :

R d \Ω 2 W 2 (u n )dx ≤ 2ǫ n J ǫn,Ω 1 ,Ω 2 ,α (u n )
and

Ω 2 \Ω 1 W (u n )dx ≤ 2ǫ n J ǫn,Ω 1 ,Ω 2 ,α (u n ).
At the limit n → ∞, the Fatou's Lemma and the continuity of W , W 1 and W 2 imply that

Ω 1 W 1 (u)dx = 0, R d \Ω 2 W 2 (u)dx = 0 and Ω 2 \Ω 1 W (u)dx = 0. Recall also that W , W 1 and W 2 , vanish respectively at s = {0, 1}, s = {1} and s = {0}. This means that u(x) ∈        {1} a.e in Ω 1 {0} a.e in R d \ Ω 2 {0, 1} a.e in Ω 2 \ Ω 1 ,
almost everywhere. Hence, we can represent u by 1l Ω for some Borel set Ω ∈ R d satisfying Ω 1 ⊂ Ω ⊂ Ω 2 . Using the Cauchy inequality, we can estimate

J ǫn,Ω 1 ,Ω 2 ,α (u n ) ≥ R d ǫ n |∇u h | 2 2 + 1 ǫ n W (u n ) dx ( because W 1 ≥ W and W 2 ≥ W ) ≥ R d ǫ n |∇u n | 2 2 + 1 ǫ n W (u n ) dx ( where W (s) = min W (s) ; sup s∈[0,1] W (s) ) ≥ R d 2 W (u n )|∇u n |dx = R d |∇[φ(u n )]|dx = |D[φ(u n )]|(R d ), where φ(s) = s 0 2 W (t)dt. Since φ is a Lipschitz function (because W is bounded), φ(u ǫ ) converges in L 1 (R d ) to φ(u). Using the lower semicontinuity of v → |Dv|(R d ), we obtain lim n→+∞ J ǫn,Ω 1 ,Ω 2 (u n ) ≥ lim inf n→+∞ |Dφ(u n )|(R d ) ≥ |Dφ(u)|(R d ).
The lim inf inequality is finally obtained remarking that φ

(u) = φ(1l Ω ) = c W 1l Ω = c W u.
Let us now prove the limsup inequality. ii) Limsup inequality : We first assume that u = 1l Ω for some bounded open set Ω satisfying Ω 1 ⊂ Ω ⊂ Ω 2 with smooth boundaries. We introduce the sequence

u ǫ (x) = q dist(x, Ω) ǫ .
and two constants c 1 and c 2 defined by

c 1 = sup s∈[0,1] {W 1 (s) -W (s)} , and c 2 = sup s∈[0,1] {W 1 (s) -W (s)} .
Note that

J ǫ,Ω 1 ,Ω 2 ,α (u ǫ ) = R d ǫ|∇u ǫ | 2 2 + 1 ǫ W (u ǫ ) dx + R d 1 ǫ q dist(x, Ω 1 ) ǫ α (W 1 (u ǫ ) -W (u ǫ dx + R d 1 ǫ q dist(x, Ω c 2 ) ǫ α (W 2 (u ǫ ) -W (u ǫ )) dx
Each of these 3 terms above is now analyzed.

1) Estimation of the first term :

I 1 ǫ = R d ǫ|∇u ǫ | 2 2 + 1 ǫ W (u ǫ ) dx.
By co-area formula, we estimate

I 1 ǫ = 1 ǫ R d q ′ (d(x, Ω)/ǫ) 2 2 + W (q(d(x, Ω)/ǫ)) dx = 1 ǫ R g(s) q ′ (s/ǫ) 2 2 + W (q(s/ǫ)) ds = R g(ǫt) q ′ (t) 2 2 + W (q(t)) dt where g(s) = |D1l {d≤s} |(R d ).
By the smoothness of ∂Ω, g(ǫt) converges to |D1l {dist(x,Ω)≤0} |(R d ) as ǫ → 0; moreover, by definition of the profile q, u ǫ converges to 1l Ω and lim sup ǫ→0

I 1 ǫ ≤ |D1l Ω |(R d ) +∞ -∞ 1 2 |q ′ (s)| 2 + W (q(s))ds .
According to remark (1), it follows that

+∞ -∞ 1 2 |q ′ (s)| 2 + W (q(s)) ds = 1 0 2W (s)ds = c W , which implies that lim sup ǫ→0 I 1 ǫ ≤ c W |D1l Ω |(R d ).
2) Estimation of the second term :

I 2 ǫ = R d 1 ǫ q dist(x, Ω 1 ) ǫ α (W 1 (u ǫ ) -W (u ǫ )) dx.
The function dist(x, Ω) is negative on Ω 1 , thus u ǫ (x) ≥ 1 2 on Ω 1 and W 1 (u ǫ (x)) = W (u ǫ (x)) for all x ∈ Ω 1 . This means that

I 2 ǫ = R d \Ω 1 1 ǫ q dist(x, Ω 1 ) ǫ α (W 1 (u ǫ ) -W (u ǫ )) dx ≤ c 1 R d \Ω 1 1 ǫ q dist(x, Ω 1 ) ǫ α dx,
where

c 1 = sup s∈[0,1] {W 1 (s) -W (s)}.
Using co-area formula, we estimate

R d \Ω 1 1 ǫ q dist(x, Ω 1 ) ǫ α dx = ∞ 0 1 ǫ g 1 (s)q s ǫ α ds = ǫ α-1 ∞ 0 g 1 (ǫ α s)q(s)ds,
where

g 1 (s) = |D1l {dist(x,Ω 1 )≤s} |(R d ).
By the smoothness of Ω 1 , g(ǫ α t) converges to |D1l dist(x,Ω 1 )≤0 |(R d ) as ǫ → 0. We then deduce that lim sup ǫ→0 I 2 ǫ = 0, as α > 1 and ∞ 0 q(s)ds is bounded.

3) Estimation of the last term :

I 3 ǫ = R d 1 ǫ q dist(x, Ω c 2 ) ǫ α (W 2 (u ǫ ) -W (u ǫ )) dx.
This estimation is similar to the second one. The function dist(x, Ω) is positive on

R d \ Ω 2 , this means u ǫ (x) ≤ 1 2 on R d \ Ω 2 and W 2 (u ǫ (x)) = W (u ǫ (x)) for all x ∈ R d \ Ω 2 .
Then, we have

I 3 ǫ = Ω 2 1 ǫ q dist(x, Ω c 2 ) ǫ α (W 2 (u ǫ ) -W (u ǫ )) dx ≤ c 2 Ω 2 1 ǫ q dist(x, Ω c 2 ) ǫ α dx,
and using co-area formula, it holds

Ω 2 1 ǫ q dist(x, Ω c 2 ) ǫ α dx = ∞ 0 1 ǫ g 2 (s)q s ǫ α ds = ǫ α-1 ∞ 0 g 2 (ǫ α s)q(s)ds, where g 2 (s) = |D1l {dist(x,Ω 2 )≤-s} |(R d ).
We deduce as before that lim sup ǫ→0

I 3 ǫ = 0.
Finally, we conclude that lim sup

ǫ→0 J ǫ,Ω 1 ,Ω 2 ,α (u ǫ ) ≤ c W |D1l Ω |(R d ).
Remark 3. This theorem is still true in the limit case α → ∞, where J ǫ,Ω 1 ,Ω 2 ,α=∞ (u) reads

J ǫ,Ω 1 ,Ω 2 ,∞ (u) = Ω 1 ǫ|∇u| 2 2 + 1 ǫ W 1 (u) dx + Ω 2 \Ω 1 ǫ|∇u| 2 2 + 1 ǫ W (u) dx + R d \Ω 2 ǫ|∇u| 2 2 + 1 ǫ W 2 (u) .

Algorithms and numerical simulations

We now compare numerically the two phase field models described previously. The first and classical model is integrated by a semi-implicit finite element method whereas our penalized Allen Cahn equation is solved by the semi-implicit Fourier spectral algorithm. In particular, we will observe that both approaches give similar solutions but,

• the algorithm used for the penalized Allen Cahn equation is more efficient and has lower complexity than the semi-implicit finite element method used for the Allen Cahn equation with Dirichlet boundary conditions.

• the convergence rate of the phase field approximation appears to behave as about O(ǫ ln(ǫ)) for the first model and as O(ǫ 2 ln(ǫ) 2 ) for our penalized version of Allen Cahn equation.

A semi-implicit finite element method for the Allen Cahn equation with Dirichlet boundary conditions

Let us give more precision about the classic semi-implicit finite element method used for the equation

u t (x, t) = ∆u(x, t) - 1 ǫ 2 W ′ (u)(x, t), on Ω 2 \ Ω 1 × [0, T ], (8)
where

u |∂Ω 1 = 1, u |∂Ω 2 = 0 and W (s) = 1 2 s 2 (1 -s) 2 .
Note that when the initial condition u 0 is chosen of the form u 0 = q (dist(Ω 0 , x)/ǫ) with Ω 0 satisfying the constraint Ω 1 ⊂ Ω 0 ⊂ Ω 2 , then we expect that the set Ω ǫ (t) defined by

Ω ǫ (t) = Ω 1 ∪ {x ∈ Ω 2 \ Ω 1 ; u(x, t) ≥ 1/2} ,
should be a good approximation to the constraint mean curvature flow t → Ω(t).

Let us introduce a triangulation mesh T h on the set Ω 2 \ Ω 1 and the discretization time step δ t . Then, we consider the approximation spaces X h,0 and X h defined by

   X h = v ∈ H 1 (Ω 2 \ Ω 1 ) ∪ C 0 (Ω 2 \ Ω 1 ) ; v |K∈T h ∈ P k (K), v |Ω 1 = 1 and v |Ω 2 = 0 X h,0 = v ∈ H 1 0 (Ω 2 \ Ω 1 ) ∪ C 1 (Ω 2 \ Ω 1 ) ; v |K∈T h ∈ P 2 (K)
where P k denotes the polynomial space of degree k. We take k = 2 in the future numerical illustrations. Then, the solution u(x, t n ) at time t n = nδ t is approximated by U h,n , defined for n > 1 as the solution on X h of

Ω 2 \Ω 1 U h,n ϕ dx+δ t Ω 2 \Ω 1 ∇U h,n ∇ϕ dx = Ω 2 \Ω 1 U h,n-1 - δ t ǫ 2 W ′ (U h,n-1 ) ϕ dx, ∀ϕ ∈ X h,0 ,
and for n = 0 by U h,0 = arg min

v∈X h v -u 0 L 2 (Ω 2 \Ω 1 ) .
This algorithm is known to be stable under the condition

δ t ≤ c W ǫ 2 ,
where

c W = sup t∈[0,1] {W ′′ (s)} -1
. More results about stability and convergence of this algorithm can be found in [START_REF] Feng | Numerical analysis of the Allen-Cahn equation and approximation for mean curvature flows[END_REF][START_REF] Brassel | Instabilité de Forme en Croissance Cristalline[END_REF].

A time-splitting Fourier spectral method for the penalized Allen-Cahn equation

We now consider the second model

u t (x, t) = △u(x, t) - 1 ǫ 2 ∂ s W ǫ,Ω 1 ,Ω 2 ,α (u)(x, t), on Q × [0, T ], (9)
with periodic boundary conditions on a given box Q, chosen sufficiently large to contain Ω 2 . In future numerical tests, we use α = 2, W (s) = 1 2 s 2 (1 -s) 2 and the potentials W 1 , W 2 are defined by

W 1 (s) = 1 2 s 2 (1 -s) 2 if s ≥ 1 2 10(s -0.5) 4 + 1/32 otherwise and W 2 (s) = 1 2 s 2 (1 -s) 2 if s ≤ 1 2
10(s -0.5) 4 + 1/32 otherwise , which clearly satisfy the assumption (H1) (see figure [START_REF] Allen | A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening[END_REF]).

The initial condition u 0 satisfies u 0 = q (dist(Ω 0 , x)/ǫ) and we will show that the set

Ω ǫ (t) = {x ∈ Q ; u(x, t) ≥ 1/2} ,
will be a good approximation of Ω(t) as ǫ tends to zero.

Our numerical scheme for solving equation ( 9) is based on a splitting method between the diffusion and reaction terms. We take advantage of the periodicity of u by integrating exactly the diffusion term in the Fourier space. More precisely, the solution u(x, t n ) at time t n = t 0 +nδ t is approximated by its truncated Fourier series :

u n P (x) = |p|∞=P c n p e 2iπp•x .
Here |p| ∞ = max 1≤i≤d |p i | and P represents the number of Fourier modes in each direction.

The step n of our algorithm writes :

• u n+1/2 P (x) = c n+1/2 p e 2iπp•x , with c n+1/2 p = c n p e -4π 2 δt |p| 2 . • u n+1 P = u n+1/2 P -δt ǫ 2 ∂ s W ǫ,Ω 2 ,Ω 1 α (u n+1/2 P
).

In practice, the first step is performed via a Fast Fourier Transform, with a computational cost O(P d ln(P )). In the simplest case of traditional Allen Cahn equation, the corresponding numerical scheme turns out to be stable under the condition

δt ≤ c W ǫ 2 ,
and the convergence of this splitting approach has been established in [START_REF] Bretin | Méthode de champ de phase et mouvement par courbure moyenne[END_REF].

Simulations and numerical convergence

We compare in this part numerical solutions obtained with these two algorithms. For each test, we took ǫ = 2 -8 and δ t = ǫ 2 . The P 2 finite element algorithm was implemented in Freefem++.

The mesh T h used in simulations are plotted in figure [START_REF] Ambrosio | Geometric evolution problems, distance function and viscosity solutions[END_REF]. The penalization method was implemented in MATLAB where we have taken P = 2 8 .

Figure 3: Mesh T h generated by Freefem++, and used in simulations plotted on figures (4) and ( 5). In both cases, ∂Ω 1 and ∂Ω 2 are respectively identified as the green and the yellow boundaries.

We first plot two situations on figures (4) and ( 5). The functions u ǫ are plotted only on the admissible set Ω 2 \ Ω 1 for the FE Dirichlet method and on all the set Q for the pseudo-spectral penalization version. We note that the solutions obtained by both methods are very similar.

In order to estimate the convergence rate of both models, we consider the case where Ω 1 and Ω 2 are two circles of radii equal to R 1 = 0.3 and R 2 = 0.4. The situation is thus very simple when the initial set Ω 0 is also a circle with radius R 0 satisfying R 1 < R 0 < R 2 . Indeed, the penalized mean curvature motion Ω(t) evolves as a circle, with radius satisfying

R(t) = max R 2 0 -2t, R 1 , that decreases until R(t) = R 1 .
The solutions of the two different models are computed for different values of ǫ with P = 2 8 , δ t = 1/P 2 and R 0 = 0.35. In both cases, the set Ω ǫ (t) appears as a circle of radius R ǫ (t). We then estimate the numerical error between R ǫ (t) and R(t). The results obtained for the first method are plotted on figure [START_REF] Bellettini | Quasi-optimal error estimates for the mean curvature flow with a forcing term[END_REF] Moreover, our approach allows us to simulate very easily and efficiently three dimensional experiments, whatever the geometry of the sets Ω 1 and Ω 2 . See for instance figure [START_REF] Brassel | Instabilité de Forme en Croissance Cristalline[END_REF] where the numerical solution of the Allen Cahn equation is plotted for different times t.

(ε) = |R(∞) -R ε (∞)| erreur(ε) ε ln(ε) ε ε 1/2
(ε) = |R(∞) -R ε (∞)| erreur(ε) ε 2 ln 2 (ε) ε 2 ε

Possible extension of the method

Another advantage of our penalization approach is that it can be easily extended for more general situations of evolving interfaces. For example, we have recently considered a mean curvature flow with an additional forcing term g and a conservation of the volume in [START_REF] Brassel | A modified phase field approximation for mean curvature flow with conservation of the volume[END_REF]. Then, the model of [START_REF] Brassel | A modified phase field approximation for mean curvature flow with conservation of the volume[END_REF] can be modified to take into account additional inclusion-exclusion constraints by simply using potential W Ω 1 ,Ω 2 instead of W in phase field equation. In this case, this leads to the following perturbed Allen-Cahn equation Two simulations obtained from this model are plotted in figure [START_REF] Brassel | A modified phase field approximation for mean curvature flow with conservation of the volume[END_REF]. We can observe that both constraints (conservation of the volume and inclusion-exclusion set) are well respected.

u t = ∆u - 1 ǫ 2 F (u), with F (u) = W ′ Ω 1 ,Ω 2 (u) -ǫg 2W Ω 1 ,Ω 2 (u) - Q W ′ Ω 1 ,Ω 2 (u) -ǫg 2W Ω 1 ,Ω 2 (u)dx Q 2W Ω 1 ,Ω 2 (u)dx 2W Ω 1 ,Ω 2 (u).

Conclusion

This paper has presented a new phase field model for the approximation of mean curvature flow with inclusion-exclusion constraints. The classical method in such situation consists to solve the Allen-Cahn equation with Dirichlet boundary conditions. Since this method appears to be not optimal while its convergence is observed with a rate about O(ǫ ln(ǫ)) only, we have introduced a new approach based on a penalized double well potential. This method was firstly motivated by a Γ-convergence result and secondly numerical tests suggesting that its convergence rate is about O(ǫ 2 ln(ǫ) 2 ). The proof of the numerical precision of the scheme is still an open problem and will be the subject of a forthcoming paper. Another advantage of our method lies in its simplicity to be implemented, since it is only based on Fourier Transforms. This simplicity allows to consider 3D Geometries, and to elaborate new strategies in more general situations, such as mean curvature flow with forcing term and conservation of volume.

Figure 2 :

 2 Figure 2: Example of potentials W , W 1 and W 2

Figure 4 :

 4 Figure 4: Numerical solutions obtained at different times t 0 = 0, t 1 = 0.022, t 2 = 0.033 and t = 0.055. The first line corresponds to the FE Dirichlet method and the second line to the penalization pseudo-spectral approach.

Figure 5 :

 5 Figure 5: Numerical solutions obtained at different times t 0 = 0, t 1 = 0.0055, t 2 = 0.0083 and t 3 = 0.016. The first line corresponds to the FE Dirichlet method and the second line to the penalization pseudo-spectral approach.

Figure 7 :

 7 Figure 7: Penalization algorithm : numerical error |R ǫ (t) -R(t)| ; Left: t → R ǫ (t) for different values of ǫ; Right: ǫ → sup t∈[0,T ] {|R(t) -R ǫ (t)|} in logarithmic scale.

Figure 8 :

 8 Figure 8: Minimal surface estimation : solution of the Allen Cahn equation at different times t. The domain Ω 1 is the union of the two red tores and Ω 2 is the box Q with N = 2 7 , ǫ = 1/N and δ t = ǫ 2

Figure 9 :

 9 Figure 9: Two numerical experiments with a forcing term (gravity force), a volume conservation and an exclusion constraint : the domain Ω 1 is empty and Ω 2 is plotted in red in each picture. We used N = 2 7 , ǫ = 1/N and δ t = ǫ 2 .

  in logarithmic scale. It clearly appears an error of O(ǫ ln(ǫ)). The same test is done for the penalization algorithm: the results are plotted on figure[START_REF] Bourdin | Design-dependent loads in topology optimization[END_REF] and we now clearly observed a convergence rate of O(ǫ 2 ln(ǫ 2 )).

	0.37									-3
								exact
								ε= 1/N
	0.36							ε= 2/N ε= 3/N	-3.2
								ε= 4/N
	0.35
										-3.4
	0.34
										-3.6
	0.33
	R(t)									-3.8
	0.31 0.32								Erreur	-4
										-4.2
	0.3
	0.29									-4.4
	0	0.005	0.01	0.015	0.02 t	0.025	0.03	0.035	0.04	-5.6 -4.6	-5.4	-5.2	-5	-4.8	ε	-4.6	-4.4	-4.2	-4	-3.8

  Figure 6: Dirichlet algorithm : numerical error |R ǫ(t) -R(t)| ; Left: t → R ǫ (t) for different values of ǫ; Right: ǫ → sup t∈[0,T ] {|R(t) -R ǫ (t)|} in logarithmic scale.

		0.36									-4.5
									exact
									ε= 1/N
		0.35							ε= 2/N ε= 3/N ε= 4/N	-5
		0.34
											-5.5
		0.33
	R(t)	0.32									-6
		0.31								Erreur	-6.5
		0.3
											-7
		0.29
		0	0.005	0.01	0.015	0.02 t	0.025	0.03	0.035	0.04	-5.6 -7.5	-5.4	-5.2	-5	-4.8	ε	-4.6	-4.4	-4.2	-4	-3.8


Jacob L Weisdorf
email: jacob.weisdorf@econ.ku.dk.

Johan Dalgaard

Oded Galor

Haim Ofek

Karl Gunnar Persson

Holger Strulik

Made for Toil: Natural Selection at the Dawn of Agriculture ♣

Keywords: Hunting, Labour Input, Malthus, Metabolism, Neolithic Revolution, Patience Capital JEL Classification Numbers: J10, J22, Q56, O10

The labour input among pre-historic foragers was normally rewarded within the same day of the effort.

For the first farmers, by contrast, labour input and its rewards could be far apart. However, the patience was worthwhile: population growth rates among early agriculturalists were up to 60 times higher than those of their foraging counterparts. It is well-known from the biological science that humans differ with respect to metabolism. This study argues that rates of metabolism well-suited for the many hours of labour input required for farming gained an evolutionary advantage with the advent of agriculture. This theory helps shedding light on the puzzles why farming was adopted despite its high labour costs, and why people of agricultural societies work more than their foraging counterparts.

Introduction

Anthropologists have long debated how much work pre-historic people had to do to achieve subsistence before the rise of agriculture. The earlier anthropological tradition assumed that hunter and gatherers lead hard lives of constant struggle to make a living. The advent of agriculture was believed to have reduced the labour input necessary to attain subsistence, which would have afforded people the extra time required to develop advanced societies.

Systematic time-budget studies done in the 1960's among contemporaneous primitive people, however, completely turned the conventional story on its head, as labour inputs in foraging communities proved to be amazingly small. A seminal time-study done by renounced anthropologist Richard B. Lee among the Dobe Bushmen, a group of hunters and gatherers located in southern Africa, showed that an average adult, despite a harsh environment, devoted only between 15 and 20 hours to work per week [START_REF] Lee | What hunters do for a living, or, how to make out on scarce resources[END_REF]. Later studies done among other groups of foragers have confirmed Lee's observation: contemporary primitive people rarely put in more than six hours of labour per day, or about two-thirds of the labour input of people living in more advanced societies [START_REF] Clark | A Farewell to Alms: A Brief Economic History of the World[END_REF]).

Lee and fellow anthropologists [START_REF] Lee | Man the Hunter[END_REF]Sahlins, 1972) became spokesmen of the view that pre-historic foragers were effectively the most leisured people in history. By contrast to the earlier anthropological tradition, the first agriculturalists were now believed to have put in more rather than less labour to attain subsistence. Accordingly, the current standpoint holds that Palaeolithic foragers possessed all the knowledge necessary to take up agriculture, but that they would not embark upon time-costly methods of food production unless there was good reason to do so. Farming, for these people, was considered a last resort [START_REF] Fernandez-Armesto | Food: A History[END_REF]. Jack R. Harlan, one of the great pioneers of historical ecology, summarizes the contemporary view as follows:

"Why farm? Why give up the 20-hour work week and the fun of hunting in order to toil in the sun? Why work harder, for food less nutritious and a supply more capricious? Why invite famine, plague, pestilence and crowded living conditions?" (Harlan, Crops and Man, 1992) This leaves a number of compelling questions: If farming required more labour input to attain subsistence than did foraging, then why was agriculture adopted in the first place? Why was it adopted by some, and not by others? And why are contemporary hunters and gatherers so indolent compared to people of more advanced societies?

This paper argues that these issues are best understood through an acceptance of the importance of heterogeneity in human metabolism. From the biological science we know that humans (like other animals) differ with respect to metabolism, and that some metabolic rates are better suited for certain environmental conditions than others. For instance, some individuals use very little energy at rest, a desirable quality in an environment where food is scarce for long periods of time, like among hunters and gatherers. A different advantage is being able to put in much labour for extended periods of time while consuming relatively little energy. Such a feature is decisive when labour input and the rewards hereof are far apart, such as in traditional agriculture. Indeed, in any environment where human physical activity is vital for obtaining foods, the efficiency with which energy is transformed into physical activity (i.e. human metabolism) plays a crucial role for survival.

If environmental conditions remain fixed for prolonged periods of time, then the composition of types of people in the economy-as is predicted by the Darwinian model of natural selection-will change in favour of those whose metabolism afford them the highest reproduction rate. And when environmental conditions change, on the other hand, so does likely the evolutionary dominant metabolism.

For millions of years of human evolution, hunting and gathering was the only food procurement strategy practised. For pre-historic foragers, the reward from putting in labour was normally obtained within the same day of the labour input. In a world entirely dominated by hunting and gathering, therefore, patience played little or no role for human existence.

However, the origins of agriculture nearly 10,000 years ago provided a major technological and environmental break with the past. By contrast to the lifestyle of foragers, patience became an essential ingredient in the productive life of early farmers. Not only does the cultivation of foods demand a labour input not required when relying on the natural stock; the labour input is also often due several months before the effort pays off [START_REF] Gill | Seasonality and Agriculture in the Developing World[END_REF]. Nonetheless, the reward to patience among early agriculturalist was staggering: their population growth rates and population densities are estimated to have been up to 60 times higher than those of contemporaneous groups of foragers [START_REF] Seabright | The Company of Strangers: A Natural History of Economic Life[END_REF]).

There has been much debate in the archaeological and anthropological literature whether the additional labour input required to take up agriculture meant that output per unit of labour input went up or down with the advent of farming [START_REF] Harlan | Crops and Man[END_REF]. Meanwhile, the high population densities observed among early agriculturalists are taken as unmistakable evidence that output per worker in farming was significantly higher than that of their foraging counterparts. This suggests that, if people differ in terms of metabolism, and if some metabolic rates are better suited for agriculture than others, then those that were well-suited for farming would have obtained an evolutionary advantage with the advent of agriculture. Not only would this help explaining why farming would be adopted despite its high labour costs. It also indicates why agriculturalists work more than hunters and gatherers-the farmers are simply 'made' for toil.

The current research links to a growing economic literature that addresses the shift from hunting to agriculture in the Neolithic period. This literature includes Hibbs and [START_REF] Olsson | Biogeography and long-run economic development[END_REF], [START_REF] Locay | From hunting and gathering to agriculture[END_REF], [START_REF] Marceau | On the Early Holocene: Foraging to Early Agriculture[END_REF], [START_REF] North | The first economic revolution[END_REF] and [START_REF] Smith | The primitive hunter culture, Pleistocene extinction, and the rise of agriculture[END_REF]. 1 In these studies, technological change is typically the underlying driving force behind the transition, a view that is perfectly consistent with the ideas proposed in the present paper. By contrast to the current study, however, the existing literature largely abstracts from the fact that the origins of farming entailed a process of work-intensification, a key element to the puzzle about agricultural adoption (e.g. [START_REF] Harlan | Crops and Man[END_REF]).

Providing a theory of natural selection, the present research also relates to a branch of literature that considers the evolution in human genetic traits over the course of history. This literature includes [START_REF] Bowles | Endogenous Preferences: The Cultural Consequences of Markets and other Economic Institutions[END_REF], [START_REF] Clark | Survival of the Richest: The Malthusian Mechanism in Pre-Industrial England[END_REF], Galor andMoav (2002, 2007), [START_REF] Hansson | Malthusian Selection of Preferences[END_REF] and Robson andKaplan (2003, 2006), all of which entertain the idea that (pre-)historic changes in economic environments (like the emergence of markets and other economic institutions) have played a crucial role in shaping human values, taste and personalities.

The present work is perhaps most closely related to that of Galor andMoav (2002, 2007).

Their main finding-that nature eventually selected human traits complementary to the process of economic growth-is identical to that of the present paper. However, the basic mechanisms driving the results are fundamentally different. In Galor andMoav (2002, 2007), a parental tradeoff between the quantity and quality of children is key for understanding the emergence of traits suitable for economic growth. In the current context, by contrast, this role is given to heterogeneity in human metabolism.

Finally, the present paper draws parallels to a recent working paper by [START_REF] Doepke | Occupational Choice and the Spirit of Capitalism[END_REF]. By contrast to the literature above, Doepke and Zilibotti address the formation of cultural, i.e. non-genetic, preference. They stress the importance of the acquisition of so-called patience capital, a cultural practise common among the middle classes of pre-industrial England, whereby immediate gratification is delayed in exchange for the achievement of skills.

According to Doepke and Zilibotti's theory, the patience capital of the bourgeoisie, by contrast to the land of the aristocracy, was highly rewarded during the Industrial Revolution, which in turn explains why fortunes would ultimately shift in favour of the bourgeoisie. The analogy to the present work comes from the idea that patience capital can be thought of as labour invested in the cultivation of foodstuff, an alternative to relying on the natural stock of foods.

The current paper continues as follows. Section 2 offers a simple framework for addressing the issues about agricultural adoption outline above. Based on this framework, section 3 performs a graphical analysis of some of the implications related to the advent of agriculture.

Finally, section 4 concludes.

The Model

Consider an economy in which a single type of good-food (or more precisely calories)-is produced. There are two potential sectors for (or methods of) production: foraging and farming.

The differences between the two are discussed in detail below. The economy consists of L individuals, who are identical from all aspects, except for rates of metabolism. Which sector (or method) an individual decides to use depends on its metabolism, as well as the earning profiles that it faces in the two sectors.

As was likely the case for all pre-industrial societies (e.g. Clark 2007), the economy under consideration is assumed to be subject to so-called Malthusian population dynamics. Such dynamics involve the combination of two elements. The first element is the law of diminishing returns in the production of output. In the presence of a fixed factor, such as land, the application of the law of diminishing returns implies that any improvement in productive potential is ultimately swallowed up by a larger population. The second element, which links directly to Malthus (1798), is a positive correlation between the resources of parents and their reproductive success. 2 The combination of these elements means that, in the long run, the economy is kept in a homeostatic, so-called Malthusian, equilibrium, where the population level is constant, and where incomes are maintained at level of subsistence.

By contrast to the Malthusian model normally used, 3 subsistence in the present context is not tied to a specific income level. Here, instead, subsistence refers to the caloric surplus (defined as an individual's income net of its own-consumption) that permits an individual exactly one surviving offspring. In this sense, the current framework provides an alternative way of thinking about subsistence in a Malthusian equilibrium. 4 Below, Malthusian population dynamics appear through the application of three basic relationships. The first is a positive relationship between an individual's labour input and its potential earnings. The second is a positive relationship between an individual's labour input and its caloric requirements. And the third is a positive relationship between an individual's caloric surplus (caloric output minus caloric requirements) and its reproductive success.

Labour Input and Income Level

The first relationship-that between labour input and income earned-is captured in the following way. Suppose that the real income (measured in units of calories) earned by an individual can be written as wi = wi (l,li,Li,Bi). The w function represents the individual's earning profile in sector i∈[H,A] (H for hunting, A for agriculture). The variable l measures the individual's labour input, which is to be found endogenously below. The variable li (l but bold) is a vector of the labour input of all other workers in sector i; Li is the total number of workers employed in sector i; and Bi measures the level of total factor productivity in sector i. 5 2 According to Malthus, the relationship between the two would be regulated either through positive (i.e. mortality) checks, or though preventive (i.e. fertility) checks. In his original work, Malthus (1798) speculated that marital, and thus birth, rates, would respond to change in the price of provisions (mostly that of wheat). However, the conclusions that arise from the so-called modern Malthusian model, which is used in the present context, are qualitatively identical. 3 See [START_REF] Clark | A Farewell to Alms: A Brief Economic History of the World[END_REF] for a detailed exposition of the modern Malthusian model and the Malthusian law of population. 4 See [START_REF] Dalgaard | Subsistence -A Bio-economic Foundation of the Malthusian Equilibrium[END_REF] for a similarly approach. 5 For simplicity, we abstract from the use of capital goods in production.

The function w is assumed to be continues in all its arguments. Diminishing returns to labour (i.e. to hours worked as well as workers employed) applies under the assumptions that (wl,wll,-wl,-wL,-wlL)>0. It is also assumed that (wB,wlB)>0, so that a higher level of total factor productivity in sector i increases the level of income as well as the value of an additional hour's work. Finally, it is assumed that wi(0,l,Li,Bi)=0, so that no work yields no pay.

As was discussed above, foraging and farming are distinguishable from at least two perspectives. Firstly, farming, unlike foraging, demands a labour investment-or requires a certain amount of patience capital-before output in this sector can be reaped. Secondly, because of the labour investment required for farming, a necessary (but not sufficient) condition for agriculture to emerge is that the earning profile for farming is steeper than that for foraging.

To capture these features, suppose that wA(l,l,LA,BA)=0 for 0≤l≤λ, and that ∂wA>/∂l>∂wH/∂l for l>λ, where λ>0 is the size of the labour investment-or patience capital-needed per farmer before any earnings in this sector can be obtained. Note also that, since foraging and farming are subject to seasonal variation [START_REF] Gill | Seasonality and Agriculture in the Developing World[END_REF], labour input in the model is measured per annum.

Specifically, this implies that λ is the total number of hours per year per farmer required in agriculture before foodstuffs can be reaped.

Labour Input and Caloric Requirements

The second relationship-that between labour input and caloric requirements-builds on observations made in the field of human biology. Following [START_REF] Leslie | Caloric requirements of human populations[END_REF], individual caloric requirements above the so-called resting metabolic rate (RMR) appear to increase with physical activities (see Table 1).

[Table 1 about here]

To capture this relationship in a non-complicated manner, suppose that the total caloric requirements of an individual, who puts in l hours of identical (i.e. physically equally demanding) labour, can be written as v = v(l). The function v is assumed to be continues and monotonic, with v(0)=RMR>0. 6 It is also assumed that vll≤0, meaning that energy-requirements increase at least at the same rate as labour input.

Caloric Surplus and Reproductive Success

The third and finally relationship-that between an individual's caloric surplus and its reproductive success-is captured in the following manner. Suppose that the reproductive success of an individual is checked (in a Malthusian sense) by the size of its caloric surplus (i.e. its caloric output minus its caloric requirements). Symbolically, the reproductive success, i.e. the number of surviving children, of an individual is simply given by the function n=n(s), where s=w-v is the individual's caloric surplus, and where n is assumed to be continues and monotonic, with n(0)=0 and n(∞)>1. Together, these assumptions imply that a so-called subsistence caloric surplus exists, defined as the caloric surplus at which an individual as able to raise exactly one surviving offspring.

Agent Heterogeneity

It has long been a well-established fact in the biological science that human metabolism is determined by a number of factors. Some of these are personally controlled (at least to some extent) by features such as weight, nutritional intake and lifestyle. Others depend on genetically determined traits, such as height, sex and race [START_REF] Benedict | Race: A Factor in Human Metabolism[END_REF].

Below, all individual are assumed to be identical, except with respect to their genetically determined rate of metabolism. 7 For the purpose of clarity, moreover, suppose that only two types of individuals (i.e. rates of metabolism) potentially exist. One type has a relative low resting metabolic rate (RMR), i.e. demands a fairly low amount of energy at rest. At the same time, 6 To simplify matters, the assumption that v(0) equals RMR implies that non-work related (i.e. leisurely) activities do not require any energy-use. In other words, people in the model spend their spare-time at rest. A change of this assumption will complicate matters without affecting the qualitative results of the model. 7 Allowing people to differ also with respect to non-genetic factors will not affect the qualitative results of the model. however, much work for this type requires relatively large amounts of energy. Therefore, this type is referred to throughout as a work-inefficient type (type IE in short).

The other type has a relatively high resting metabolic rate, i.e. requires comparatively many calories while resting. But, by contrast to the work-inefficient type, this so-called work-efficient type (type E in short) is capable of supplying much labour using comparatively modest amounts of energy. 8 Expressed in symbolic terms, the difference between the two types of individuals imply that RMRIE=vIE(0)<RMRE=vE( 0), and that ∂vIE>/∂l>∂vE/∂l>0 for l>0. Figure 1 provides an illustration of the two types of people and their respective metabolism (cf. the vE curve and vIE curve in the graph).

[Figure 1 about here]

Transmission of Traits

Suppose that children in the model are born to a single parent, and that the metabolism of a parent is passed on to its offspring according to the following rule. With probability σ1, the offspring obtains the same trait as its parent. However, with probability σ0, the offspring obtains the trait of the opposite type. By this construction, even if only type IE people exist in the economy, then there is always a tiny possibility that a type E person (e.g. by mutation) will emerge, and vice versa.

8 From an empirical point of view, it is highly likely that the work-efficient types coincide with a type of humans who were comparatively well-suited for tolerating, i.e. digesting and derive energy from, wheat and other sorts of grains-especially those containing gluten which is believed to have consisted of a major problem for human digestion and nutritional intake (e.g. [START_REF] Cordain | The Paleolithic Diet[END_REF]. The diet of hunters and gatherers mainly consist of roots and lean meats, both of which are excellent for providing fibres and protein; this kind of a diet dampens the appetite, but is not particularly energy-dense. By contrast, the various sorts of grains produced by farmers all contain high levels of carbohydrates, i.e. sugars and starches, and, compared to the diet of foragers, therefore, are highly energydense. What is more, the shift from a reliance on roots and meats to a reliance on grains also helps explaining the deterioration in the health observed among early farmers [START_REF] Cohen | Paleopathology at the Origins of Agriculture[END_REF].

Optimisation and Labour Input

The final element of the model concerns the optimal labour input. Faced with information about its metabolism (captured by the v function) and its earning potential in the two sectors (captured by the w functions), an individual obtains its maximal caloric surplus when a number of hours is allocated to work, so that the marginal revenue of the labour input equals the marginal cost, i.e. so that wl = vl. In Figure 2 below, the symbol l* indicates the optimal labour input any type of metabolism.

Analysis

In the following, the framework presented above is used to perform a graphical analysis of some of the implications related to the advent of agriculture. The main aim is to try to shed light on the questions formulated in the introduction section. Namely, why was farming ultimately chosen over foraging when it needed more labour input to attain subsistence? Why was it chosen by some, and not by others? And why are contemporary hunters and gatherers so indolent compared to people of more advanced societies?

Before the Origins of Agriculture

Figure 2 provides an illustration of the optimal labour input of the two potential types of individuals-type IE and type E-before the rise of farming. Note that the existence of an earning profile for farming (the wA-curve in Figure 2) suggests that agricultural is theoretical possible, but that, at the moment, it is economically unviable for both types. That farming is economically unviable means simply that it demands more units of energy to use than it is able to provide.

Nonetheless, the presence of an earning profile for farming is consistent with the belief among anthropologists and archaeologists that knowledge about how to practise agriculture can exist without being used (e.g. [START_REF] Harlan | The Living Fields: Our Agricultural Heritage[END_REF].

[Figure 2 about here] It is also evident from Figure 2 that the work-inefficient type (type IE)-because its metabolism is better suited for foraging than that of type E-is capable of generating a larger caloric surplus (cf. the vertical distance between the w-curve and the v-curve) than its workefficient counterpart. Type IE, therefore, has a greater rate of reproduction than type E. Hence, as long as the earning profiles for foraging and farming remain fixed, the composition of types of people in the economy will gradually shift in favour of type IE individuals. This process continues until type E people ultimately disappear, either entirely (if their reproduction rate is below that of replacement) or as a share of the population (if their reproduction rate is above or equal to replacement).

Malthusian Equilibrium

Suppose that the economy depicted in Figure 2 is in a Malthusian equilibrium. In this equilibrium, the caloric surplus of type IE (the evolutionary optimal type under the current conditions) is at subsistence level. In figure 2, this is indicated by the bold double-arrow, stretching between the income level of a type IE individual and its caloric use. Since the metabolism of type E people is less well-suited for foraging, their caloric surplus is below the level of subsistence, indicated by the dotted double-arrow in Figure 2. In the Malthusian equilibrium portrayed in Figure 2, therefore, only type IE individuals exist, and an occasional occurrence of type E individuals, for reasons explained above, will only be temporary.

How does a positive shock to productivity in the foraging sector affect the economy? As is shown in Figure 3, an increase in total factor productivity in hunting (from BH to BH'), by shifting up the wH-curve, in the short run will increase the labour input and the caloric surplus of both types. However, as is evident from illustration, it will not affect the evolutionary optimal trait, which is still that of type IE. That is, even though type E individuals may (in principle) raise their caloric surpluses above the level of subsistence after the shock to productivity, type IE people are still generating the larger caloric surplus of the two.

[Figure 3 about here] Since, starting from a Malthusian equilibrium, an increase in the caloric surplus improves the reproductive success of individuals, a positive shock to productivity leads eventually to growth in the size of the population. In the long run, therefore, and because of diminishing returns to labour, more people and more hours worked drive down earnings (i.e. shift down the wH-curve), until the economy is back in a Malthusian equilibrium. This is illustrated in Figure 4. Hence, prior to the origins of agriculture, the metabolism of type IE people remains the dominant trait in the population.

[Figure 4 about here]

The Origins of Agriculture What factors that ultimately made farming economically viable is not an important issue in the present context. For this matter, the reader is referred to the existing theories about agricultural adoption (see Weisdorf 2005 for a survey). The important thing is that farming-at least for some types of metabolism-was eventually capable of providing more energy than was demanded for its use. In other words, the earning profile for farming would need to shift upward and/or to the left for agriculture to be practised.

There are two analytically interesting scenarios for agricultural adoption. Both of these are examined in the following. In the first scenario, only one type of metabolism (that of type E people) finds it advantageous to make the transition to farming. In the second scenario, both types will to shift, at least in the short run. For each scenario, the analysis will demonstrate the effects, in the short as well as in the long run, of agricultural adoption on labour input, population growth, the dominant trait, and the composition of traits in the two sectors and in the economy as a whole.

Note that the origins of agriculture-both in the short and the long run-entail a process of work-intensification. That is, the labour input-at least for some-has increased. What is more, some people have taken up agriculture while others remain in foraging. In the short run, there is population growth among farmers; foragers, by contrast, face a constant population level both over the short and the long run. Finally, regardless of the time-horizon, farmers, as has been pointed out by archaeologists and anthropologists [START_REF] Harlan | The Living Fields: Our Agricultural Heritage[END_REF], put in more labour to attain subsistence than their foraging counterparts.

Agricultural Adoption: Scenario 2

In the second scenario worth analyzing-the one in which both types find it advantageous to shift to farming-the economy ends in a Malthusian equilibrium qualitatively identical to that of the previous scenario. To see how, suppose that the earning profile for farming shifts further up and/or more to the left than was previously the case. With Figure 2 as a starting point, Figure 7 provides an illustration. As before, the work-efficient type (type E) shifts from a caloric surplus below subsistence (dotted double-arrow) to one which is above (dashed double-arrow).

Accordingly, under the new conditions, and as the illustration shows, type E people are able now to increase in number.

[Figure 7 about here] By contrast to the previous scenario, however, type IE people in this case will also make the transition to agriculture. For the purpose of clarity, the wA-curve in the current scenario shifts just enough, so that type IE individuals are indifferent between using foraging and farming. That is, the caloric surpluses of type IE individuals afford them a constant population level (cf. the bold double-arrows), regardless of the sector in which they choose to be employed. 11 As in the previous scenario, therefore, type E people now carry the evolutionary dominant trait, at least in the short run.

Since, in the short run, the labour force in agriculture is increasing (as type E people increase in number), it follows that the earning profile for farming (i.e. the wA-curve) will gradually shift downward. In the long run, therefore, which is illustrated in Figure 8, the caloric surplus of both types of individuals is depressed. Hence, as the caloric surplus of type IE individuals eventually (in the illustrated case: immediately) falls below the level of replacement, this type is better off returning to foraging.

[Figure 8 about here] Like in the previous scenario, therefore, the Malthusian equilibrium of this scenario has only type E people in farming, as type IE people ultimately resort to foraging. As above, the caloric surplus of both types is at subsistence level, so that population growth is absent altogether.

Accordingly, the Malthusian equilibrium of the present scenario, which is illustrated in Figure 8, is qualitatively identical to that of the previous scenario, which was illustrated in Figure 6.

Note that agricultural adoption also in the second scenario leads to a process of workintensification (i.e. the labour input needed to attain subsistence-at least for some-has increased). In the long run, some are employed in agriculture while others remain in foraging. In the short run, there is population growth, but only among farmers. Finally, as in the previous scenario, foragers put in less effort to attain subsistence than do agriculturalists.

Henceforth, whenever productivity growth emerges, shifting up the earning profile in the sector in which it occurs, it permits a higher population density for the people employed in that sector. Meanwhile, since, from a historical perspective, productivity growth has appeared more frequently in farming than in foraging, the composition of types of individuals in the economy has gradually shifted in favour of the work-efficient type, i.e. people who were 'made' for toil.

Conclusion

The advent of agriculture has made possible an enormous increase in the population of humans.

Over a period of ten millennia, the world's population level has risen from about six millions by the time of the first transition to farming, to some six billions by the present day. Apart from a few thousand people still depending on foraging for existence, the majority of the earth's population now rely (nearly) entirely on agriculture for subsistence.

This study demonstrates that people, whose rate of metabolism was well-suited farming, would gain an evolutionary advantage over less work-efficient people once agriculture became economically viable. Since most present-day people on the earth descend from pre-historic farmers, the current theory predicts that today's world is populated by people, who are genetically made for toil. It will prove the theory wrong, however, if the metabolism among people of modern agricultural societies do not differ significantly from those of contemporary foragers. Weisdorf, J. ( 2005

Table 1 . Rates of Energy Expenditure at Model Activity Levels

 1 ), "From Foraging to Farming: Explaining the Neolithic Revolution", Journal of Economic Surveys 19:4, pp. 561-586.

	Level	Cost above RMR	Example of activity
		(kcal/kg/hr)
	Resting	0.00	Lying still
	Sitting	0.25	Sitting quietly
	Standing	0.50	Standing "at ease"
	Walking	2.50	4.8 km/hr, no load
	Heavy	4.50	Hoeing, ploughing
	Very heavy	7.00	Tree felling
	Extreme	10.00	Near maximum
			exertion
	Source: Leslie et al. (1984, Table 1)

A Malthusian equilibrium before the rise of agriculture

	Caloric needs (v)	> subsistence = subsistence < subsistence		v IE (l)
	Income
	(w)
		s(l IE *)	s(l E *)	w H (l,l H ,L H ,B H )	v E (l)
				w A (l,l A ,L A ,B A )
	RMR E
	RMR IE
			λ	labour input
				(l)
	Figure 3: A

positive shock to productivity before the rise of agriculture (short-run effect)

  (l,l H ,L H ,B H ) v E (l) w A (l,l A ,L A ,B A ') w A (l,l A ,L A ',B A ') A (l,l A ,L A ,B A ') w A (l,l A ,L A ',B A ') w H (l,l H ,L H ,B H ) v E (l)

	v IE (l) w H (l,l H ,L H ,B H ) w H (l,l H ,L H ,B H ') v E (l) w A (l,l A ,L A ,B A ) labour input v IE (l) w A (l,l A ,L A ,B A ) s(l E *) labour input RMR v IE (l) RMR IE RMR E Caloric needs (v) Income (w) l IE * l E * λ = subsistence > subsistence < subsistence s(l IE *) s(l E *) l E * l IE * RMR IE RMR E Caloric needs (v) Income (w) = subsistence > subsistence < subsistence s(l IE *) l IE * RMR IE l E * Caloric needs (v) Income (w) λ' w H (l,l H ,L H ,B H ) v E (l) λ w A (l,l A ,L A ,B A ) w A (l,l A ,L A ,B A ') = subsistence > subsistence < subsistence s(l E *) s(l IE *) l E * l IE * labour input RMR IE RMR v IE (l) Caloric needs (v) Income (w) l IE * labour input RMR IE l E * λ' λ l IE * labour input RMR IE l E * s(l IE *) λ' w H λ' = subsistence > subsistence < subsistence s(l E *) Income (w) s(l IE *) w A (l,l A ,L A ,B A ') w H (l,l H ,L H ,B H ) v E (l) w A (l,l A ,L A ,B A ) RMR E s(l E *) Income (w) RMR v IE (l) Caloric needs (v) = subsistence > subsistence < subsistence Caloric needs (v) v IE (l) > subsistence = subsistence < subsistence
	(l) (l) (l) (l) (l)

Figure 4: A positive shock to productivity before the rise of agriculture (long-run effect)

λ w H (l,l H ,L H ,B H ') w H (l,l H ,L H ',B H ') v E (l) E

Figure 5: The rise of agriculture when only type E shifts (short-run effect) w H (l,l,L,B H ) E Figure 6: The rise of agriculture when only type E shifts (long-run effect) E Figure 7: The rise of agriculture when type IE and E both shift (short-run effect) Figure 8: The rise of agriculture when type IE and E both shift (long-run effect) w

Weisdorf (2005) provides a survey of most of these contributions.

Effectively, of course, for the type IE people to actually shift to agriculture, their caloric surplus in farming should be (at least slightly) higher than that in foraging. However, for the arguments made below, this is not important.

Agricultural Adoption: Scenario 1

The first scenario is illustrated in Figure 5. The point of departure is the Malthusian equilibrium demonstrated in Figure 2. In this scenario, the earning profile for farming shifts upward and to the left, but just enough so that the work-efficient type (type E) is willing to take up agriculture. 9 As is evident form the illustration, type E finds it advantageous to shift to farming, because its metabolism enables it to obtain a larger caloric surplus when using farming compared to foraging.

As the illustration shows, this is clearly not the case for type IE people.

[Figure 5 about here]

As is also demonstrated in Figure 5, the caloric surplus of type E individuals now exceeds the level of subsistence (the dashed double-arrow compared to the bold one), meaning that type E people now increase in number. 10 Type IE people, on the other hand, do not. Hence, the metabolism of type E individuals, at least in the short run, have replaced that of type IE people as the dominant trait in the economy.

What will happen in the long run? The first thing to note is that type E individuals, according to the rule for transmitting traits, give birth predominantly to type E individuals. Hence, since the reproductive rate of type E people exceeds that of replacement, the labour force of the farming sector increases over time. Due to diminishing returns to labour, therefore, the earning profile for farming (i.e. the wA-curve) gradually shifts downward, until the caloric surpluses of type E individuals eventually reach subsistence. Thus, in the long-run, as illustrated in Figure 6, the population level of both types is stagnant (none of the two types of metabolisms are dominant in an evolutionary sense), and the caloric surplus for both types is at subsistence level. Hence, the economy has returned to a Malthusian equilibrium, but, this time, one in which foraging and farming co-exists.

[Figure 6 about here] 9 Note that, if only type IE people exist at this point in time, then farming will not emerge (even though it is indeed economically viable for type E) before type E individuals emerge themselves. 10 If this was not the case, i.e. if type E individuals were incapable of producing a caloric surplus sufficient to keep a constant population, then they would ultimately disappear from the population, and the farming sector would subsequently disintegrate.